Assessment of Spatio-Temporal Landscape Changes from VHR Images in Three Di ﬀ erent Permafrost Areas in the Western Russian Arctic

: Our study highlights the usefulness of very high resolution (VHR) images to detect various types of disturbances over permafrost areas using three example regions in different permafrost zones. The study focuses on detecting subtle changes in land cover classes, thermokarst water bodies, river dynamics, retrogressive thaw slumps (RTS) and infrastructure in the Yamal Peninsula, Urengoy and Pechora regions. Very high-resolution optical imagery (sub-meter) derived from WorldView, QuickBird and GeoEye in conjunction with declassified Corona images were involved in the analyses. The comparison of very high-resolution images acquired in 2003 / 2004 and 2016 / 2017 indicates a pronounced increase in the extent of tundra and a slight increase of land covered by water. The number of water bodies increased in all three regions, especially in discontinuous permafrost, where 14.86% of new lakes and ponds were initiated between 2003 and 2017. The analysis of the evolution of two river channels in Yamal and Urengoy indicates the dominance of erosion during the last two decades. An increase of both rivers’ lengths and a significant widening of the river channels were also observed. The number and total surface of RTS in the Yamal Peninsula strongly increased between 2004 and 2016. A mean annual headwall retreat rate of 1.86 m / year was calculated. Extensive networks of infrastructure occurred in the Yamal Peninsula in the last two decades, stimulating the initiation of new thermokarst features. The significant warming and seasonal variations of the hydrologic cycle, in particular, increased snow water equivalent acted in favor of deepening of the active layer; thus, an increasing number of thermokarst lake formations.


Introduction
Over the past three decades, arctic regions have experienced various changes in land cover [1] and landscape structure [2] due to the rapid increase of air temperature [3] and human development [4]. In high-latitude regions, the recent atmospheric warming trend was considerably stronger than elsewhere, with an average increase of surface air temperature of about 0.6 • C/decade over the last 30 years [3]. Arctic regions are extensively underlain by permafrost and thus very sensitive to rising The URG site is located in the Urengoy region (67°45′N, 76°90′E), about 150 km north of Novy Urengoy city in the discontinuous permafrost zone within a tributary of the Pur River basin. This site has a size of 105.09 km 2 , with a maximum altitude of 30-35 m developed on a marine terrace (Drozdov et al. 2015). URG has a MAAT of −7.4 °C and a MAP of 510 mm. At this site, the permafrost thickness generally reaches 100-300 m [40]. The landscape is dominated by low shrub tundra and contains The YAM site is located in the western part of the Yamal Peninsula (70 • 40 N, 68 • 48 E), 40 km from the coast of the Kara Sea, on the Bovanenkovo gas field and it covers an area of 122.67 km 2 . YAM lies within the continuous permafrost zone and based on ERA-Interim data has a mean annual air temperature (MAAT) of −8.3 • C and mean annual precipitation (MAP) of 444 mm . This area has a maximum elevation of 94 m and is dominated by alluvial-lacustrine-marine plains and terraces dissected by rivers (Seyyakha river crosses the area from NE to SW). Several hillslope processes, such as gully erosion [37] and landslides [38], characterize this region. The landscape is dominated by shrub tundra, thermokarst lakes and ponds, as well as poorly vegetated areas affected by landslides (RTS) mostly near the Seyyakha river or lake shores. In this continuous permafrost zone, the permafrost thickness reaches at least 500 m on the marine and coastal-marine plains and reduces to 100-150 m at the younger river terraces [39].
The URG site is located in the Urengoy region (67 • 45 N, 76 • 90 E), about 150 km north of Novy Urengoy city in the discontinuous permafrost zone within a tributary of the Pur River basin. This site has a size of 105.09 km 2 , with a maximum altitude of 30-35 m developed on a marine terrace (Drozdov et al. 2015). URG has a MAAT of −7.4 • C and a MAP of 510 mm. At this site, the permafrost thickness generally reaches 100-300 m [40]. The landscape is dominated by low shrub tundra and contains extensive peatlands [41] and high lake densities [42]. Isolated larch trees also occur, especially near the Pur River tributary crossing the URG site from west to east.
The PEC site is located on the eastern side of the Pechora delta (68 • 12 N, 54 • 30 E), at the south of Bolvansky Cape station. It covers an area of 78.43 km 2 and lies in the vicinity of the Pechora river in a sporadic permafrost zone. The topography is very flat with elevation up to 40 m and low shrub tundra is the dominant land cover. MAAT at this site is −3.1 • C, whereas MAP is 550 mm. This rolling coastal plain belongs to the southern tundra subzone. The surface lithology consists of Quaternary loams and sands with inclusions of pebbles [43].

Data Sources
VHR multi-sensor remotely sensed imagery from different years courtesy of DigitalGlobe (http: //www.digitalglobe.com), along with high-resolution Corona images, were used to assess landscape changes in the study sites. The selected images have cloud cover of less than 5% and were acquired during the peak vegetation growing season from July to August [44]. All the images used in this study are listed in Table 1. Along with satellite images from 2003/2004 and 2016/2017, in two of the three analyzed sites, Corona images were used for the assessment of river dynamics. The Corona images are declassified, Cold War era satellite surveillance system images taken from 1963 to 1980, these being panchromatic photographs with large swath widths and varying ground resolutions (0.6-15 m). These images were geo-referenced and co-registered with the newest Quickbird imagery for each site to allow for change detection.
Moreover, for the analysis of retrogressive thaw slumps from the YAM site, a very high-resolution digital elevation model (ArcticDEM) at 2 m spatial resolution (https://www.pgc.umn.edu/data/arcticd em/) and derived parameters were used for the delineation of these features.

Climate Data and Active Layer Thickness
As in the study of Chet , an et al. [45], climate conditions at geolocations of the three case study areas are inferred from ERA-Interim atmospheric reanalysis [36]. Previous studies [45,46] showed good agreement between climate data from meteorological stations operated by the Russian hydrometeorological observation network (Russian Institute of Hydrometeorological Information-World Data Center, RIHMI-WDC, available at http://meteo.ru/) and ERA-Interim data in our region of interest. This also confirms the results of Lindsay et al. (2014) [47] indicating that ERA-Interim is well suited for high northern latitude applications. Therefore, for this study, we use ERA-Interim data rather than nearby hydrometeorological station data.
Atmospheric reanalysis is a technique of assimilation of meteorological observation data using numerical circulation models. Atmospheric reanalysis intends to remove eventual physical inconsistencies due to deficiency of meteorological observations (e.g., sparse distribution of meteorological stations or shortcomings of observational techniques) and it is considered the most accurate reconstruction of past climate conditions. The horizontal resolution of ERA-Interim is~75 km.
Following results from Georgievski et al. [46], we have complemented our analysis with the active layer thickness from ESA Permafrost climate change initiative (CCI). The map of ESA CCI active layer thickness [35] is based on the transient permafrost model CryoGrid CCI driven by land surface temperature (MODIS LST/ ESA LST CCI) and land cover (ESA LC CCI). ESA Permafrost CCI maps are available for period 2003-2017, and the spatial resolution is~1 km. Considering the small-scale heterogeneity medium resolution, ESA Permafrost CCI active layer thickness has been verified against in situ observations documented in the CALM database. The CALM network is a coordinated attempt to standardize measurements of active layer properties in circumpolar regions. Comparison of ESA Permafrost CCI active layer thickness with corresponding CALM sites measurements shows that datasets are either in good agreement (Vaskiny Dachi) or qualitatively good at capturing trends (Urengoy Gas Field GP15 and Bolvansky Cape).
Since CALM data are not available for all the areas covered by VHR imagery, in further analysis we use ESA Permafrost CCI active layer thickness spatially averaged for YAM, URG and PEC regions.
We note that the resolution of the ERA-Interim data is rather coarse compared to the scale of VHR images and may not fully represent the exact conditions at the three sites. However, on the one hand, they can be used to identify large-scale climate drivers of environmental change that also act at small scale and, hence, may contribute to permafrost degradation. On the other hand, it has been shown that ERA-Interim data generally agree well with observations in the high northern latitudes [45,47,48].

Land Use Land Cover (LULC) Changes
For the analysis of land cover change, a supervised classification of each satellite image was conducted through an object-based image analysis (OBIA) approach. OBIA is now a widely accepted and used alternative, demonstrated as superior to pixel-based methods [49][50][51]. For comparison purposes, each site was classified into 6 classes: water, shrub tundra, grassland and sparse vegetation, disturbed inundated areas and barren including artificial surfaces.
The first step of the supervised classification consists of delineating training and validation samples. The polygons (samples) have been delineated for the analyzed years, within each study area. Each polygon was manually delineated on the true color satellite image, aided by the panchromatic band, covering a unique land cover class. The areas were chosen such that the land-cover class could be clearly identified on the image. To independently assess the accuracy of land cover maps, the total number of samples was randomly split into training samples (70%) and validation samples (30%), before the classification procedure.
The OBIA land cover classification was conducted using eCognition Developer software, version 8.9, starting with the image segmentation process and resulting in automated delineated polygons (objects) by merging contiguous pixels from the four spectral bands. Subsequently, the resulting objects have been classified into land cover classes using the Random Forest (RF) machine learning method [52]. The RF method has proved to be among the most accurate and widely used methods for land cover classification, since it is robust, has capabilities to work with various data types and it is not affected by statistical data distribution [53,54]. Mean and standard deviation of the pixel values within training objects, for each of the four spectral bands, were used as input features for the RF classification process with 500 decision trees. The RF classification resulted in a land cover map for each year and study area, with the five pre-defined classes.
The accuracy of each land cover map was independently assessed based on the validation sample dataset, regarding the widely used metrics of overall accuracy (OA), Cohen's kappa index of agreement (KIA), producer's accuracy (PA) as omission error and user's accuracy (UA) as commission error [55].

Water Bodies and Fluvial Dynamics
Based on the results from the landcover classification using object-based image analysis and random forest, connected pixels of class water were aggregated into objects. From all resulting water bodies, the rivers were manually excluded, and the remaining objects were exported as polygons and classified as lakes and ponds based on a threshold areas of 100 m 2 [25]. The generated polygons were then used to assess changes in lake area and lake abundance in each case study region during the observation period. Moreover, the lakes/ponds were classified based on the magnitude of changes into disappeared, new, relatively stable (up to 5% areal extent changes), shrinking and growing for each site.
The planform changes of two different river channels in YAM and URG were analyzed. We used four scenes (1961, 1976, 2004 and 2016) to detect subtle changes in the Seyyakha river channel and three (1967, 2003 and 2017) for the tributary of Pur. Digital methods that parse the data derived from remote sources (e.g., satellite images, orthophotos) are very suitable for detecting the changes that occur in river channels. To quantify the rate of spatio-temporal changes in watercourses, vector polygons were generated for each analyzed year.
For each polygon, the centerline that follows the configuration of the minor riverbed was extracted using the Polygon to centerline toolbox in ArcGIS 10.4.1. This helps to evaluate the total length of the river and, at the same time, is the base file for the channel investigations.
To assess channel shifts and the rates of lateral migration, we used the Channel Migration Toolbox [56]. Using the Transect tool that automatically generates perpendicular lines on the centerline of the river course, channel width at every 100 m was calculated. The mean values of flow direction correspond to the azimuth of transect lines generated by Digital Shoreline Analysis System (DSAS) version 5.0 [57]. Likewise, to establish the channel type, the sinuosity index was calculated using the Sinuosity tool by dividing the length of the centerline by the straight line distance between the end points of the channel.
Based on each channel polygon, the fingerprint of the main geomorphic processes was calculated. The total amount of land surface erosion and deposition was determined as the difference between the pairs of channels for the available imagery (1961-1976, 1976-2004, 2004-2016, for the Seyyakha river sector from Yamal;and 1967-2017, for the tributary of the Pur River).

Landslides and Infrastructure
We examined the occurrence of RTS in the YAM site on the multispectral images from 2004, 2016 and available digital elevation model (DEM) and manually digitized all the features. Both polygon datasets were counted and their corresponding area calculated. The headwall of each landslide was manually digitized in both scenes and then the longitudinal change in each case was measured in ArcGIS. Longitudinal transects were drawn at every 1 m between the 2004 and 2016 headwalls and the mean annual value of retreat was calculated for each feature.
The infrastructure components were also manually digitized on two different datasets (2004 and 2016) in the Yamal Peninsula using ArcGIS 10.4.1. Roads and buildings were digitized as polygons and pipelines as lines in both datasets.

Climate Drivers of Environmental Change
Annual mean 2 m air temperature increased in the period from 2003 to 2017 at the geolocations of all three sites ( Figure 2a). The linear trend analysis using the least-squares method shows the most substantial increase at PEC (~0.6 • C/decade), a slightly weaker increase at YAM (~0.5 • C/decade) and the lowest increase is detected at URG (~0.3 • C/decade). Georgievski et al. (in preparation) analyzed climate patterns relevant for permafrost degradation and identified abrupt warming events (years with annual mean temperature warmer than the 5-year running mean temperature) have a high impact on deepening of active layer thickness. The comparison of Figure 2a,b shows that this is also a valid assumption for the three sites; as the warm year events coincide with the strongest deepening of the active layer. Furthermore, the corresponding linear trends show the strongest deepening at the PEC site (~51.5 cm/decade), while at the URG and YAM sites, active layer thicknesses increased bỹ 22.7 cm/decade and~15.3 cm/decade, respectively. Although the YAM site is warming at a faster rate than URG (~0.5 • C/decade vs.~0.3 • C/decade) during the considered period, the latter site is characterized by a quicker deepening of the active layer (~22.7 cm/decade vs.~15.3 cm/decade). This probably reflects the impact of changes in the hydrological cycle. Though the annual mean precipitation shows a notable declining trend during the VHR observation period (2003-2017) at all sites (YAM: −3.2 cm/decade, URG: −1.3 cm/decade, PEC: −6.9 cm/decade), there is no significant change in the precipitation trend during the whole ERA-Interim period (1979-2018) except for a slight increase at the PEC site. However, seasonal variations in precipitation, in particular, the increasing trend during the transition months between winter and summer (March, April and May), may lead to changes in snow cover that can have a critical impact on the active layer thickness. A trend analysis of the snow water equivalent (SWE) shows almost the same increase in the annual mean for both the YAM and URG sites (~3 cm/decade), where at the latter site it shows a much stronger increasing linear trend in the month with maximum SWE (~11 cm/decade vs.~5 cm/decade). Therefore, it is very likely that the active layer at the URG site does not completely freeze during the months with maximum SWE as the deeper snow cover leads to better isolation of the soil from the cold surface than at the YAM site. In a warming climate, this can compensate for a slower warming of the air temperature and, hence, result in a stronger trend in the deepening of the active layer. The PEC site shows a fragile increasing trend for both annual SWE (~0.8 cm/decade) and maximum SWE (1.5 cm/decade). Therefore, the primary driver of active layer deepening at the PEC site during the considered period (2003-2017) is the increase of mean annual air temperature, while at the YAM and URG sites changes in hydrologic cycle seem to modify the impact of warming on the deepening of active layer thickness. Trends of increasing SWE at our sites support findings of Bulygina et al. [48] for the previous period based on snow depth observations. They also found increasing trends of maximum snow depth in the Western Russian Arctic.

Land Use Land Cover (LULC) Changes
The comparison of major land cover fractions between the VHR scenes shows a general decrease of grassland and an increase of shrubs and water classes in all sites ( Figure 3). Here, the increase in water classes is most evident at the discontinuous permafrost URG site.

Land Use Land Cover (LULC) Changes
The comparison of major land cover fractions between the VHR scenes shows a general decrease of grassland and an increase of shrubs and water classes in all sites ( Figure 3). Here, the increase in water classes is most evident at the discontinuous permafrost URG site. Overall, the generated land cover maps recorded excellent accuracy values with OA ranging from 0.91 to 0.95, while KIA ranges between 0.87 and 0.93 ( Table 2). As expected, the high accuracy values (close to 1) were achieved for the water class due to their notable spectral separability. Shrub tundra, including sparse forest in the Urengoy site, also yielded very high accuracies ranging from 0. 85    Overall, the generated land cover maps recorded excellent accuracy values with OA ranging from 0.91 to 0.95, while KIA ranges between 0.87 and 0.93 ( Table 2). As expected, the high accuracy values (close to 1) were achieved for the water class due to their notable spectral separability. Shrub tundra, including sparse forest in the Urengoy site, also yielded very high accuracies ranging from 0. 85   For the Pechora site, the shrub tundra class showed the largest dynamics between 2004 and 2016, when it extended its area by 24.4% (4 km 2 ) homogeneously across the study area. Here, the replaced classes mostly comprise grassland/sparse vegetation (72.5%, 2.9 km 2 ) and disturbed, inundated areas (26%, 1.04 km 2 ). On the contrary, grassland/sparse vegetation, and barren, including artificial surfaces, areas decreased by 10.3% (3.4 km 2 ), concentrated in the south-east and western part, and 10.5% (0.8 km 2 ) concentrated in the north-east, respectively. Grassland/sparse vegetation, besides shrub tundra, lost area due to the expansion of disturbed, inundated areas (1.1 km 2 ), while barren class, was replaced by grassland of 0.51 km 2 and disturbed, inundated areas of 0.31 km 2 . Water surface area increased in the analyzed interval by 1.7% (0.23 km 2 ), with expansions mainly localized around small lakes and the northwestern shore. The disturbed, inundated areas class showed the lowest net change with only −0.04 km 2 ( Figure 4).
The same tendencies have been observed for the Yamal study area, except barren and artificial surfaces, which recorded a positive net change of 14.5% (0.91 km 2 ). They mainly replaced shrub tundra (on 1.28 km 2 ) and grassland (0.07 km 2 ). Here, new roads were identified on the map in the northern and western part of the study area, accompanied by new buildings occurring during the analyzed interval. A net gain of 3.86% (2.25 km 2 ) was observed for shrub tundra, mostly replacing grasslands (2.42 km 2 ) and disturbed, inundated areas (1.12 km 2 ), predominantly in the northern and western part of the study area. However, part of this gain was compensated for by a loss of 1.28 km 2 to artificial surfaces. The surface water area increased by 5.17% (0.95 km 2 ) through lake expansion replacing disturbed, inundated areas (0.56 km 2 ), barren, including artificial surfaces (0.29 km 2 ) and grasslands (0.1 km 2 ). Grasslands lost 15.88% (4.08 km 2 ) due to shrub development of 2.42 km 2 and the expansion of disturbed, inundated areas of 1.49 km 2 , localized mainly in the south-western part of the study area ( Figure 5). Similar to Pechora, the disturbed, inundated areas showed the lowest net change (−0.04 km 2 ).
For Urengoy, positive net changes have been observed for shrub tundra and sparse forest (2.94%) and water bodies (10.53%), negative net changes for grassland (9.83%) and disturbed, inundated areas (9.83%). In contrast, the barren area (including artificial surfaces) did not change during the analyzed interval. Shrub tundra and sparse forest developed mainly at the expense of disturbed, inundated areas (1.13 km 2 ), grassland (1.03 km 2 ) and barren land (0.17 km 2 ), with new patches predominantly appearing in the south-eastern part. Water bodies mainly replaced disturbed, inundated areas (0.4 km 2 ), shrub tundra and sparse forest (0.51 km 2 ) and grassland (0.08 km 2 ), mainly due to changes in the river watercourse and appearance of new lakes. Grasslands and disturbed, inundated areas primarily developed into shrubs, with the latter also changing to pristine lakes ( Figure 6).
Remote Sens. 2020, 12, x 12 of 28 inundated areas primarily developed into shrubs, with the latter also changing to pristine lakes ( Figure 6).

Water Bodies (Lake/Ponds) Dynamics
The total area covered by lakes and ponds larger than 0.01 ha slightly increased in each case study region. The total water bodies surface area increased by 0.45% and 4.72% in PEC and YAM, respectively, whereas in URG, the area covered by water bodies larger than 100 m 2 expanded by 2.38% during the observation period (Figure 7a). Half of the lakes and ponds remained predominantly stable in PEC, whereas in URG and YAM the majority of the water bodies show an increasing trend. In URG, 317 new lakes and ponds were initiated between 2003 and 2017, whereas in PEC, only 20 newly-formed water bodies appeared during the 12-years interval (Figure 7c).

Water Bodies (Lake/Ponds) Dynamics
The total area covered by lakes and ponds larger than 0.01 ha slightly increased in each case study region. The total water bodies surface area increased by 0.45% and 4.72% in PEC and YAM, respectively, whereas in URG, the area covered by water bodies larger than 100 m 2 expanded by 2.38% during the observation period (Figure 7a). Half of the lakes and ponds remained predominantly stable in PEC, whereas in URG and YAM the majority of the water bodies show an increasing trend. In URG, 317 new lakes and ponds were initiated between 2003 and 2017, whereas in PEC, only 20 newly-formed water bodies appeared during the 12-years interval (Figure 7c). The number of lakes and ponds larger than 0.01 ha increased by 0.74% in PEC, 4.46% in YAM and 14.86% in the URG region during the observation period. Ponds are the dominant type of water bodies in each study site (e.g., YAM-73%, URG-88.5% and PEC-73.2%). Both lakes and ponds increased in number in URG and YAM, whereas in PEC, only ponds experienced an increase, whereas lakes number remained stable (Figure 7b).
The mean lake area changed slightly in PEC (−0.28%) and YAM (+0.1%), and decreased considerably in URG (−10.75%). Despite the lake density in URG (10.2/km 2 ) is twice as large as in PEC (5.2/km 2 ) and YAM (3.6/km 2 ), the average size of inventoried lakes is smaller in URG (8 ha vs. 40 and 20 ha, respectively). The increase in the number of ponds in URG is a result of thermokarst processes but is also related to the partial drainage of larger lakes.
Although the majority of water bodies are growing, the fraction of change per individual lake (%) is more significant for shrinking lakes (Figure 7d). By far, the most remarkable changes over the analyzed interval occurred in discontinuous permafrost area (Figure 8). The number of lakes and ponds larger than 0.01 ha increased by 0.74% in PEC, 4.46% in YAM and 14.86% in the URG region during the observation period. Ponds are the dominant type of water bodies in each study site (e.g., YAM-73%, URG-88.5% and PEC-73.2%). Both lakes and ponds increased in number in URG and YAM, whereas in PEC, only ponds experienced an increase, whereas lakes number remained stable (Figure 7b).
The mean lake area changed slightly in PEC (−0.28%) and YAM (+0.1%), and decreased considerably in URG (−10.75%). Despite the lake density in URG (10.2/km 2 ) is twice as large as in PEC (5.2/km 2 ) and YAM (3.6/km 2 ), the average size of inventoried lakes is smaller in URG (8 ha vs. 40 and 20 ha, respectively). The increase in the number of ponds in URG is a result of thermokarst processes but is also related to the partial drainage of larger lakes.
Although the majority of water bodies are growing, the fraction of change per individual lake (%) is more significant for shrinking lakes (Figure 7d). By far, the most remarkable changes over the analyzed interval occurred in discontinuous permafrost area (Figure 8).

Fluvial Dynamics
For at least six months per year, snow and ice cover the drainage networks in these regions, thereby decreasing the intensity of in-channel geomorphological processes [58]. The erosion is dominant in both analyzed cases, except from 1961 to 1976 for the Seyyakha river (YAM) (Figure 9a). The annual erosion rates have greater values for the tributary of the Pur river (URG). Still, between 2004 and 2016 the erosion was barely twice as high as the accumulation ratio for the Seyyakha river.
The length of the Seyyakha river increased by 0.65% in the first period, 1.71% in the second period and 0.42% in the third. The length of the river in URG decreased in the first phase (3.4%) but increased (1.79%) between 2003 and 2017 after new meander bends were initiated (Figure 9d). The sinuosity index of the Seyyakha river increased in all periods as the meanders continuously developed (Figure 9f). In URG, the initial decrease of the sinuosity was followed by a slight increase in the last years.
Channel migration is a natural phenomenon as rivers try to reach a dynamic equilibrium state. The migration rates for the Seyyakha river indicate a mature stage and a quasi-stable channel, with annual migration values between 1.

Fluvial Dynamics
For at least six months per year, snow and ice cover the drainage networks in these regions, thereby decreasing the intensity of in-channel geomorphological processes [58]. The erosion is dominant in both analyzed cases, except from 1961 to 1976 for the Seyyakha river (YAM) (Figure 9a). The annual erosion rates have greater values for the tributary of the Pur river (URG). Still, between 2004 and 2016 the erosion was barely twice as high as the accumulation ratio for the Seyyakha river.
The length of the Seyyakha river increased by 0.65% in the first period, 1.71% in the second period and 0.42% in the third. The length of the river in URG decreased in the first phase (3.4%) but increased (1.79%) between 2003 and 2017 after new meander bends were initiated (Figure 9d). The sinuosity index of the Seyyakha river increased in all periods as the meanders continuously developed (Figure 9f). In URG, the initial decrease of the sinuosity was followed by a slight increase in the last years.
Channel migration is a natural phenomenon as rivers try to reach a dynamic equilibrium state. The migration rates for the Seyyakha river indicate a mature stage and a quasi-stable channel, with annual migration values between 1.8 m (1961-1976)

Retrogressive Thaw Slumps
In the Yamal Peninsula, the majority of retrogressive thaw slumps (RTS) were identified in sloped terrain along lakeshores and the Seyyakha River ( Figure 10).
RTS generally occur in clusters where favorable site-specific conditions for their development exist. For Yamal, RTS are abundant in the western part and were initiated by thermo-erosion and mechanical erosion due to fluvial processes and waves. In a few cases, the development of RTS is related to the melting of ice-wedge polygons, known as "baydjarakhs" [59]. The features cover small spatial footprints ranging from 0.01 to 1.32 ha (mean area: 0.26 ha). The majority of the features have lengths/widths lower than 100 m and in some cases present evidence of coalescence. Furthermore, 51% of the RTS have a length to width ratio below 1. The number of RTS increased from 24 in 2004 to 37 in 2016 whereas the total area of thaw slumps had grown by about 69% for the same interval. The mean rate of slump headwall retreat ranged between 0.3 m/year and 4.9 m/year uring the period from 2004 to 2016. Based on the VHR satellite images from 2004 and 2016, the mean rate of RTS growth was between 11 and 679.8 m 2 /y ( Figure 11).

Retrogressive Thaw Slumps
In the Yamal Peninsula, the majority of retrogressive thaw slumps (RTS) were identified in sloped terrain along lakeshores and the Seyyakha River ( Figure 10).
RTS generally occur in clusters where favorable site-specific conditions for their development exist. For Yamal, RTS are abundant in the western part and were initiated by thermo-erosion and mechanical erosion due to fluvial processes and waves. In a few cases, the development of RTS is related to the melting of ice-wedge polygons, known as "baydjarakhs" [59]. The features cover small spatial footprints ranging from 0.01 to 1.32 ha (mean area: 0.26 ha). The majority of the features have lengths/widths lower than 100 m and in some cases present evidence of coalescence. Furthermore, 51% of the RTS have a length to width ratio below 1. The number of RTS increased from 24 in 2004 to 37 in 2016 whereas the total area of thaw slumps had grown by about 69% for the same interval. The mean rate of slump headwall retreat ranged between 0.3 m/year and 4.9 m/year uring the period from 2004 to 2016. Based on the VHR satellite images from 2004 and 2016, the mean rate of RTS growth was between 11 and 679.8 m 2 /y ( Figure 11).

Infrastructure
In the Yamal Peninsula, enormous unexploited gas resources were discovered in the 1980s. Since then, an extensive network of roads, drilling sites, buildings and gas pipelines appeared around the significant gas deposits. In the investigated site within the Bovanenkovo gas field, the total area covered by infrastructure constructions (roads and buildings) increased from 107.08 to 308.02 ha between 2004 and 2016 ( Figure 12). The surface covered by roads expanded by 3.7 times, whereas the total area of buildings increased by 2.2 times during the 12-years interval. Moreover, the entire length of gas pipelines increased by 7.4 times in the same period. The expansion of the infrastructure in YAM was related to notable changes in hydrological features ( Figure 13). In YAM, no footprint of roads, buildings or pipelines is found in the 1961 and 1976 Corona images. In PEC and URG, the buildings are sparse and concrete roads are entirely missing.
Remote Sens. 2020, 12, x 19 of 28 In the Yamal Peninsula, enormous unexploited gas resources were discovered in the 1980s. Since then, an extensive network of roads, drilling sites, buildings and gas pipelines appeared around the significant gas deposits. In the investigated site within the Bovanenkovo gas field, the total area covered by infrastructure constructions (roads and buildings) increased from 107.08 to 308.02 ha between 2004 and 2016 ( Figure 12). The surface covered by roads expanded by 3.7 times, whereas the total area of buildings increased by 2.2 times during the 12-years interval. Moreover, the entire length of gas pipelines increased by 7.4 times in the same period. The expansion of the infrastructure in YAM was related to notable changes in hydrological features ( Figure 13). In YAM, no footprint of roads, buildings or pipelines is found in the 1961 and 1976 Corona images. In PEC and URG, the buildings are sparse and concrete roads are entirely missing.

Discussion
VHR sensors (e.g., QuickBird, WorldView, GeoEye) have the potential to reveal subtle landscape changes in the permafrost dominated landscape and to overcome the limitations of medium and lowresolution satellite images [60]. Although the latter satellite images are preferred for extensive mapping of broad-scale vegetation or water bodies, these datasets are inappropriate to capture local high spatial heterogeneity of land use classes in Arctic regions [61]. The use of multispectral VHR data proved to be ideal for detecting different categories of land cover and for assessing a wide variety of disturbances. These images are also essential to enhance the fuzzy boundaries delineation of various types of land cover (e.g., water, sparse vegetation, shrubs).

Regional Climate Impact on Permafrost Degradation
Identifying climate changes on the regional scale that drive the local environmental change in the high northern latitudes remains a challenge [62]. Here we explore the role of regional climate drivers obtained from ERA-Interim (~75 km horizontal resolution) on the processes that are hypothesized to affect processes of the ecosystem and environmental change on the local scale. The active layer deepening indicated by data from ESA Permafrost CCI (~1 km horizontal resolution) serves as a link between the low-resolution climate data and the VHR remote sensing data. In the considered period at the three study sites, we find that warming is the primary driver of deepening of the active layer, while changes in the hydrological cycle can modify the impact of warming. Warm air and ground temperatures, changes in precipitation and deep snow have been recognized as major climatic drivers affecting permafrost degradation [63], for example, expressed by deepening of the active layer, RTS and thermokarst lake formation. At our three study sites, warming and increased SWE have been identified to act in favor of deepening of the active layer, which consequently results in an increasing number of RTS and thermokarst lake formation. Changes in vegetation can further

Discussion
VHR sensors (e.g., QuickBird, WorldView, GeoEye) have the potential to reveal subtle landscape changes in the permafrost dominated landscape and to overcome the limitations of medium and low-resolution satellite images [60]. Although the latter satellite images are preferred for extensive mapping of broad-scale vegetation or water bodies, these datasets are inappropriate to capture local high spatial heterogeneity of land use classes in Arctic regions [61]. The use of multispectral VHR data proved to be ideal for detecting different categories of land cover and for assessing a wide variety of disturbances. These images are also essential to enhance the fuzzy boundaries delineation of various types of land cover (e.g., water, sparse vegetation, shrubs).

Regional Climate Impact on Permafrost Degradation
Identifying climate changes on the regional scale that drive the local environmental change in the high northern latitudes remains a challenge [62]. Here we explore the role of regional climate drivers obtained from ERA-Interim (~75 km horizontal resolution) on the processes that are hypothesized to affect processes of the ecosystem and environmental change on the local scale. The active layer deepening indicated by data from ESA Permafrost CCI (~1 km horizontal resolution) serves as a link between the low-resolution climate data and the VHR remote sensing data. In the considered period at the three study sites, we find that warming is the primary driver of deepening of the active layer, while changes in the hydrological cycle can modify the impact of warming. Warm air and ground temperatures, changes in precipitation and deep snow have been recognized as major climatic drivers affecting permafrost degradation [63], for example, expressed by deepening of the active layer, RTS and thermokarst lake formation. At our three study sites, warming and increased SWE have been identified to act in favor of deepening of the active layer, which consequently results in an increasing number of RTS and thermokarst lake formation. Changes in vegetation can further amplify the impact of increasing SWE. For all three regions, an increase in shrub tundra is documented. This type of vegetation is prone to trap snow, which keeps the ground insulated so that warmer soil temperature can be maintained during the cold season.

LULC Changes and River Dynamics
Our analysis revealed accelerated land surface dynamics across three different permafrost zones in western Siberia and northeastern European Russia. A pronounced increase in the extent of tundra shrubs was found in all cases. This is in agreement with previous findings, which considered this type of Arctic disturbance as a critical driver of tundra greening [64][65][66][67][68][69]. In the three investigated Arctic study sites, recent expansions in shrublands were associated with the loss of the grassland and sparse vegetation class. Based on experimental and observational research, previous studies revealed that climate warming is the driving factor of the so-called Arctic "shrubification" [70]. Compelling pieces of evidence for the latitudinal shrub lines advancing in response to warming atmosphere temperatures were presented in different studies [71][72][73][74]. Furthermore, the increase in productivity of Arctic tundra is frequently related to recent summer warming [45,75,76].
The transition from land to water is another spatially extensive change taking place in each of the three study sites. Water area gain is subject to thermokarst activity, which is responsible for new pond formation, lake basin expansions and river widening in YAM and URG. In discontinuous permafrost (URG) the growth of water surfaces was twice as large as in continuous permafrost (YAM), whereas, in sporadic permafrost (PEC), we noticed only a slight increase. In all three study sites, ponds (below 1 ha) are the dominant water body type, as elsewhere in the Arctic lowlands [21,25].
Unfortunately, the dynamics of Arctic ponds have received comparatively little attention due to the limited availability of VHR images. The existing global inventories of lakes derived from Landsat data [77] have large inaccuracies in permafrost environments [20]. By using VHR images (spatial resolution < 5 m), ponds as small as 0.01 ha can be accurately delineated [25,78]. A direct comparison with the PeRL database [25] revealed good correspondence between our findings and similar analyses using VHR data. Thus, for the YAM site, we mapped 445 water bodies, whereas, in PeRL database, 443 features larger than 0.01 ha were inventoried [25]. The increase in number and surface area of lakes and ponds in the three Arctic sites agrees with other similar analyses in high latitude environments [79,80].
The variation of channel characteristics relates to transport and deposition of load amount, and with specific fluvial and thermal erosion [81]. Our results indicate the dominance of erosion during the last two decades, whereas in-channel accumulation characterized the 1961-1976 period in the Seyyakha River. Both considered rivers revealed a tendency of a slight increase in their lengths between 2003 and 2017. The imbalance between erosion and accumulation also led to changes in width, sinuosity and river flow direction [81]. Because erosion of banks is almost twice as large as the accumulation at Seyyakha, a significant widening of the river was observed in the last decades. Here, the fast retreat of ice-rich river banks is initiated by a combination of mechanical and thermal erosion [82]. An increase of the fluvial thermal erosion due to rising water temperature during the flood season was observed in the last decades [83]. This increase in water temperature might explain the recent widening of Seyyakha.
In the case of Seyyakha, the annual lateral migration of the channel is similar to other results reported in the Mackenzie Delta [81,84]. The river channel mobility in the Urengoy is twice as high, mainly due to the rapid evolution of small meandering channels. However, since the annual rates of fluvial channel changes are subtle in both YAM and URG, the available coarse-resolution products are inappropriate to resolve local scale river-banks variation.

RTS Landslide
Because thaw slumps in the Arctic regions are commonly small-sized features, they are below the detection limit of medium resolution satellite images (e.g., Landsat) [20]. VHR satellite images allowed identifying thaw slumps as small as 10 m in diameter within the Bovanenkovo field. The morphology and distribution of identified RTSs do not differ from similar features described in Central Yakutia [18], Banks Island [85], Yamal Peninsula [86] or Ellesmere Island [87]. We have identified small features in the vicinity of Seyyakha River with areas generally lower than 1 ha. Only two RTSs reach a maximum size of 13,000 m 2 . In the Yamal Peninsula, the RTS are preferentially localized along the lake/river shores. The abundance of these features in ice-rich conditions along rivers and lake shores was described in many other regions in the Arctic [18,88,89]. Our analysis revealed an evident increase in the occurrence of RTSs between 2004 and 2016. Similar findings were reported by other recent studies in the Arctic regions [85,87,90,91].
Retrogressive thaw slumps and "baydjarakhs" are generally considered as indicators of permafrost degradation in continuous permafrost zones [18,58,92]. Once the scarp of the RTS is created by ice-rich permafrost thawing, a further thaw is expected as permafrost is more exposed in this way [58]. Thus, the thaw slump expands by retrogressive thermo-erosion with rates usually ranging from centimeters to tens of meters per year [93]. The analysis of 22 RTSs in the Yamal Peninsula during the period from 2004 to 2016 provided a mean headwall retreat rate for all RTSs of 1.86 m/y. In comparison, Ward-Jones et al. [87] reported a mean retreat of 6.2 m/y for Ellesmere Island and Lantuit and Pollard [94] found 9.6 m/y for the coast of Herschel Island (Canada). On the other hand, our headwall retreat rates appear to be higher than the values provided by Lantz and Kokelj [88] for the Mackenzie Delta (0.5 to 1 m/y), but similar to Séjourné et al. [18] for Central Yakutia (0.5 to 3.16 m/y).

Infrastructure
An expansive industrial infrastructure occurred in the Yamal Peninsula between 2004 and 2016. In YAM, the intense human pressure is associated with hydrological disturbances and vegetation disappearance. The increased human pressure stimulated the initiation of new thermokarst features (Figure 13), such as in other ice-rich permafrost landscapes [95]. The observed land cover disturbances are usually accompanied by the modification of the ground thermal regime [63], causing substantial adverse impacts on fragile and sensitive Arctic ecosystems.
The main disadvantages of VHR satellite imagery are their high price and limited availability. However, previous studies have shown the capabilities of both VHR and medium-resolution satellite images. Thus, the surface area of water extracted from Landsat-5 TM amounted to 44-95% of the total surface area [21], whereas Sentinel-2 show promising results when mapping water bodies larger than 350 m 2 (Freitas et al., 2019). However, the delineation of small water bodies' boundaries is problematic when the pixel footprint does not lie entirely within the water body [96]. Landsat satellite images allowed the identification of only 25% of the total number of landslides in Himalaya [97] while SPOT images prevent the identification of about 30% of the small landslides in Hong Kong [98]. In the Western Russian Arctic the majority of the thaw slumps are too small to be detected by medium resolution satellite imagery, whereas the calculation of headwall retreat rate is feasible only if using VHR satellite or aerial images with sub-metric accuracies. The river's water surface delineation revealed that Planet imagery (3 m spatial resolution) allows accuracy of only 3% greater than Sentinel-2 data [99]. Because the size of individual features of infrastructure (e.g., roads and buildings) is very small specific spatial resolution imagery is required for their identification. The inventories of build-up areas based on Landsat are incomplete and lack subcategories [100]. In particular cases, Sentinel-2 allows identification of 3 m wide roads [101].

Conclusions
Our study demonstrates the usefulness of VHR satellite data for the detection of smaller features and quantifying subtle changes in Arctic environments. These data allow assessing a wider variety of changes in the permafrost environments compared with medium and low-resolution satellite images. In addition, we linked the scale of regional climate drivers inferred from global ERA-Interim reanalysis via the intermediate resolution ESA Permafrost CCI product to the scales at which the processes hypothesized to affect local vegetation and landscape dynamics operate. The three investigated sites in the Western Russian Arctic experienced a range of land cover changes. Our results indicate a pronounced increase in the extent of tundra shrubs in all cases. The proportion of land covered by water also increased significantly in the Yamal and Urengoy sites. The total area covered by lakes and ponds and the number of water bodies increased in all sites, especially in Urengoy, where 317 new lakes and ponds were initiated between 2003 and 2017. The most obvious consequence of the interplay between fluvial accumulation and erosion is the recent increase in river lengths. Relatively high annual rates of lateral channel migration were documented in Urengoy, whereas Seyyakha River experienced a significant widening of the channel. The number of retrogressive thaw slumps in the Yamal Peninsula increased by 54% over the observation period, whereas the mean annual headwall retreat rate is 1.86 m/y. The impact of human activities on the hydrology and land cover in the Yamal Peninsula has intensified in the last two decades.
Since the image resolution is a crucial factor for the detection of subtle disturbances in the Arctic there is a great need for accessible VHR data over large circumpolar areas to improve our understanding of their response to climate changes. Currently, the low availability and the high costs of the images impede more comprehensive research at larger spatial scales.