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

Abstract

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 cloud- and snow-affected pixels. Climatic changes observed within the study area’s natural environments show a statistically significant steady greening (~21,000 km² transition from tundra to taiga) and a slight decrease (~700 km²) 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.

1. 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 primary factors cited as facilitating the northward advance of shrub species in the Arctic (e.g., [5,13,14,15,16,17,18]). Several studies have also detected shifts in treelines both northward and at higher elevations in arctic and subarctic regions (e.g., [19,20,21,22,23]).

Changes in arctic surface hydrology are also well documented. Regional-scale appearances and disappearances of lakes in permafrost-affected areas have been observed throughout the circumarctic region (e.g., [24,25,26,27,28]). Widespread reductions in lake extent and complete drainage have occurred primarily within the discontinuous permafrost zone, while lake expansion and new lake formation are occurring within the continuous permafrost zone. This spatial pattern can be attributed to climate-induced changes in permafrost conditions: permafrost degradation and increased groundwater storage in the discontinuous permafrost zone (e.g., [29]) and intensification of surface thermokarst processes and ponding in continuous permafrost (e.g., [24,30,31]).

Satellite-based remote sensing is the primary tool for assessing broad-scale land cover change. The Landsat suite of sensors (Thematic Mapper, Enhanced Thematic Mapper Plus, and Operational Land Imager) provides the longest record (1972 to the present) of continuous, consistent, and freely available imagery for the entire globe (e.g., [32]). However, the use of optical satellite imagery is complicated for arctic applications due to issues including consistently high cloud cover, data gaps, and polar night (e.g., [33,34,35]). Similar challenges for remote sensing–based studies of land cover and land use changes in lower-latitude regions were successfully addressed by utilizing dense stacks, or dense time stacks, of Landsat imagery. The dense stacking methodology has been applied to a variety of spatial change studies, including forestry (e.g., [36]) and urban development (e.g., [37]). By ignoring the established use of an anniversary date for change detection, dense stacking utilizes all available scenes, including images normally rejected due to data gaps and/or high amounts of cloud cover. Stacking all available imagery for a time period of interest allows substitution of areas obscured by cloud cover and/or data gaps in a particular image by coverage from another image within the stack. This approach creates “cleaner” coverage within a dataset by compensating for clouds and data gaps that typically introduce error. Additional spatial information such as digital elevation models (DEMs), maps of known features, and/or calculated indices (e.g., NDVI) suitable for the development of a more nuanced classification scheme can be incorporated [37].

The dense stacking methodology is also starting to be applied in arctic environments (e.g., [38,39,40,41]). Applying the dense stacking methodology in Google Earth Engine (GEE) presents an opportunity to overcome the challenges posed by high-latitude environments and significantly increase the data capacity of Landsat-based arctic impact studies. In this study, we applied the Landsat dense stacking technique to the detection and characterization of land cover change across a relatively large area of the Siberian Arctic. The objectives of this study were: (1) to test the applicability of the dense stacking methodology for arctic land cover change assessments using GEE, and (2) to analyze more than 30 years of previously undocumented spatial vegetation and surface hydrology changes in the Lower Yenisei River Region of central Siberia.

2. Data and Methods

2.1. Study Area

The study area in the central Siberian Arctic was defined by 2 overlapping Landsat scene extents (path 153, rows 11 and 12) for an approximately 60,750 km² area representative of the lower Yenisei River region. This area is of particular interest as it includes several major arctic cities and is a physiographic, vegetative, and permafrost transition zone undergoing significant climatic warming [29,40]. This major river system contributes the largest annual discharge into the Arctic Ocean [42] and separates the West Siberian Plane and the Central Siberian Plateau, starting at the Putorana Mountains and Plateau included within the study area (Figure 1).

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 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.

2.2. Data Acquisition and Processing

2.2.1. 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. Images comprising each three-year composite.

Satellite	Years	Paths	Rows	No. of Images
Landsat 5	1985–1987	151–156	10–13	146
Landsat 5	1988–1990	151–156	10–13	148
Landsat 5	1991–1993	151–156	10–13	108
Landsat 5	1994–1996	151–156	10–13	122
Landsat 5	1997–1999	151–156	10–13	107
Landsat 7	2000–2002	150–156, 230	10–13, 231	143
Landsat 7	2003–2005	150–156, 226	10–13, 234	116
Landsat 7	2006–2008	150–156	10–13	126
Landsat 7	2009–2011	150–156	10–13	140
Landsat 7	2012–2014	150–156	10–13	200
Landsat 8	2015–2017	150–156, 222–231	10–13, 231–233	285

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).

2.2.2. 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 (1985–2017) 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. Confusion matrices for the four classified image examples shown in Figure 4.

a. 1985–1987	Validation Data
Classification Data	Barren	Closed Forest	Tundra/Open Forest	Water	Snow	User Accuracy
Barren	1146	0	4	0	9	98.9%
Closed Forest	0	4564	0	0	0	100.0%
Tundra/Open Forest	0	14	1048	0	0	98.6%
Water	91	0	87	8109	137	96.2%
Snow	0	0	0	397	362	47.0%
Producer Accuracy	92.6%	99.6%	92.0%	95.3%	71.3%
Overall accuracy = 95.34%, kappa statistic = 92.6%.
b. 1994–1996	Validation Data
Classification Data	Barren	Closed Forest	Tundra/Open Forest	Water	Snow	User Accuracy
Barren	1159	0	0	0	0	100%
Closed Forest	0	4535	29	0	0	99.4%
Tundra/Open Forest	0	33	1029	0	0	96.9%
Water	76	0	0	8348	0	99.1%
Snow	12	0	0	0	747	98.4%
Producer Accuracy	92.9%	99.3%	97.3%	100.0%	100.0%
Overall accuracy = 99.1%, kappa statistic = 98.5%.
c. 2000–2002	Validation Data
Classification Data	Barren	Closed Forest	Tundra/Open Forest	Water	Snow	User Accuracy
Barren	846	0	0	0	313	73.0%
Closed Forest	0	4564	0	0	0	100.0%
Tundra/Open Forest	0	5	1039	0	0	99.5%
Water	0	0	0	8424	0	100.0%
Snow	6	0	0	249	442	63.4%
Producer Accuracy	99.3%	99.9%	100.0%	97.1%	58.5%
Overall accuracy = 96.4%, kappa statistic = 94.2%.
d. 2015–2017	Validation Data
Classification Data	Barren	Closed Forest	Tundra/Open Forest	Water	Snow	User Accuracy
Barren	954	0	0	0	205	82.3%
Closed Forest	0	4564	0	0	0	100.0%
Tundra/Open Forest	0	13	2278	0	0	99.4%
Water	133	0	0	8158	133	96.8%
Snow	0	0	0	7	752	99.1%
Producer Accuracy	87.8%	99.7%	100.0%	99.9%	69.0%
Overall accuracy = 96.8%, kappa statistic = 95.82%.

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 10.5. This postclassification conditional statement allowed for the inclusion and tracking of built-up areas in the final 3-year land cover classifications of the study site (Figure 3, step 8).

3. 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². Barrens also expanded by 34% to occupy approximately 2000 km² 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²) and by a slight 4% reduction (approximately 700 km²) 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. Total area change (km²) observed by land cover class in the study site from the 1985–1987 and 2015–2017 Landsat composites.

		To (2015–2017)
	Area (km2)	Barren	Closed Forest	Tundra/Open Forest	Urban	Water	Snow
From (1985–1987)	Barren	x	57.73	795.10	0	11.93	322.94
	Closed Forest	106.19	x	3 459.89	2.80	236.13	46.63
	Tundra/Open Forest	1865.28	24,856.37	x	20.92	631.89	235.51
	Urban	0	1.84	10.73	x	0.10	25.45
	Water	607.12	88.61	721.83	1.82	x	201.79
	Snow	591.25	16.08	228.28	10.40	36.73	x

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² 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%).

To further examine these changes across the 32-year time period, the total areas occupied by different land cover classes were graphed for each of the 11 classified composites (Figure 6). Linear regressions of these totals, with significance (p) indicated by Spearman’s rank correlation coefficient (R), show that barren areas and closed forests increased approximately 78 km² per decade (R = 0.74 and p-value ≤ 0.01) and 5460 km² per decade (R = 0.55 and p-value = 0.08), respectively. Tundra/open forests and water bodies decreased by 6,120 km² per decade (R = 0.62 and p-value = 0.04) and 200 km² per decade (R = 0.09 and p-value = 0.79), respectively. Trends in barrens, closed forests, and tundra/open forests displayed moderate to strong significance, with the most pronounced shifts in these three classes occurring over the first two decades observed (1985–2005) (Figure 6).

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.

4. 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² 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., [44,29]).

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.

5. 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² 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²) and 34% (2000 km²), respectively, with a corresponding 25% decrease in tundra/open forests (22,000 km²). Water bodies also slightly decreased by 4% (approximately 700 km²). 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.