Strong shrub expansion in tundra-taiga, tree infilling in taiga and stable tundra in central Chukotka (north-eastern Siberia) between 2000 and 2017

The study area is in central Chukotka, NE Siberia, Russia (figure 1). Large parts of the region are characterised by mountain complexes. Central Chukotka has a continental type of climate: cold winters with average January temperatures of −30 °C, summers with average July temperature of +13 °C, short growing seasons (100 d yr^–1) and low precipitation (200 mm yr^–1; Menne et al 2012a, Menne et al 2012b).

Vegetation data were collected during the Russian-German expedition in summer 2016 (Kruse and Stoof-Leichsenring 2016). Vegetation properties were investigated in four areas along a gradient (figure 2) from treeless mountainous tundra (16-KP-04, Lake Rauchuagytgyn area), via the tundra–taiga transition (16-KP-01, Lake Ilirney area; 16-KP-03, Nutenvut lakes area) to the northern taiga (16-KP-02, Bolshoy Anyuy river area). In total, 57 sites were investigated (Shevtsova et al 2019). The sites cover a large Landsat NDVI gradient (0.3–0.8) and differ in elevation (100–900 m a.s.l.), slope angle (0–54°) and aspect.

Zoom In Zoom Out Reset image size

The 16-KP-04 focus area is in the Anadyr Mountains. It represents mostly mountainous terrain with elevation up to 1500 m a.s.l. and open tundra. The area around 16-KP-01 at Ilirney Lake, represents the tundra–taiga ecotone, which is characterised by diverse vegetation: open tundra, shrublands and forest tundra. The highest elevation is 1000 m a.s.l. 16-KP-03 is characterised mostly by open and forest tundra and mountainous areas (700–1000 m a.s.l) like 16-KP-01. 16-KP-02 belongs to the northern taiga biome with low mountains reaching an elevation of 300–500 m a.s.l. in the east. The river lowlands are covered by wetlands and uplands of forests and open land with visible fire scars.

For every field site, a detailed description of the vegetation was made (Shevtsova et al 2019). We used 52 sites (V01–V16, V18–V24, V26–V51, V54, V57–V58) from 57 investigated field-sites. Some sites were excluded due to their location in mountainous terrain being affected by shadows in the satellite images. During fieldwork we investigated vegetation in circular plots of 30 m in diameter. Within each plot, we sampled five 2 × 2 m subplots, with one subplot placed at the centre of the circle and the other four placed 7.5 m away in each of the cardinal directions. Foliage projective cover of tall (>100 cm) shrubs and trees was estimated for the entire plot, whereas foliage projective cover for all major taxa of the ground layer was estimated in each of the five 2 × 2 m subplots and averaged for each plot. To accommodate the absence of the same species at compared sites and to standardise data, we applied the Hellinger transformation (Legendre and Legendre 2012) to the foliage projective cover data.

Landsat collection 1 Level 2 imagery, processed by the United States Geological Survey (USGS) to surface reflectance with a 30-m spatial pixel ground resolution, provides the best available remote-sensing product covering long-time series at a landscape-scale spatial resolution. From the Landsat archive we extracted cloud-free acquisitions for the four focus areas from an early (years 2000/2001/2002) and late time-slice 2016/2017 (appendix A, table A1). We used only peak-summer acquisitions collected from mid-July to mid-August. We found only one cloudless or partly cloudless peak-summer acquisition from 2002 and one from 2017 for the 16-KP-04 area, one each from 2001 and 2016 for 16-KP-01 and 16-KP-03 and one from 2000 and 2016 for 16-KP-02. Apart from summer acquisitions, we selected cloudless snow-covered springtime acquisitions for the corresponding year and covering the same area as the peak-summer acquisitions.

Due to sensor problems on Landsat-7 since 2003 and the fact that Landsat-8 was launched in 2013 we needed to use both Landsat-7 and Landsat-8 data. To provide high consistency between the spectral bands of different Landsat missions, we applied the standard Landsat-7 to Landsat-8 transformation from Roy et al (2016), but we still observed a sensor-related bias. Therefore, we empirically derived new transformation coefficients. We used concurrent Landsat-7 (12.08.2013 and 14.08.2013) to Landsat-8 (13.08.2013 and 15.08.2013) acquisitions and conducted random sampling of 500 000 points (only land, avoiding noisy shadowed areas) to build linear regression models for each spectral band we used in the processing:

$$\begin{align*}L8ETM & + like - GREEN{\text{ }}surface{\text{ }}reflectance = 1.03 \nonumber\\ & \!\!\!\!\!*{\text{ }}L8{\text{ }}OLI - GREEN{\text{ }}surface{\text{ }}reflectance + 0.001\end{align*}$$

$$\begin{align*}& L8{\text{ }}ETM + like - RED{\text{ }}surface{\text{ }}reflectance \nonumber\\ & = 1.05*{\text{ }}L8{\text{ }}OLI - RED{\text{ }}surface{\text{ }}reflectance + 0.0005\end{align*}$$

$$\begin{align*}& L8{\text{ }}ETM \pm like - NIR{\text{ }}surface{\text{ }}reflectance \nonumber\\ & = 0.91*{\text{ }}L8{\text{ }}OLI - NIR{\text{ }}surface{\text{ }}reflectance–0.009\end{align*}$$

$$\begin{align*}& L8{\text{ }}ETM + like - SWIR1{\text{ }}surface{\text{ }}reflectance \nonumber\\ & = 1.04*{\text{ }}L8{\text{ }}OLI - SWIR1{\text{ }}surface{\text{ }}reflectance + 0.008\end{align*}$$

Where $ETM + - like$ is the transformed Landsat-8 OLI to Landsat-7 ETM +, $L8{\text{ }}OLI$ is the original Landsat-8 OLI, $GREEN$ is the green, $RED$ the red, $NIR$ the near-infrared and $SWIR1$ the first short-wave infrared Landsat surface reflectance [0–1].

To the Landsat-7 and transformed Landsat-8 data we applied the topographical c-correction (Riano et al 2003) to compensate for differences in solar illumination due to topography. We derived solar-geometry related parameters (solar zenith angle and solar azimuth angle) from the Landsat metadata. The slope and aspect angles of slopes were calculated from the ASTER Global Digital Elevation (GDEM, Tachikawa et al 2011) with the same spatial 30 × 30 m pixel ground resolution as Landsat. Additionally, we masked clouds and water bodies in the Landsat spectral-band set using an optimised threshold in the NIR band.

These pre-processed Landsat spectral bands were used to calculate Landsat spectral indices (appendix C, table C1) from summer acquisitions.

The length of the first detrended correspondence analysis axis is equal to 2.55 standard deviation units meaning the field data in our case is distributed on a short linear gradient rather than a unimodal distribution (Borcard et al 2011). A suitable constrained ordination method for short gradients assuming linear species responses along the environmental gradients is redundancy analysis (RDA, Legendre and Legendre 2012). Using RDA, we related foliage projective cover data to the Landsat spectral indices by extraction of pixel values using the field sites' coordinates. The best parameters for the RDA model were chosen according to stepwise selection and analysis of variance.

We built up the classification based on RDA scores using two clustering approaches: k-means (hard classification, appendix D) and fuzzy c-means (soft classification; FCM; Bezdek 1981; appendix D). In both cases of clustering, compositionally similar data were arranged into groups according to similar behaviours of the explanatory Landsat data. In FCM classification we set a threshold of 0.6 and values exceeding this remained unclassified and were marked as uncertain (discussed in appendix D). The optimal number of statistically significant clusters for k-means classification was chosen according to the gap statistics method (Tibshirani et al 2001).

Based on hard classification we predicted cluster numbers for Landsat data from four focus areas for both early and late time slices. We obtained difference maps for each area reflecting changes in vegetation composition over the 15-year period. We interpreted the uncertainties that could appear in the land-cover classification in the study region using FCM results from 2016/2017. We validated the accuracy of the derived land-cover maps of 2016/2017 using field data from 2018 (appendix E).

We observed dense forests, formed by larch (Larix cajanderi Mayr) often with characteristic fire scars in the 16-KP-02 area (Shevtsova et al 2019). Sparse forest tundra stands are common for the 16-KP-01 and 16-KP-03 areas. The sparse larch stands are commonly accompanied by low and dwarf shrubs: Salix spp., Betula nana L., Ledum palustre L., Vaccinium spp., Empetrum nigrum L. Elevational transitions between forest tundra and open tundra are often covered by patches of dwarf pine (Pinus pumila (Pall.) Regel). Open tundra on gentle slopes at lower elevations (540–700 m a.s.l) is graminoid rich and accompanied by dwarf Salix and Vaccinium. Intermediate elevations in mountain areas of 16-KP-01 (700–900 m a.s.l.) and 16-KP-03 (700–800 m a.s.l.) as well as lower elevations in 16-KP-04 (500–600 m a.s.l.) are dominated by Dryas octopetala L. accompanied by a selection of herbs (Poaceae, Fabacea, Astraceae etc) and become barren at higher elevations. Mosses are common everywhere, having highest relative abundance in open tundra.

Landsat spectral indices NDVI, NDWI and NDSI together explained 33% of the variance in the field vegetation data. All three of them explain a significant unique portion of the variance (table 1, figure 3). Adding further indices as explanatory variables to the RDA model did not improve the explained variance. Most of the variance (21%) of the field compositional data is explained by NDVI; NDWI and NDSI, explain around 7% and 5%, respectively. The Landsat spectral indices are distributed in the two RDA-axes space, where NDVI is negatively associated with the 1st RDA axis, NDWI positively with the 1st RDA axis and NDSI positively with the 2nd RDA axis. As axes 1 and 2 carry most information (29%) only these two RDA axes were retained for further analysis.

Zoom In Zoom Out Reset image size

The k-means classification (figure 4) yielded four reasonable land-cover classes named 1) larch closed-canopy forest, 2) forest tundra and shrub tundra, 3) graminoid tundra, 4) prostrate herb tundra and barren areas.

Zoom In Zoom Out Reset image size

The larch closed-canopy forest is characterised by the highest NDVI values ($\overline {0.75}$) and low NDSI (<0.6, figure 5). Despite the high biomass content, larch closed-canopy forests have generally poor biodiversity. The vegetation composition consists mostly of abundant L. cajanderi (20–70%) with moss or grass undergrowth (figure 6).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The forest tundra and shrub tundra class is characterised by higher NDVI values (0.68–0.78), but in contrast to larch closed-canopy forest has a higher NDSI (>0.6). In the landscape it represents dense Alnus viridis ssp. fruticosa (Rupr.) Nyman or Pinus pumila shrub communities and open Larix cajanderi forest (<20% cover) with relatively high dwarf birch (Betula nana) abundance (5–25%), with heathers (Ledum palustre L., Vaccinium vitis-idaea L., V. uliginosum L.) and a lichen–moss ground layer.

The graminoid tundra is characterised by moderate to higher NDVI (0.5–0.68), higher NDWI than the first two classes and a higher NDSI (>0.6). This class represents hummock (Eriophorum vaginatum L.) and non-hummock (Carex spp.) graminoid tundra with moss prevailing in the ground layer.

The prostrate herb tundra and barren areas class is characterised by low NDVI (<0.5), very high NDWI (>0.5) and higher NDSI (>0.6). Prostrate herb tundra is dominated by Dryas octopetala and lichens accompanied by some herbs (Fabaceae, Orobanchaceae, Poaceae, Rosaceae, Saxifragaceae etc) and prostrate Salix spp. shrubs.

In the RDA model, larch closed-canopy forest is characterised by a combination of low axis-1 and axis-2 scores indicating high NDVI and low NDSI. NDSI increases along axis 2, indicating highly reflective surfaces covered by snow. High axis-2 scores in combination with low or moderate axis-1 scores indicate moderate to high NDVI, reflecting graminoid tundra or forest tundra and shrub tundra, respectively. The prostrate herb tundra and barren areas class is characterised by very high axis-1 scores originating from high NDWI and low NDVI caused by low biomass and dry barren land surface.

The validation of land-cover maps showed 66% agreement between class numbers predicted for validation sites in tundra and tundra-taiga from projective cover and from Landsat Indices (appendix E). The highest disagreement (19%) was obtained for prostrate herb tundra and barren areas that were misclassified as graminoid tundra when predicted from Landsat data.

Overall, the land cover in the treeless tundra (16-KP-04) has not significantly changed (<5%) from 2000 to 2017. In contrast, the forest tundra and shrub tundra class shows a major increase in area of up to 20%–25% in the tundra–taiga (16-KP-01, 16-KP-03). Larch closed-canopy forest expanded significantly in the northern taiga (16-KP-02) accompanied by an increase of forest and shrub tundra at higher elevations.

The northernmost area (16-KP-04) is today dominated by prostrate herb tundra and barren areas (59% in 2017) and graminoid tundra (36% in 2017, appendix F, figure F1). The change map (figure 7(a) between 2002 and 2017 indicates only a slight transition of graminoid tundra to forest and shrub tundra (4%).

Zoom In Zoom Out Reset image size

In the typical tundra–taiga transition (16-KP-01), a large part of the area is represented by open and forest tundra: graminoid tundra occupies 92% of the area in 2001 and 68% of the area in 2016; forest tundra and shrub tundra occupies only 3% of the area in 2001, whereas by 2016 it covered 27% (appendix F, figure F2). The change map of this area (figure 7(b)) indicates the expansion of forest and shrub tundra by 21% replacing the graminoid tundra from 2001 to 2016. It occurs mostly in the drain valleys and gentle slopes around Lake Ilirney and smaller lakes. South of Lake Ilirney we observe a change from graminoid tundra to prostrate herb tundra and barren areas due to road construction.

The southernmost location (16-KP-03) represents the tundra–taiga transition as well. In contrast to 16-KP-01, it has more prostrate herb tundra and barren areas (34% in 2001 and 29% in 2016, appendix F, figure F3). The change map (figure 7(c)) shows a transition from graminoid tundra to forest tundra and shrub tundra of 19% between 2001 and 2016. Almost no change in the prostrate herb tundra and barren areas is detected. The larch closed-canopy forest slightly expanded by 4% from 2001 to 2016.

The northern taiga 16-KP-02 area is covered by forest tundra and shrub tundra (39% in 2000 and 69% in 2016, appendix F, figure F4) and larch closed-canopy forest (5% in 2000 and 13% in 2016). The change map (figure 7(d)) from 2000 to 2016 indicates a 9% transition from forest tundra and shrub tundra to larch closed-canopy forest, as well as a 40% transition from graminoid tundra to forest tundra and shrub tundra.

Our approach yielded useful information on land-cover change despite facing several methodological challenges. We focussed on changes of major zonal vegetation types. Accordingly, azonal vegetation types were not investigated during fieldwork including polygonal tundra in lake depressions and floodplain vegetation, which are classified as graminoid tundra or forest and shrub tundra.

The vegetation survey plots were placed within an as homogenous as possible environment to represent broad land-cover classes and allow upscaling. However, heterogeneity in tundra and tundra–taiga landscapes is often higher than the spatial resolution of most of the remote-sensing products (Myers-Smith et al 2011, Frost et al 2014). Consequently, the spectral signal of a Landsat pixel with a spatial resolution of 30 × 30 m is mixed if different vegetation types are present within the pixel. For instance, the spectral signal of a green forest canopy with low red reflectance due to photosynthetic pigment absorption and high multiple NIR scattering and thus high NDVI values can be mixed with the high red and low NIR reflectance of a lichen understorey reducing NDVI for the plot (appendix G, figure G1). Also, Pinus pumila shrubland, a characteristic vegetation type in the investigated region, could hardly be distinguished by Landsat because it is usually found alongside sparse vegetation or bare ground and covers small patches of 25–800 m². For example, dense healthy green P. pumila patches occur in V23 (16-KP-03), but a large percentage of the cover in the site is of white-coloured lichen. This causes a significant decrease in NIR and an increase in red reflectance, which results in lower NDVI values for the 30 × 30 m pixel (appendix G, figure G2).

Landsat-derived NDVI and NDWI selected for building the RDA model are known for their potential to reflect vegetation-related properties. The combination of red and NIR spectral bands used for the NDVI calculation highlights plant biomass and its chlorophyll content (Tucker 1979, Boelman, 2003). NDWI utilises green and NIR reflectance and reflects the water content (Lara, 2018). NDSI has been developed to indicate snow cover (Volovcin, 1976), although it has also been used for identifying dense forests without leaves, for example, in central Alaska (Hall et al 1998). In our data, spring-NDSI values of snow-covered landscapes helped to separate closed-canopy forests that are not covered by snow due to high standing wood from other vegetation types.

Our transformation coefficients allowed a high consistency between Landsat-7 and Landsat-8 spectral bands. A comparison of Landsat-7 data with Landsat-8 data without our transformation or with application of the standard transformation (Roy et al 2016) resulted in inaccuracy due to sensor differences of the same strength as vegetation changes. We observe less variability in the spectral bands of Landsat-7 data than Landsat-8 data due to the much lower radiometric sensitivity of the Landsat-7 sensor (USGS, 2019a, 2019b).

We significantly corrected Landsat data for shadowed areas in our classification by combining two approaches. First, we used the standard way of avoiding topographic effects and variations in the sun illumination angle by employing normalised difference spectral indices instead of single spectral band values or single band ratios (Mróz and Sobieraj 2004). Second, we applied topographical corrections (Riano et al 2003), because in the mountainous parts of Chukotka, the quality of the spectral indices was still affected by steep slopes and very strong shadows, especially for the snow-covered acquisitions.

We validated the derived land-cover maps with new field data from the expedition in 2018 (appendix E). 66% of validation sites were correctly assigned to one of the four land-cover classes. The highest misclassification, namely systematic overestimation, we found for prostrate herb tundra and barren areas over graminoid tundra. However, the derivation of relative changes compensates for the misclassification error.

4.2.1. Tree infilling with no evidence of substantial treeline advance

The expansion of the sole tree species Larix cajanderi is attributed to a change from forest tundra and shrub tundra to larch closed-canopy forest. The most prominent results are obtained for the northern taiga area (16-KP-02). The larch closed-canopy forest increase resulted in an absolute L. cajanderi cover increase from 2000 to 2016 of about 5% (or an area of 15.8 km²) with respect to its abundance in each of the accounted land-cover classes (appendix H, tables H1, H2). The spatial pattern of larch expansion derived from the land-cover classes' differences suggests tree infilling behaviour rather than treeline advance since the larch closed-canopy forest only significantly increased in the northern taiga where it was previously present at the landscape scale. Contrarily, we found no evidence of significant larch forest advance in the tundra–taiga ecotone.

We compared the increase of larch tree cover in our study region to other Siberian subarctic regions where the treeline is formed by Larix spp. (figure 8). Both areas of the tundra–taiga ecotone (16-KP-01 and 16-KP-03) have similar increases to north-western Siberia as explored by Frost and Epstein (2013). Increase of larch cover in northern taiga (16-KP-02) is twice as high as in the tundra–taiga ecotone. Despite these similarities, one needs to be careful in such comparisons since tree cover is increasing nonlinearly and one needs to consider the period of investigation (Frost & Epstein: 1965–2009; this study: 2000–2017) and the different spatial resolutions of the studies.

Zoom In Zoom Out Reset image size

It is likely the ecological niche of L. cajanderi will shift further north with recent warming. However, larch migration to the north might be limited by seed dispersal and reproduction rates (Wieczorek et al 2017a). Our inferences of strong infilling in southern areas and a simultaneously slow treeline advance in northern areas is in accordance with simulation results of an individual-based larch model (Kruse et al 2019a). Climate-driven infilling might be overestimated in areas recently heavily impacted by fires. However, we selected a region in the northern taiga area with no fire dynamics within the investigated period.

4.2.2. Shrub expansion and graminoid decline

Throughout the explored region forest tundra and shrub tundra increased in area from 13% in the early time-slice to 32% in the late time-slice (appendix H, table H1). This increase in vegetation biomass is also captured by the significant greening trend between 2002 and 2008 observed in our study area, interpreted as an increase in growing season NDVI (MODIS and AVHRR, in Guay et al 2014) and increase in peak-summer MODIS NDVI from 2000 to 2018 (appendix B).

We interpret the areal increase as expansion of the forest tundra and shrub tundra into areas formerly covered by graminoid tundra. The forest and shrub tundra class has an average shrub cover of 60%, which is much higher than the shrub cover in all other vegetation classes (figure 9). The main taxa contributing to shrubification are deciduous dwarf shrubs Vaccinium vitis-idaea, Betula nana (strong contribution only in tundra–taiga), Ledum palustre and V. uliginosum (tundra–taiga, northern taiga), Empetrum nigrum (stronger contribution in northern taiga, less strong in tundra–taiga), Cassiope tetragona and Salix spp. (minor contribution in both tundra–taiga and northern taiga) and Alnus viridis ssp. fruticosa (strong contribution only in tundra–taiga, figure 10).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The shrubification goes along with a significant decrease in the graminoid tundra class from an average of 60% in 2001%–41% in 2016 throughout the explored region (appendix H, table H1). Dominant taxa of the graminoid tundra are Poaceae and Cyperaceae, which form 37% of the cover of this class, much higher than for other classes.

Myers-Smith et al (2011) describe how shrub expansion in the tundra–taiga ecotone in recent years can be related to various environmental drivers including climate (earlier snowmelt alongside higher temperatures), natural disturbances (permafrost degradation, fire) and anthropogenic impact. Many studies report tall (mostly Alnus, rarely Pinus pumila), low and dwarf (Salix, Betula) shrub expansion all over the arctic region (Sturm et al 2001, Tape et al 2006, Myers-Smith et al 2011, Frost and Epstein 2013, Frost et al 2014). Most of these studies investigate different and/or longer time periods or focus on local-scale changes (floodplains, river slopes etc). The rates of shrubification in the central Chukotka mountainous area appear to be much lower than in the other studied circum-Arctic regions that represent lowland landscapes. We suspect that most of the shrubification in the tundra–taiga transition area in Chukotka could be related to climate changes, because anthropogenic impact is low except in the 16-KP-01 region, where road construction along Lake Ilirney is indicated as a decrease in vegetation cover (figure 7).

Trends in the annual mean air temperature suggests an increase of 0.25 °C–1°C per decade (1979–2009 years trend) in our study area, as well as a decrease in snow cover duration (1–2 d per decade, 1979–2009 years trend; Liston and Hiemstra 2011). The increase in air temperature has caused warming of permafrost in recent years (Biskaborn et al 2019), which may also support a deepening of the active layer and thus can benefit plants with long root systems. According to Fyodorov-Davydov (2008), active-layer thickness significantly increased in eastern Siberia (Northern Yakutia)) from 1999 to 2006 by about 35–70 cm in typical tundra, 7 cm in the tundra–taiga ecotone and 13 cm in the northern taiga.

Despite modern climate conditions in the tundra–taiga ecotone being favourable for both shrub and forest expansion, at first, fast reproduction rates allow low and dwarf shrubs to rapidly colonise the area. These areas will be subsequently replaced by generally slower expanding forests (cf Kharuk et al; 2006, Montesano et al; 2016, Wieczorek et al 2017b, Kruse et al 2019a).

Overall, our findings of shrub expansion and decline of graminoid cover agree with modelling results. For example, a simulation experiment by Carlson et al (2018) for the subarctic claims an increase in the abundance of deciduous shrubby Salix and a decline of the graminoid Carex under a warming treatment.

4.2.3. Stable treeless tundra