Land Cover Change in the Lower Yenisei River Using Dense Stacking of Landsat Imagery in Google Earth Engine

Climate warming is occurring at an unprecedented rate in the Arctic due to regional amplification, potentially accelerating land cover change. Measuring and monitoring land cover change utilizing optical remote sensing in the Arctic has been challenging due to persistent cloud and snow cover issues and the spectrally similar land cover types. Google Earth Engine (GEE) represents a powerful tool to efficiently investigate these changes using a large repository of available optical imagery. This work examines land cover change in the Lower Yenisei River region of arctic central Siberia and exemplifies the application of GEE using the random forest classification algorithm for Landsat dense stacks spanning the 32-year period from 1985 to 2017, referencing 1641 images in total. The semiautomated methodology presented here classifies the study area on a per-pixel basis utilizing the complete Landsat record available for the region by only drawing from minimally cloudand snow-affected pixels. Climatic changes observed within the study area’s natural environments show a statistically significant steady greening (~21,000 km2 transition from tundra to taiga) and a slight decrease (~700 km2) in the abundance of large lakes, indicative of substantial permafrost degradation. The results of this work provide an effective semiautomated classification strategy for remote sensing in permafrost regions and map products that can be applied to future regional environmental modeling of the Lower Yenisei River region.


Introduction
The high-latitude regions of Eurasia are warming at approximately 0.12 • C per year, significantly faster than the global average (e.g., [1][2][3]). Environmental changes associated with a warming climate have significant impacts on arctic and subarctic ecosystems, including surface and subsurface hydrology. Numerous observational studies indicate that such changes related to the structure, composition, and functioning of tundra and boreal forest biomes are ongoing. For example, increased photosynthetic productivity under warming climatic conditions, derived from remotely sensed normalized difference vegetation index (NDVI) data, frequently referred to as "arctic greening," has been reported for several Eurasian arctic regions (e.g., [4][5][6][7][8]). One of the major drivers of the observed greening trend is the increased abundance of shrub species in tundra ecosystems (e.g., [9][10][11][12][13]). Lengthening growing season, reduced seasonality, thickening of the active layer (the uppermost permafrost-affected soil layer that thaws seasonally), and increased annual snow depths are the Figure 1. Study area map. The study area is defined by two Landsat scene extents (black boxes). The area observed encompasses a transitional zone between physiographic provinces (West Siberian Plain and Central Siberian Plateau) and biomes (taiga and tundra) centered on the Lower Yenisei River just before its outlet into the Kara Sea.
The region has a subarctic climate, classified as Dfc in the Köppen climate classification system. Few weather stations with adequate observational records are located within this large and diverse study area. Based on weather stations in the cities of Norilsk and Igarka within the study area, though, the mean annual temperature of the region is between −8.2 °C and −10.2 °C, based on 30-year climatic normals. The region has just 2 months per year when mean monthly air temperature (MMAT) is greater than 10 °C and 4 months when MMAT is greater than 0 °C. The hottest month of the year (July) has MMAT between 8.5 °C and 20 °C and the coldest month (January) ranges between -22.8 °C and −33.5 °C. Annual precipitation totals range from 399 mm to 474 mm, with greater amounts typically falling during the warm season. Mean annual air temperatures have been increasing since 1985 at an average rate of 0.05 °C/year in Igarka, located in the southern part of the study area, and 0.045 °C/year in Norilsk, located in the northern part. The total annual precipitation has also been increasing, in Igarka at an average rate of 2.4 mm/year and in Norilsk at an average rate of 2.7 mm/year. The area observed encompasses a transitional zone between physiographic provinces (West Siberian Plain and Central Siberian Plateau) and biomes (taiga and tundra) centered on the Lower Yenisei River just before its outlet into the Kara Sea.
The region has a subarctic climate, classified as Dfc in the Köppen climate classification system. Few weather stations with adequate observational records are located within this large and diverse study area. Based on weather stations in the cities of Norilsk and Igarka within the study area, though, the mean annual temperature of the region is between −8.2 • C and −10.2 • C, based on 30-year climatic normals. The region has just 2 months per year when mean monthly air temperature (MMAT) is greater than 10 • C and 4 months when MMAT is greater than 0 • C. The hottest month of the year (July) has MMAT between 8.5 • C and 20 • C and the coldest month (January) ranges between -22.8 • C and −33.5 • C. Annual precipitation totals range from 399 mm to 474 mm, with greater amounts typically falling during the warm season. Mean annual air temperatures have been increasing since 1985 at an average rate of 0.05 • C/year in Igarka, located in the southern part of the study area, and 0.045 • C/year in Norilsk, located in the northern part. The total annual precipitation has also been increasing, in Igarka at an average rate of 2.4 mm/year and in Norilsk at an average rate of 2.7 mm/year. The lower Yenisei River basin is occupied by permafrost of variable continuity. Based on the USSR permafrost map by Yershov et al. [43], discontinuous permafrost occupies 10.3% of the study area, while the remainder of the area, except for large water bodies, is underlain by continuous permafrost (Figure 2). Permafrost temperature at the zero annual amplitude depth (~10 m) is highly variable and ranges from 0 • C to −10 • C depending on the location and site-specific edaphic conditions within the region. The lower Yenisei River basin is occupied by permafrost of variable continuity. Based on the USSR permafrost map by Yershov et al. [43], discontinuous permafrost occupies 10.3% of the study area, while the remainder of the area, except for large water bodies, is underlain by continuous permafrost (Figure 2). Permafrost temperature at the zero annual amplitude depth (~10 m) is highly variable and ranges from 0 °C to −10 °C depending on the location and site-specific edaphic conditions within the region.

Figure 2.
Permafrost distribution within the study area (digitized from Yershov et al. [43]). Discontinuous permafrost occupies 10.32% of the study area, primarily in the southernmost limits. The remaining area is underlain by continuous permafrost, with the exception of large lakes. Yershov et al. further defined continuous permafrost by its temperature at the depth of zero annual amplitude (~10 m), which varies from 0 °C to −10 °C.
The study area represents the transition from Siberian taiga, primarily occupying areas in the southern portion, and tundra biomes, primarily in the north. Field surveys conducted throughout the area during the summer from 2010 to 2015 indicate that taiga forests are dominated by Siberian larch (Larex sibirica), Siberian spruce (Pecea ovobata), and silver birch (Betula pendula). Tundra landscapes Figure 2. Permafrost distribution within the study area (digitized from Yershov et al. [43]). Discontinuous permafrost occupies 10.32% of the study area, primarily in the southernmost limits. The remaining area is underlain by continuous permafrost, with the exception of large lakes. Yershov et al. further defined continuous permafrost by its temperature at the depth of zero annual amplitude (~10 m), which varies from 0 • C to −10 • C.
The study area represents the transition from Siberian taiga, primarily occupying areas in the southern portion, and tundra biomes, primarily in the north. Field surveys conducted throughout the area during the summer from 2010 to 2015 indicate that taiga forests are dominated by Siberian larch (Larex sibirica), Siberian spruce (Pecea ovobata), and silver birch (Betula pendula). Tundra landscapes are characteristic of elevated peatlands covered by mosses (Sphagnum genus), lichens (Cetraria islandica), grasses (Poa spp.), and low shrubs such as blueberries (Vaccinium uliginosum), dwarf birch (Betula nana), and alder (Alnus fruticosa). Critical within this analysis is the relationship between specific vegetative covers and permafrost continuity established in the southern portion of the study area by Tyrtikov [44] and Rodionov et al. [45]. Here, tundra developed on lacustrine clays with thick organic horizons is underlain by 10-40 m thick near-surface Holocene-age permafrost. Near-surface permafrost is, however, absent under forests. In forested areas, the table of older Pleistocene-age permafrost is located at 3-5 m depth. This vegetation-near-surface permafrost relationship allows us to use land cover change as an indicator of permafrost change.

Landsat Imagery and Composite Generation
Using GEE, an image collection, or data stack, was produced for 3-year periods comprising all images intersecting the study area from July 1 to September 30 to produce a cloud-free composite of scenes with minimal snow cover. The entirety of the Landsat 5, 7, and 8 archives available for this area was included in this analysis using the tier 1 top of atmosphere (TOA) reflectance product, described in [46]. The reflectance product is preferred over the TOA radiance because the reflectance algorithm removes the exoplanetary effects associated with variable solar irradiance as a function of variability in (1) solar zenith angles, (2) spectral band differences, and (3) Earth-to-Sun distance at different times of the year. The constants for all Landsat sensors (TM, ETM+, OLI) are tabulated in [46], where an overview of the calibration procedure is provided.
The use of alternative surface reflectance from the Landsat Ecosystem Disturbance Adaptive Processing System (LEDAPS) product was avoided, since users are cautioned against applying it in high latitudes (>65 • ) [47]. This method includes pixels from overlapping images from adjacent acquisitions, including the addition of ascending orbits for Landsat 8, with specific paths and rows for each 3-year period detailed in Table 1 and the processing chain shown in Figure 3. The defined study size comprises the combined footprint of Landsat path 153, rows 11 and 12 (within the Landsat WRS-2 orbit pattern) ( Figure 3, step 2). The total period of scene acquisition dates ranged from 1985 through 2017, resulting in 11 3-year image composites with corresponding dates, provided in Table 1. are characteristic of elevated peatlands covered by mosses (Sphagnum genus), lichens (Cetraria islandica), grasses (Poa spp.), and low shrubs such as blueberries (Vaccinium uliginosum), dwarf birch (Betula nana), and alder (Alnus fruticosa). Critical within this analysis is the relationship between specific vegetative covers and permafrost continuity established in the southern portion of the study area by Tyrtikov [44] and Rodionov et al. [45]. Here, tundra developed on lacustrine clays with thick organic horizons is underlain by 10-40 m thick near-surface Holocene-age permafrost. Near-surface permafrost is, however, absent under forests. In forested areas, the table of older Pleistocene-age permafrost is located at 3-5 m depth. This vegetation-near-surface permafrost relationship allows us to use land cover change as an indicator of permafrost change.

Landsat Imagery and Composite Generation
Using GEE, an image collection, or data stack, was produced for 3-year periods comprising all images intersecting the study area from July 1 to September 30 to produce a cloud-free composite of scenes with minimal snow cover. The entirety of the Landsat 5, 7, and 8 archives available for this area was included in this analysis using the tier 1 top of atmosphere (TOA) reflectance product, described in [46]. The reflectance product is preferred over the TOA radiance because the reflectance algorithm removes the exoplanetary effects associated with variable solar irradiance as a function of variability in (1) solar zenith angles, (2) spectral band differences, and (3) Earth-to-Sun distance at different times of the year. The constants for all Landsat sensors (TM, ETM+, OLI) are tabulated in [46], where an overview of the calibration procedure is provided.
The use of alternative surface reflectance from the Landsat Ecosystem Disturbance Adaptive Processing System (LEDAPS) product was avoided, since users are cautioned against applying it in high latitudes (>65°) [47]. This method includes pixels from overlapping images from adjacent acquisitions, including the addition of ascending orbits for Landsat 8, with specific paths and rows for each 3-year period detailed in Table 1 and the processing chain shown in Figure 3. The defined study size comprises the combined footprint of Landsat path 153, rows 11 and 12 (within the Landsat WRS-2 orbit pattern) ( Figure 3, step 2). The total period of scene acquisition dates ranged from 1985 through 2017, resulting in 11 3-year image composites with corresponding dates, provided in Table  1. An initial classification test revealed that residual striping caused by the scanline corrector failure onboard Landsat 7 after May 31, 2003, impacted the classification. However, the inclusion of adjacent scenes mitigated the number of pixels removed from the analysis (note the scene extents indicated in Figure 3, step 2). After the cloud and shadow mask were applied, NDVI was calculated for each image in the collection, producing an additional image collection of NDVI values (Figure 3, step 5).
Two image composites were then created by flattening the image collections, Landsat reflectance, and NDVI. For the 3-year period, the median reflectance values for all reflected optical bands was extracted, comprising bands 1-5 and 7 for Landsat 5 and 7, and bands 1-9 for Landsat 8 (for a description of the wavelengths of each band, visit https://landsat.usgs.gov). Simultaneously, the maximum NDVI value for each pixel was extracted from the NDVI image collection to identify the "greenest," or most lush, vegetation observed during the time period to provide further evidence of the land cover type. The subsequent NDVI image composite was appended to the reflectance composite as a band. This resulted in an image composite in which each pixel represented the median reflectance and maximum NDVI value for the 3-year period ( Figure 3, step 6). In total, 11 3-year composites were produced, spanning the years 1985 to 2017, selected based on the availability of images within each image collection sufficient for data gaps, cloud cover compensation, and satisfactory land cover classification. Five composites from 1985 to 1999 use acquisitions from Landsat 5, 5 composites from 2000 to 2014 use those from Landsat 7, and a single composite from 2015 to 2017 use Landsat 8 acquisitions (Table 1).

Landsat Classification and Accuracy Assessment
The 3-year composites of Landsat imagery were then classified within GEE using the random forest classification method with the following 5 land cover categories: (1) barren ground, (2) taiga or closed forest, (3) tundra and open forest, (4) water, and (5) snow.
The random forest classification method is a learning algorithm that builds an ensemble of classification and regression-tree (CART) classifiers that use bagging (otherwise known as bootstrapping) of pixels within training polygons representing the spectral signatures of various land cover classes [48]. The algorithm uses random vectors to search through the input variables to establish the splits in nodes of trees. Classes are derived through a voting process of the ensemble of trees, whereby the most popular class at the node is selected. In the random forest algorithm, the trees are not pruned because of the bagging process (which lowers the potential for overfitting), thereby reducing the computational requirements and improving temporal efficiency. In this study, 30 trees were trained with approximately 100,000 randomly sampled points from user-defined training polygons.
An in situ dendrochronology study conducted in the nearby Putorana Mountains located within the study area [49] demonstrated that the determining factor in distinguishing between "open" and "closed" forests is the density of individual trees. A closed forest has abundant single-tree growth in-filled by shrubs, while an open forest has sparse single-stem or clustered multistem trees growing 10 to 20 m apart with herbaceous and moss vegetation in between. At the 30 m Landsat resolution, the open forest land cover class is indistinguishable from tundra and was therefore categorized as a single class.
The training sites/polygons representing the land cover types of interest were interpreted from high-resolution imagery available on Google Earth that was consistent across the time range of the study  and verified by ground validation and primary knowledge of conditions at field sites visited in July 2010-2015. To avoid accuracy inflation, sample points for ground validation sites were produced from an additional set of user-defined polygons separate from those used for training the classification algorithm. The accuracy reports for each classified time period included a confusion matrix and estimates of overall accuracy, kappa coefficient, user accuracy, and producer accuracy for each land cover class. These metrics were calculated using equations provided by Foody [50] and Congalton and Green [51] and are shown in Table 2 for selected classification periods, with the full time series provided in Appendix A.
The urban land cover class was omitted from the supervised random forest classification due to the many roads and structures resembling the spectral signature of barren land, resulting in misclassification. Open source vector data from Open Street Map (OSM) (www.openstreetmap.org) was incorporated to help distinguish urban footprints, which cover a very small portion of the study area (approximately 0.1%). Local administrative unit polygons were downloaded from OSM and a conditional statement was applied postclassification to convert pixels classified as barren ground to "built-up" when located within the municipalities of Igarka, Dudinka, Norilsk, and smaller settlements in the study area. The conditional statement was applied using the Raster Calculator tool in ArcMap

Results
The random forest classification yielded land cover maps representative of 11 three-year periods spanning the 32 years from 1985 to 2017 for the Lower Yenisei study area. Four selected examples of these classified composites, one per decade (1985-1987, 1994-1996, 2000-2002, 2015-2017), are shown in Figure 4. The accuracy assessments for all 11 classified maps yielded kappa statistics ranging from 90.83% to 98.5%. The confusion matrices (Table 2) indicate that the three categories snow, barren, and tundra/open forest were more frequently misclassified. The distinction between tundra, open forest, and barren classes was difficult to identify, especially for areas at higher elevations, such as on the plateaus in the northeastern portion of the study area. Although the lower producer accuracies are indicative of errors related to the classification algorithm, the kappa statistic is greater than 90% for all classification periods, indicating strong agreement between the four classifications and their corresponding validation datasets [52]. Overall, the accuracies of the period-specific land cover maps are sufficient for assessing regional land cover pattern and analyzing changes.
While a total of six land covers were classified, this discussion will focus on the barren, closed forest, tundra/open forest, and water classes. Areas classified as snow are late-lying snowbanks persisting late into the summer season in extremely small portions of the study area and can vary significantly between classified composites. Urban areas are also excluded from the discussion of results, which will focus on natural drivers of change. Future permafrost change studies in the Lower Yenisei River could use these additional classes to examine anthropogenic and snow-related processes.
Dense taiga, or closed forest, abundant in the southern portion of the region, gradually transitions to open forest and tundra moving northward. The majority of the closed forest class is confined to lower elevations of the Yenisei River Valley and in proximity to large water bodies, while the tundra and barren landscapes dominate higher latitudes and elevations. However, localized factors greatly complicate generalized regional vegetation trends. For example, throughout the study area, the closed forest and tundra landscapes are interwoven in a complex pattern. Such spatial heterogeneity is controlled by highly variable permafrost conditions, reflective of a complex history of sedimentation and landscape development [29].
Over the observed 32 years, the total area occupied by closed forests expanded by 38% to occupy an additional approximately 21,000 km 2 . Barrens also expanded by 34% to occupy approximately 2000 km 2 more of the study area. Expansion of closed forests was accompanied by a 25% decrease in the total areal extent of the tundra/open forest class (approximately 22,000 km 2 ) and by a slight 4% reduction (approximately 700 km 2 ) in area occupied by water bodies.
To analyze the spatial pattern of vegetation change between the two-year periods, maps depicting shifts in land cover classes were compiled. Figure 5 shows the total changes (expansion and reduction) in barrens, closed forests, tundra/open forests, and water bodies between the earliest and latest composites classified, 1985-1987 and 2015-2017. The difference maps in Figure 5 indicate that the tundra/open forest and closed forest land cover classes are closely linked through their spatial patterns of change. Throughout the study area, the reduction in tundra/open forest landscapes was accompanied by a corresponding expansion of closed forests. The observed pattern is additionally corroborated by comparing the change in classes on a pixel-by-pixel basis. Using the Combine ArcMap tool, an image was produced with an appended attribute table containing unique "from-to" sequences, which contains the class of pixels in the 1985-1987 image composite classification, and the class that the pixels transition to in the 2015-2017 image composite. Table 3 provides the total area represented by the unique class transition sequences, with trends similar to Figure 5, where the greatest difference is in the transition from open forest to closed forest, totaling 24,857.37 km 2 in the 32-year period.
The change from open to closed forest is by far the most common transition, comprising 70.6% of the total areal extent of change observed between the two classifications, followed by the transition from closed forest to open forest (9.8%) and open forest to barren (5.3%).    Surface hydrology compared to other land cover change remained relatively consistent throughout the 32 years observed, and the linear trend was not significant (Figure 6). To further analyze the hydrologic changes detected, only relatively large lakes were considered. The dense stacking methodology, which utilizes all available imagery to compile a best composite image, can be problematic for evaluating long-term changes in streams and small lakes subject to seasonal fluctuations. For example, a given image composite with the majority of its images from spring and/or early summer might result in artificially inflated areas occupied by water due to snow-melt flooding [26][27][28]. To address this potential problem, the following water bodies were eliminated from the analysis following guidelines provided by Smith et al. [24]: all rivers, streams, and lakes smaller than 40 ha, and lakes connected to river systems, as these are either highly sensitive to seasonality or ephemeral in nature compared to larger closed lakes [53]. The total area occupied by large (>40 ha) closed lakes and the number classified in each three-year composite period are presented in Figure   Figure 6. Bar graphs show total areas occupied by the natural land cover classes of interest separated into continuous (black bars) and discontinuous (gray bars) permafrost zones of the study area. Each bar represents one of the classified three-year composites. Lines represent linear trends in the total areas occupied by each land cover class.
Surface hydrology compared to other land cover change remained relatively consistent throughout the 32 years observed, and the linear trend was not significant (Figure 6). To further analyze the hydrologic changes detected, only relatively large lakes were considered. The dense stacking methodology, which utilizes all available imagery to compile a best composite image, can be problematic for evaluating long-term changes in streams and small lakes subject to seasonal fluctuations. For example, a given image composite with the majority of its images from spring and/or early summer might result in artificially inflated areas occupied by water due to snow-melt flooding [26][27][28]. To address this potential problem, the following water bodies were eliminated from the analysis following guidelines provided by Smith et al. [24]: all rivers, streams, and lakes smaller than 40 ha, and lakes connected to river systems, as these are either highly sensitive to seasonality or ephemeral in nature compared to larger closed lakes [53]. The total area occupied by large (>40 ha) closed lakes and the number classified in each three-year composite period are presented in Figure 7, broken down into continuous and discontinuous permafrost zones within the study area. Over the 32-year period, there were 100 fewer large lakes in the entire study area, and their total area was reduced by 6.6%. Within the continuous permafrost zone, a linear regression shows large lakes increasing by 93.2 ha per decade (R = 0.12 and p-value = 0.73), but decreasing within the discontinuous zone by 6.7 ha per decade (R = 0.25 and p-value = 0.45). There has been significant variation in both of these zones, resulting in weak observed relationships. While the majority of these lakes have expanded, and in some instances new lakes have formed, there are also a number of lakes that have shrunk and/or disappeared in both permafrost zones (Figure 8). However, the degree of variation over the last two decades of the observed period makes definitive spatial patterns of lake area reduction difficult to determine.
Remote Sens. 2018, 10, x FOR PEER REVIEW 12 of 21 7, broken down into continuous and discontinuous permafrost zones within the study area. Over the 32-year period, there were 100 fewer large lakes in the entire study area, and their total area was reduced by 6.6%. Within the continuous permafrost zone, a linear regression shows large lakes increasing by 93.2 ha per decade (R = 0.12 and p-value = 0.73), but decreasing within the discontinuous zone by 6.7 ha per decade (R = 0.25 and p-value = 0.45). There has been significant variation in both of these zones, resulting in weak observed relationships. While the majority of these lakes have expanded, and in some instances new lakes have formed, there are also a number of lakes that have shrunk and/or disappeared in both permafrost zones ( Figure 8). However, the degree of variation over the last two decades of the observed period makes definitive spatial patterns of lake area reduction difficult to determine.

Discussion
A land cover change analysis of the large and diverse region in the Lower Yenisei River basin indicates significant shifts in vegetation over the last 32 years. The most pronounced change detected in this study was the infilling of vegetation in tundra/open forest landscapes and their gradual transition to closed forests. At the 32-year scale, such changes are most likely due to the proliferation of relatively fast-growing shrubs in tundra and open forests. This finding is indirectly supported by numerous studies conducted throughout the circumpolar Arctic (e.g., [7,11,17,[54][55][56][57]). Specifically, within this study area a comparative analysis was done between high-resolution Gambit imagery from the 1960s and contemporary GeoEye-1 imagery over a 58 km 2 area in the vicinity of Dudinka (Figure 1). This comparison across a longer time series revealed a 25.9% increase in shrub cover in a tundra/open forest landscape [8]. Additionally, a transect study in the Putorana Mountains to the northeast of the study area indicated an ongoing filling-in, or densification, of open forests and an upslope shift in tree-line position [49]. Increased air temperatures prolonging the growing season, increased snow accumulation, and a thicker active layer are the primary drivers identified that likely promote these vegetation shifts (e.g., [11,[15][16][17]).
The study area's limited and highly dispersed weather stations preclude a detailed analysis of observed climatic changes at a regional scale. However, World Meteorological Organization station records from Igarka and Norilsk show significant air temperature increases since 1985 [29,40]. These trends, representative of the northern and southern portions of the study area, suggest that the detected regional vegetation change can at least partially be associated with regional climatic warming. However, the complex pattern of land cover change within the study area is likely controlled by highly heterogeneous snow and permafrost conditions. Several studies (e.g., [16,49]) suggest that increasing snow accumulation can be a primary driver of an altitudinal rise in the occurrence of shrubs and treelines. An analysis of 25 × 50 km spatial

Discussion
A land cover change analysis of the large and diverse region in the Lower Yenisei River basin indicates significant shifts in vegetation over the last 32 years. The most pronounced change detected in this study was the infilling of vegetation in tundra/open forest landscapes and their gradual transition to closed forests. At the 32-year scale, such changes are most likely due to the proliferation of relatively fast-growing shrubs in tundra and open forests. This finding is indirectly supported by numerous studies conducted throughout the circumpolar Arctic (e.g., [7,11,17,[54][55][56][57]). Specifically, within this study area a comparative analysis was done between high-resolution Gambit imagery from the 1960s and contemporary GeoEye-1 imagery over a 58 km 2 area in the vicinity of Dudinka (Figure 1). This comparison across a longer time series revealed a 25.9% increase in shrub cover in a tundra/open forest landscape [8]. Additionally, a transect study in the Putorana Mountains to the northeast of the study area indicated an ongoing filling-in, or densification, of open forests and an upslope shift in tree-line position [49]. Increased air temperatures prolonging the growing season, increased snow accumulation, and a thicker active layer are the primary drivers identified that likely promote these vegetation shifts (e.g., [11,[15][16][17]).
The study area's limited and highly dispersed weather stations preclude a detailed analysis of observed climatic changes at a regional scale. However, World Meteorological Organization station records from Igarka and Norilsk show significant air temperature increases since 1985 [29,40]. These trends, representative of the northern and southern portions of the study area, suggest that the detected regional vegetation change can at least partially be associated with regional climatic warming. However, the complex pattern of land cover change within the study area is likely controlled by highly heterogeneous snow and permafrost conditions. Several studies (e.g., [16,49]) suggest that increasing snow accumulation can be a primary driver of an altitudinal rise in the occurrence of shrubs and treelines. An analysis of 25 × 50 km spatial resolution Modern-Era Retrospective analysis for Research and Applications (MERRA) reanalysis data indicates a 30-year positive trend in winter precipitation at higher elevations in the northeastern part of the study area [40]. There, the shifts in land cover from tundra to dense shrubs were detected primarily in stream valleys, which can be attributed to higher snow accumulation occurring in topographic depressions. Similar results were obtained in the air photo-based study of shrub expansion on the North Slope of Alaska [56]. In that study, it was found that in complex topography shrubs grow preferentially in areas that have a greater potential for snow and moisture accumulation, resulting in shrub expansion upslope along drainage channels. A similar pattern of change is evident from the maps in Figure 5. Increasing snow depth can have a profound warming effect on permafrost by offering more thermal insulation to the ground, preventing permafrost aggradation in cold months.
Warming permafrost and increasing active-layer thickness, which can potentially promote shrub expansion, were reported for several regions in Russia (e.g., [29,[58][59][60]). There are no long-term direct permafrost observations on undisturbed natural landscapes within the study area currently available through international open access sources such as the Global Terrestrial Network on Permafrost. However, the Igarka weather station has documented soil temperature at 3.2 m depth increasing from −0.5 • C in the late 1970s to +2.5 • C in 2014, and a corresponding lowering of the permafrost table to about 5 m depth was reported for an artificially grass-seeded surface nearby [29]. Similar observations from the Norilsk weather station revealed a progressive increase in active-layer thickness at a rate of 0.05-0.06 m/year over the 1999-2013 period [29]. Although these changes are probably exaggerated due to anthropogenic influences of increased development and activity in the region, they are indicative of a more general regional permafrost warming trend. The direct active-layer observations at Circumpolar Active Layer Monitoring (CALM) sites located in the tundra/open forest landscapes in the vicinity of Norilsk (CALM site R32) and Igarka (CALM site R40) report a 0.87 cm/year increase in active-layer thickness over the 2005-2017 period in Norilsk and 2.63 cm/year over the 2008-2016 period in Igarka (http://www.gwu.edu/~calm/). This evidence suggests that climate-induced thawing of near-surface permafrost underlying tundra/closed forest landscapes in the study area can partially explain the detected shrub expansion and transition to closed forests. The previously documented relationship between land cover and permafrost within the study suggests that remote sensing-based assessments of vegetation change can potentially be used as an indicator of changing permafrost conditions (e.g., [29,44]).
Changes in permafrost conditions can greatly affect surface hydrology (e.g., [24,[26][27][28]30]). Thaw propagation into ice-rich near-surface permafrost can be accompanied by ground subsidence (e.g., [60]) and/or formation or expansion of thermokarst depressions (e.g., [61]). Lowering of the permafrost table may result in the formation of taliks, enhancing underground drainage and resulting in dryer surface conditions (e.g., [24,25]). To investigate whether the change in number and areal extent of large (>40 ha) closed lakes detected in this study can potentially be attributed to permafrost degradation, associations between lake changes and permafrost conditions were also examined. The lack of a clear regional relation between lake change and permafrost-controlled land cover suggests that changes in climatic factors (e.g., precipitation, evaporation) rather than subsurface drainage might be responsible for detected lake change. While permafrost conditions can play a significant role in the surface hydrology of the study area (e.g., [29]), they are more likely to be manifested through changes in smaller lakes and ponds. More detailed study using high-resolution imagery is required to detect those changes.

Conclusions
The Landsat dense stacking methodology as performed in GEE proved to be a successful means of quantifying land cover change in an arctic environment by effectively minimizing the effects of cloud cover and data gaps by utilizing the entirety of the available Landsat record. We tested and applied the method to analyze 32 years of spatial changes in vegetation and surface hydrology for an area of more than 60,000 km 2 in the Lower Yenisei River region. The goal of this mapping effort was to relate land cover changes to underlying permafrost conditions, with the direct connection between permafrost continuity, depth, and age and overlying vegetation previously established in this region by Tyrtikov [44], then Rodionov et al. [45], and most recently by Streleskiy et al. [29].
The Lower Yenisei is undergoing significant warming and increasing total annual precipitation. Increased air temperatures and snow accumulation together thicken the active layer and therefore are likely the primary drivers of arctic vegetation shifts (e.g., [11,[15][16][17]). Within the study area over the observed 32 years, closed forests and barren ground expanded by 38% (21,000 km 2 ) and 34% (2000 km 2 ), respectively, with a corresponding 25% decrease in tundra/open forests (22,000 km 2 ). Water bodies also slightly decreased by 4% (approximately 700 km 2 ). The significant expansion of closed forests indicates a thickening of the active layer and, specifically, degradation of large areas underlain by younger Holocene-age, near-surface permafrost, likely from the infilling of large shrubs and trees in tundra/open forest landscapes (e.g., [8,49]).
Spatial variation in large lake reduction and expansion have also been linked to permafrost degradation (e.g., [24,[26][27][28]30]). However, there was no significant statistical relationship detected between changes in large lake extent and changes in other vegetative cover in the lower Yenisei River region. The thermokarst processes that affect surface hydrology are likely highly localized in this area [29] and are occurring at a finer scale than the 30 m resolution from Landsat imagery. Future work should consider applying the method developed and described here to higher-resolution imagery to examine and quantify finer-scale land-cover changes in order to detect more of these edaphic permafrost conditions.

Conflicts of Interest:
The authors declare no conflict of interest.