Periglacial vegetation dynamics in Arctic Russia: decadal analysis of tundra regeneration on landslides with time series satellite imagery

Changes in vegetation productivity based on normalized difference vegetation index (NDVI) have been reported from Arctic regions. Most studies use very coarse spatial resolution remote sensing data that cannot isolate landscape level factors. For example, on Yamal Peninsula in West Siberia enhanced willow growth has been linked to widespread landslide activity, but the effect of landslides on regional NDVI dynamics is unknown. Here we apply a novel satellite-based NDVI analysis to investigate the vegetation regeneration patterns of active-layer detachments following a major landslide event in 1989. We analyzed time series data of Landsat and very high-resolution (VHR) imagery from QuickBird-2 and WorldView-2 and 3 characterizing a study area of ca. 35 km2. Landsat revealed that natural regeneration of low Arctic tundra progressed rapidly during the first two decades after the landslide event. However, during the past decade, the difference between landslide shear surfaces and surrounding areas remained relatively unchanged despite the advance of vegetation succession. Time series also revealed that NDVI generally declined since 2013 within the study area. The VHR imagery allowed detection of NDVI change ‘hot-spots’ that included temporary degradation of vegetation cover, as well as new and expanding thaw slumps, which were too small to be detected from Landsat satellite data. Our study demonstrates that landslides can have pronounced and long-lasting impacts on tundra vegetation. Thermokarst landslides and associated impacts on vegetation will likely become increasingly common in NW Siberia and other Arctic regions with continued warming.


Introduction
Rapid change in Arctic vegetation productivity since the 1980's-generally referred to as 'Arctic greening'-has been linked to climate warming and to changes in sea ice extent (Bhatt et al 2010, Dutrieux et al 2012, Myers-Smith et al 2020, Berner et al 2020. The strongest responses observed have been increases in shrub cover and height, especially within the Low Arctic tundra zone, where vegetation is highly sensitive to air and soil temperature during the growing season (Epstein et al 2004, Macias-Fauria et al 2012. Increased shrub cover lowers albedo of the land surface, which may trigger a positive climatic feedback (Blok et al 2011) as more warmth is absorbed by the darker surface. However, the magnitude and direction of vegetation response to changing climate is highly variable across space and time (Epstein et al 2013, Berner et al 2020. Decreasing vegetation greenness-so called 'browning'-has been linked, for example, to increasing surface water and tundra wetness in some regions (Lara et al 2018) as well as to earlier snowmelt (Gamon et al 2013).
Most remote sensing assessments of Arctic 'greening' or 'browning' have been based on changes in the normalized difference vegetation index (NDVI) derived from satellite sensors such as the Advanced Very-High Resolution Radiometer (AVHRR, 1981present) (Walker et al 2009, Bhatt et al 2010, Forbes et al 2010 or the Moderate Resolution Imaging Spectroradiometer (MODIS, 2000-present) (Dutrieux et al 2012, Miles and Esau 2016). The coarse spatial resolution of these sensors (250 m to ∼12.5 km) does not allow distinguishing different landscapelevel drivers of NDVI change (Myers-Smith et al 2020). For instance, Walker et al (2009) found little change in AVHRR GIMMSg NDVI from 1982 to 2007 on Yamal Peninsula, West Siberia, Russia (figure 1), but while in the field observed many local factors affecting vegetation, such as thermal erosion gullies, landslides, reindeer herding and infrastructure. Similarly, Bhatt et al (2010) found that GIMMSg NDVI was relatively stable on Yamal from 1982 to 2008, which they connected to the lack of positive trend in summer warmth index (sum of positive monthly mean temperatures) and/or to widespread reindeer herding.
Long-term moderate and high-resolution satellite image datasets, such as Landsat (1972-present), are needed to address the issue of spatiotemporal heterogeneity in vegetation productivity. Evaluation of Landsat NDVI trends 1999-2014 across the Arctic Coastal Plain of northern Alaska showed that the dominant factors controlling greenness were geomorphology and regional climate change (Lara et al 2018).
Trend analysis of Landsat imagery shows strong potential for detecting landscape change and attributing changes to specific factors such as shrub expansion (Fraser et al 2012, Brooker et al 2014. However, the short growing season, persistent cloud cover and data loss from the scan line corrector failure on Landsat ETM+ (Enhanced Thematic Mapper Plus) (USGS 2019a) limit the amount of usable imagery and may result in too sparse time series for robust trend analysis in many Arctic regions. In addition, the differences in designated spectral band widths of Landsat TM (Thematic Mapper) and ETM+ and , even more so between band widths of Landsat OLI (Operational Land Imager) in comparison to previous Landsat sensors (table 1), present further challenges to assessing multi-decadal trends in Landsat NDVI (Flood 2014).
Slope processes including landslides and subsequent enhancement of erect willow growth are important factors affecting NDVI on Yamal (Walker et al 2009). Cryogenic landslides, such as active layer detachments (ALD) and (retrogressive) thaw slumps, expose marine sediments (Ukraintseva and Leibman 2007) which leads to slightly higher groundwater salinity relative to background values (Khitun et al 2015). This provides colonizing plants with nutrients, thus enhancing revegetation and re-formation of soils. Furthermore, this fosters expansion of willows (Salix sp.), which also benefit from deeper active layer (seasonally thawing and freezing layer above the permafrost) and thicker snow cover in depressions formed by landslides (Leibman et al 2015). When landslides remove the surface layer they expose permafrost and, in the process, alter soil temperature (Loranty et al 2018), hydrology (Kokelj and Lewkowicz 1998) and nutrient content (Ukraintseva and Leibman 2007). The exposed permafrost can then thaw and release previously frozen soil organic matter for decomposition (Abbott andJones 2015, Olefeldt et al 2016). Changed soil conditions may stimulate both photosynthesis and respiration, thus impacting the carbon balance within the affected area (Abbott and Jones 2015). As landslides-especially thaw slumps-often occur next to water bodies (Abbott et al 2015), and due to the multitude of lakes, ponds and rivers on central Yamal Peninsula, the sediments and organic matter are likely to be transported further away from the actual landslide.
The last massive landslide event on Yamal took place during the end of summer 1989 when thousands of landslides, mostly ALDs, occurred in the north-west part of the peninsula. These landslides are believed to have resulted from high pore pressure caused by large winter and summer precipitation events combined with rather warm summer temperatures (Leibman and Egorov 1996). In addition to natural triggers, landslides on Yamal may result from anthropogenic activities, such as infrastructure construction and off-road vehicle traffic. Russia's largest natural gas deposits have been found on the Yamal Peninsula, and exploitation of this resource rapidly expanded during the past several decades (Kumpula et al 2010. Cryogenic landslides also create a potential threat to hydrocarbon extraction infrastructure. The peninsula is also home for about 1000 fully nomadic Nenets households with ca. 330 000 reindeer in total (Administration of Municipal Education of Yamalskii district 2017). Reindeer are the main large herbivore in the area and are attracted to wide, windexposed bare surface areas during periods of high insect harassment (Skarin et al 2010(Skarin et al , 2020. These animals can impede regeneration of vegetation after the disturbance via grazing and trampling (Kumpula et al 2010), yet the magnitude and degree of reindeer impact on landslide evolution is unknown. Khitun et al (2015) studied vegetation regeneration on ALDs in Central Yamal in situ. Their focus was on botanical surveys and phytosociological studies on landslides that occurred in 1989, the mid 1970's and the 1950's, and ca. 1000-year old ancient landslides, as well as their surroundings. According to their observations, the first vegetation was dominated by graminoids that occur on shear surfaces during the first 10-15 years. The first willows appear approximately 15 years after the surface failure and the transition to sedge-willow vegetation cover occurs after another 35-40 years. Cannone et al (2010) followed vegetation colonization of active layer detachments on Ellesmere Island in Canadian High Arctic. There it required over 50 years to develop floristic composition similar to surrounding communities, but the time elapsed was not yet sufficient enough to resemble the original vegetation cover. In Canada, Tasseled Cap Trend Analysis from Landsat image stack was also used by Brooker et al (2014) to study the  Commercial very high-resolution (VHR; 1 to 5 m) optical satellite imagery for example from IKONOS, QuickBird and WorldView satellites can provide more detailed insight into landscape dynamics; however, it has been used to a lesser extent due to the high costs. VHR imagery has been used widely to assess revegetation after forest fires (Mitchell and Yuan 2010, Mitri and Gitas 2013, Chu et al 2016, Fang et al 2019, but regarding slope disturbances the applications are mainly in the detection and mapping of landslides (Lu et al 2011, Rudy et al 2013, Murillo-García et al 2015, Amatya et al 2019 rather than revegetation analyses. To the best of our knowledge, no landscape-scale remote sensing techniques have been used to connect revegetation dynamics on landslides to local and regional NDVI-patterns. The aim of this research was to assess vegetation regeneration using NDVI from moderate-and very high spatial resolution satellites in an area affected by a large number of landslides. We evaluate the applicability of Landsat satellite image time series for assessing revegetation of landslide shear surfaces over a 30year period. We used Google Earth Engine (Gorelick et al 2017) in combination with novel approach for cross calibration of data from different sensors (Berner et al 2020) to create robust NDVI time series. We also used three VHR satellite images for more detailed analysis of NDVI change within our study area in Central Yamal. In particular, we looked at (1) how vegetation regeneration on landslides affects landscape level NDVI dynamics in low Arctic tundra; and (2) how spatial resolution affects the utility of optical multispectral satellite data in detection of landscape level factors affecting regional NDVI.

Study area
The Yamal Peninsula, West Siberia (figure 1) is about 700 km long and 150 km wide with continuous, often ice-rich, permafrost. The prevailing landscape topography is rather flat with elevations only up to ca. 80 m above sea level (a.s.l.), and is dominated by sandy ridges, mires and wetlands, ponds, lakes, rivers and river valleys. Figure 1 shows the Mordy-Yakha study area (70 • 10 ′ N, 68 • 31 ′ E) of ca. 35 km 2 that was chosen for detailed investigation. Our team has conducted extensive field work and interviews with local Nenets reindeer herders in this region over the past couple decades ( The shape of the study area was determined by the extent of two ridge/landslide complexes deemed representative of Central Yamal (Leibman and Egorov 1996). The area consists of one large, ca. 9 km long north-south ridge, and one smaller, 3 km long ridge. The maximum altitude is ca. 60 m a.s.l. The deposits are fine-grained, sandy to clayey, marine and alluvialmarine deposits. The landslides occur: (1) on slopes adjacent to sandy ridge tops where the clay above the permafrost is overlain by sand, which comprises most of the active layer; and (2) on slopes adjacent to clayey ridge tops with a more clayey composition of the entire active layer (Leibman and Egorov 1996).
The study area belongs to bioclimatic subzone D of the Raster Circumpolar Arctic Vegetation Map, while G3 (non-tussock sedge, dwarf-shrub, moss tundra) and S1 (erect dwarf-shrub tundra) are the main map units within the study area (Raynolds et al 2019). Dwarf shrub-moss-lichen vegetation occupies the well-drained sandy ridge tops and slopes that are not affected by landslides. Clayey parts of the ridge tops are mostly covered with sedge and sphagnum mires and flat-topped polygonal peatlands. Almost 30 years after the initial landslide event, the vegetation on affected slopes is dominated by grasses and low willow shrubs about 10 cm high. Valleys and old landslide cirques, where snow accumulates during winter, are typically dominated by tall willows up to 150 cm high.

Climatic conditions
The nearest weather station with available data was Marre-Sale (figure 1), which is located on the coast of the Kara Sea about 80 kilometres south-west from our study area. The data were first accessed through the National Climatic Data Center's Daily Summaries from (NOAA 2018). Some of the gaps in the precipitation data were filled with data for the same station that is available from rp5.ru (https://rp5.ru/Weather_archive_in_Marresale). Figure 2 shows the mean annual air temperature (MAAT), which increased 0.06 • C yr −1 between 1960 and 2017 (P < 0.001, r 2 = 0.29), rising from an average of −8.4 • C (1961-1990) to −6.2 • C (1988-2017). The average annual precipitation for the same period was 312 ± 54 mm yr −1 (±1 SD), of which approximately half fell as snow and half as rain (Leibman et al 2015). Linear regression shows that annual precipitation increased by 1.1 mm yr −1 or approximately 11 mm per decade (P = 0.011, r 2 = 0.124), though varied notably between years (figure 2).

Landslide delineation
To delineate the landslide events of 1989, we searched the U.S. Geological Survey Earth Explorer Server (https://earthexplorer.usgs.gov/) for contemporaneous Landsat scenes that were free of clouds, shadows and artefacts. We used two high-quality images from around the time of peak vegetation greenness (Landsat 4 TM from July 1988 and 5 TM from July 1990 (table 2)) and downloaded these as Level-1 products (USGS 2019b). The pixel values of all images were converted from radiance to Top-of-Atmosphere (TOA) reflectance. ArcMap version 10.5 (ESRI, Inc. USA) and ERDAS IMAGINE 2016 (Hexagon Geospatial, USA) software were used in satellite data processing.
We delineated the extent of the 1989 landslides based on changes in NDVI between the Landsat TM images from 1988 and 1990. The NDVI is a widely used index of vegetation greenness (Pettorelli et al   Laidler and Treitz 2003). Pixels with change of 0-(−0.35) were converted to vector polygons, from which visually-determined undisturbed (non-landslide) pixels were deleted and remaining polygons adjusted based on a 1993 Satellite Pour l'Observation de la Terre (SPOT) image, which had higher spatial resolution than Landsat TM (table 2). These polygons were used in the subsequent analysis of NDVI change on the shear surfaces and surrounding unaffected areas.

Landsat NDVI time series
We assessed landslide impacts and subsequent revegetation using time series of Landsat NDVI from both disturbed and adjoining undisturbed tundra. This first involved generating 500 random sample sites split evenly between disturbed and nearby undisturbed areas as identified with the landslide polygons. To avoid generating sample sites that straddled the border between these areas, we converted the landslide polygons to polylines, buffered the polylines by 30 m in each direction, and then erased the buffered areas from the broader area of interest. We generated the sample sites using the spatialEco package (Evans 2019) in the statistical software R (R Core Team 2018), and then uploaded as a shapefile to Google Earth Engine (GEE), a cloud-computing platform for planetary-scale remote sensing analyses (Gorelick et al 2017).
Landsat surface reflectance data were acquired for each sampling site and then used to assess how summer NDVI varied through time in landslide and nonlandslide affected areas. Specifically, we downloaded all May through September (1985-2018) Landsat 5, 7, and 8 Collection 1 (Tier 1) surface reflectance data (Masek et al 2006, Vermote et al 2016 for each sampling site using the GEE Python API. We then excluded observations affected by snow, water, clouds, and cloud shadows using the CFmask quality  (Pekel et al 2016). We further excluded observations from scenes with high cloud cover (>80%) or that were acquired under high solar zenith angle (>60 • ). Next, we computed NDVI for each remaining clear-sky observation and then cross-calibrated NDVI from Landsat 5 and 8 with Landsat 7 using a machine learning random forest algorithm (Breiman 2001) previously trained with Landsat measurements from tundra ecosystems around the Arctic (Berner et al 2020). We observed that NDVI at these sites typically peaked each summer in early July and was then stable through August before declining in autumn. We therefore computed yearly summer NDVI for each site by averaging NDVI of observations that were acquired during July and August.

Very high-resolution satellite data
The VHR-dataset included QuickBird-2 (2004), WorldView-2 (2013) and WorldView-3 (2017) imagery, hereafter referred as QB-2, WV-2 and WV-3. These data were received as Ortho Ready Standard Imagery 2A products (DigitalGlobe Inc. of Longmont, CO, USA). All VHR images were orthorectified based on TanDEM-X (German Aerospace Center, DLR) 12 m Digital Elevation Model. The georeferencing of QB-2 was corrected to align with more geometrically accurate WV images (RMSE 1.46 m). After the geometric corrections, the pixel values of all three images were converted from radiance to TOA reflectance. The acquisition dates and spatial resolutions are shown in table 2. The NDVI layers from QB-2, WV-2 and WV-3 were divided into 20 classes using natural breaks to determine values representing bare or semi-bare land surface within 1989 landslide polygons. Such NDVI values were 0.01-0.25, 0.04-0.29 and 0.03-0.22 for 2004, 2013 and 2017, respectively. The values differ slightly between years due to differences in sensors (table 2), and possibly also due to varying temperature and moisture conditions between years. The NDVI changes 2004-2013, 2013-2017and 2004-2017 were also calculated. New landslides, mainly thaw slumps along the shores of small lakes and ponds, were manually digitised from all VHR images by visual interpretation.

Landslide delineation and impacts assessed with landsat
Approximately 65 new landslides appeared within our study area in 1989, which triggered a pronounced decline in tundra greenness (NDVI) (figures 3 and 4). These landslides varied in extent from ca. 0.2 ha up to ca. 16 ha, with most (80%) being . The landslide event in 1989 caused a pronounced decline in tundra greenness (NDVI) which then gradually increased over subsequent decades, though it remained lower than non-landslide affected areas in 2018. Solid lines depict the annual average of mean summer NDVI for random sample sites (n = 250) located in either landslide or non-landslide affected areas, with shaded bands encompassing one standard deviation in NDVI among sampling sites each year. Years with observations from fewer than 80% of sample sites were excluded. smaller than five hectares. In total, the surface disturbance affected ca. 300 hectares within the study area. Landsat mean summer NDVI time series for the 500 sample sites shows the decrease in NDVI within landslide polygons after 1989 and regeneration of vegetation during the following decades (figure 4). Sites that subsequently experienced landslides had slightly lower mean NDVI from 1985 to 1988 (0.572 ± 0.034 NDVI) than undisturbed sites (0.597 ± 0.033 NDVI) (t-test; t = 8.05, P < 0.001). Following the landslides, the difference in mean NDVI between affected and unaffected areas was 0.338 ± 0.009 in 1990, with the difference gradually narrowing to 0.065 ± 0.003 by 2018. Mean NDVI in disturbed areas increased at a rate of 0.009 NDVI yr −1 from 1990 to 2013 (P < 0.001, r 2 = 0.77), though declined slightly, albeit not significantly, from 2013 to 2018 (P = 0.23, r 2 = 0.16). Undisturbed areas showed no trend in mean NDVI from 1985 to 2018 (P = 0.58), though NDVI declined at a rate of −0.008 NDVI yr −1 from 2013 to 2018 (P = 0.025, r 2 = 0.78).

Landslide activity and vegetation recovery mapped with VHR data
Landslide shear surfaces were still readily detectable from QuickBird-2 and even from WorldView-2 and WorldView-3 images despite the amount of time that had passed since the surface disturbance occurred (15, 24 and 28 years, respectively). Yet, it was difficult to evaluate revegetation from VHR optical satellite data by means of visual interpretation, as the relative difference in spatial resolution was considerable (table 2). Ten new landslides, mostly slope failures adjacent to lakes and ponds, were identified from 2004 QB-2 imagery, while 44 new landslides were identified from 2013 WV-2 imagery and nine from 2017 WV-3 imagery (figure 5). These landslides were considerably smaller than the 1989 ALDs; ca. 65% of landslides were smaller than 0.1 ha, while the largest was about 1.6 ha. They affected around 0.05%, 0.23% and 0.22% of the land area within the study region in 2004, 2013 and 2017, respectively. Nine of the new landslides evident in 2004 had visually increased in area by 2013 (median = 87.3%, min = 2%, max = 644.8%), but only four of them also increased between 2013 and 2017. Twenty-four of the 44 new landslides from 2013 increased in area by 2017 (median = 50.7%, min = 8.3%, max = 714.4%). Of those landslides that occurred between 2004 and 2013 on shore of a small lake, pond or stream, 87% had increased in surface area by 2017, while only 20% of those that occurred 'inland' had increased during the same period.
Based on the classification of VHR NDVI data, landslides within the study area had little or no vegetation cover across 1.03 km 2 in 2004, 0.58 km 2 in 2013 and 0.48 km 2 in 2017 (the areas within the landslide polygons with NDVI 0.01-0.25, 0.04-0.29 and 0.03-0.22, respectively). Based on these values, vegetation regeneration had advanced 0.45 km 2 during the nine-year period 2004-2013 (0.05 km 2 yr −1 ) between the QB-2 and WV-2 images, and 0.1 km 2 (0.025 km 2 yr −1 ) between WV-2 and WV-3 images. The revegetation of landslide shear surfaces clearly stands out in the NDVI-change between 2004 QB-2 and 2013 WV-2 (figure 6), despite the differences in spectral and spatial resolution between the two respective sensors (tables 1 and 2). In 2017, the NDVI-values within the study area were generally lower than in 2013, especially in thermoerosion gullies and small channels on the slopes affected by the 1989 landslides. The greatest increase of NDVI values between 2013 and 2017 occurred within the same or adjacent areas to those that displayed the most pronounced decrease during the period 2004-2013. Some of these places were apparently reindeer herders' campsites, indicated by a small, heavily trampled area and the presence of herders' dwellings and reindeer corralling sites (figure 7). This heavy trampling affected around 0.22% in 2013 and 0.25% in 2017 of total land area within our study region. The meteorological records for these three years (table 3) show that precipitation prior to the acquisition of WV-2 imagery was almost double the amount of that prior to the QB-4 and WV-3 images, while 2017 had the lowest precipitation and highest number of thaw degree days (sum of positive temperatures) prior to the acquisition of the image.

Discussion
Our analysis revealed that massive landslides in 1989 caused a pronounced decline in local vegetation cover and that these areas gradually regained vegetation but did not fully regenerate to pre-disturbance levels during the following three decades. The study by Khitun et al (2015) revealed that during the first 10-15 years, the shear surfaces of landslides in Central Yamal become populated by pioneer vegetation dominated by halophytic graminoids, bryophytes and herbs. The first willow seedlings appeared after 15-17 years, and sedge-willow vegetation began to dominate ca. 50-55 years after a landslide event. Our study demonstrates that extensive active-layer detachment slides remained clearly visible in Landsat imagery for approximately a decade, though after 20 years the landslide shear surfaces were rendered nearly indistinguishable from the surrounding landscape (figure 3), especially if the interpreter has no ground-level knowledge. This means that even the pioneer vegetation and young willows are enough to make the landslides barely detectable by visual analysis of Landsat data. However, the NDVI was still ∼10% lower than the background values almost 30 years after the initial disturbance, which allows us to distinguish between areas affected and not affected by the 1989 landslides. As of 2018, the mean NDVI on ALD shear surfaces was 0.533 ± 0.048, compared to 0.589 ± 0.035 in undisturbed areas. The difference in NDVI between shear and undisturbed surfaces did not change much during the final decade of our time series. However, vegetation regeneration has progressed enough to support low-erect willows, as observed in the field. The Landsat time series revealed Table 3. Thaw degree days (sum of positive air temperatures) and precipitation at the time of acquisition of QB-2 (2004), WV-2 (2013) and WV-3 (2017). Maximum values in bold. Precipitation was calculated from the first day with temperature above 0 • C in 2013 (DOY 131, May 11th). 15.7.2004 21.7.2013 21.7.2017 Thaw degree days ( • C) 281.6 265.1 327.4 Cumulative precipitation from DOY 131 (mm) 20.9 44.7 15.7 Figure 6. NDVI change between the respective acquisition periods of very high-resolution imagery covering the years 2004-2013, 2013-2017, and 2004-2017. no trend in NDVI in the surrounding area from 1985 to 2018; however, NDVI declined on both landslide surfaces and surrounding areas between 2013 and 2018. The larger inter-annual variability in NDVI on landslide surfaces could be due difference in surface temperature and hydrology. The VHR imagery portrayed the revegetation of landslide surfaces in much greater detail than Landsat and allowed us to not only identify small landslides, but also to quantify gradual vegetation succession from bare and semi-bare soil surfaces to more closed vegetation cover. However, one must be careful when investigating the NDVI dynamics from only a few images with very high spatial resolution, as the potential effect of annual climatic variations is higher than in more temporally dense time series. Based on the NDVI-values of VHR data, bare and semi-bare landslide surfaces within the key study area decreased by 44% (0.45 km 2 ) within nine years between 2004 and 2013. July of 2017 was much drier than in 2013, and the NDVI values dropped between 2013 and 2017 especially where they had increased 2004-2013. Walker et al (2009 noted that the NDVI dynamics in Yamal are affected by interactions among ecological and societal factors. In our study area, some of the NDVI-change 'hot-spots' were reindeer herders' campsites, where hundreds of reindeer gathered in a small area and heavily trampled the vegetation, on which graminoid cover regenerated during the following year(s).
The VHR imagery revealed abundant green vegetation in 2013 relative to 2004 and 2017. Furthermore, NDVI during summer of 2017 was generally lower than in 2013 both in the WV and Landsat imagery. Yamal Peninsula experienced a record long, warm and dry period in July of 2016, and a similar warm, dry period preceded the acquisition of WV-3 and Landsat images in July 2017. The heat wave and drought likely led to particularly dry conditions on the clayey shear surfaces and contributed to drainage channels displaying a stronger decrease in NDVI in comparison to the surrounding tundra. This suggests that the extreme intra-seasonal and inter-annual climatic fluctuations may facilitate or delay vegetation growth. Hot and dry conditions have also contributed to declining shrub growth in Greenland (Gamm et al 2018) and to episodes of low tree growth along the forest-tundra ecotone in northeastern Siberia (Berner et al 2013). The importance of soil moisture for shrub growth has also been highlighted by Myers-Smith et al (2015). If summer heat waves and droughts become more common, then the likely result is that the shear surfaces, and regional tundra vegetation in general, could become less productive. Further  al 2003)). Also, a clustered reindeer herd is clearly visible. research including remote sensing, field studies and local (indigenous) knowledge is required to determine if the decline in NDVI observed in our study area is a signal of tundra 'browning' as has been observed in, for example, northern Alaska (Lara et al 2018), and whether this shift will have consequences for regional carbon uptake and reindeer forage.
While no new large-scale ALDs were observed within our study area since 1989, many new thaw slumps occurred especially between 2004 and 2013, while between 2013 and 2017 fewer new slides occurred. However, the extent of most pre-existing thaw slumps, especially along the shores of small lakes, ponds and streams, increased as the thawing of ground ice continued. This is consistent with observations by Khomutov et al (2017). They noted that the number and area of thaw slumps have increased in Central Yamal since 2012, while the prevailing meteorological conditions have most likely prevented the occurrence of additional ALDs. Detachment of the active layer usually requires an icy transient layer to develop at the base of the active layer. Warming conditions have contributed to deeper thaw, which has thawed ground ice within the transient layer, thus preventing the ALDs from occurring (Khomutov et al 2017). The relatively warm summers of 2012 and 2013 likely resulted in seasonal thaw reaching the top layers of tabular ground ice, causing an increase in the number of thaw slumps, which we were able to observe from the WV-2 image. Regular monitoring of landslide activity is important not only for understanding ongoing changes in the Arctic environment, but also to minimize potential impacts on infrastructure (Nelson et al 2001), which can have large-scale negative (socio-)ecological implications in this sensitive tundra environment (Forbes et al 2014). The region has experienced rapid expansion of hydrocarbon extraction activities and related infrastructure since 1988 . For example, giant gas field Bovanenkovo is located only about 20 km north from our study area.
The processes analyzed here are not unique to Northern Russia. For example, widespread ALDs have been reported from Ellermere Island, Canada (Lewkowicz 1990). Other studies have found increasing degradation of ice-rich permafrost (thermokarst) in parts of northern Canada (Lantz and Kokelj 2008, Lantuit and Pollard 2008, Kokelj et al 2017, Lewkowicz and Way 2019 and Alaska (Jorgenson et al 2006) during last half century. In North-America, permafrost degradation, especially in a form of retrogressive thaw slumping, has been related to widespread buried glacial ice (Segal et al 2016) and massive ice wedges (Jorgenson et al 2006), which have been affected by warming air temperatures since 1980's. Widespread ice wedge degradation across the Arctic was also reported by Liljedahl et al (2016). Continuing warming of the Arctic (Post et al 2019) will likely lead to an increasing number and volume of thaw slumps. Altered micro-climate and soil conditions coupled with warming climate may enhance vegetation growth in areas affected by thaw slumps (Lantz et al 2009), as has similarly been observed for ALDs in Central Yamal (Khitun et al 2015). The extent to which permafrost degradation influences satellite observations of greening and browning across the Arctic remains an important unanswered question.

Conclusions
The aim of this study was to assess vegetation recovery on landslide slopes in the Low Arctic tundra zone of Yamal, West Siberia, where active layer detachments and retrogressive thaw slumps are widespread. The analysis of Landsat time series revealed that the regeneration of vegetation on shear surfaces rapidly progressed during the first two decades after the last major landslide event in 1989; however, NDVI was still lower in the landslide affected areas as of 2018, which indicates that vegetation has not fully regenerated after ∼30 years. Regeneration to original levels of composition, cover and productivity may require several more decades. The Landsat time series analysis also revealed that NDVI declined from 2013-2018 on both landslide and unaffected areas, possibly due to high summer temperatures and drought stress. These recent declines in NDVI warrant further investigation.
The very high-resolution imagery used here represents meteorologically contrasting growing seasons, which makes it challenging to separate the general NDVI trend from inter-annual variability. The use of such highly detailed imagery, however, allowed the detection of NDVI change 'hot-spots' indicating temporary decreases in vegetation productivity, as well as detection of new and expanding thaw slumps, which are too small to be detected from moderate-resolution satellite data like Landsat. Thaw slumps are expected to become more common as ice rich permafrost degrades while Arctic climate continues to warm, but they are less likely to substantially affect regional NDVI due to their relatively small size in comparison to more extensive, often massive, active layer detachments. Cryogenic landslides are widespread in permafrost regions and have pronounced effects on local vegetation that could contribute to the satellite observations of both greening and browning in the Arctic.

Data availability statement
The data that support the findings of this study are available upon reasonable request from the authors.