Spatial heterogeneity of greening and browning between and within bioclimatic zones in northern West Siberia

The climate of the Arctic is changing, and there is strong interest in understanding how vegetation will respond to and contribute to this change. The ongoing changes in plant productivity will affect many aspects of northern systems including changes in active-layer depth, permafrost, biodiversity, wildlife and human use of the region (Bhatt et al 2010). However, the direct impact of the climate warming factors on boreal and arctic vegetation cannot be established unambiguously.

Satellite observations of vegetation dynamics have shown that this warming leads to increased vegetation growth and a marked 'greening' trend in non-forested tundra (Walker et al 2009). Reflecting temperature change, a widespread 'greening' indicated by the increasing growing-season maxima of normalized difference vegetation index (NDVI) has been reported north of 64°N (e.g., Walker et al 2009, Epstein et al 2012). However, since 2003, the greening in tundra has slowed down (Bhatt et al 2013) and the surface temperature trends have become negative or at least more spatially fragmented in the tundra zone (Comiso and Hall 2014).

The region south of 64°N in the boreal forests (taiga) experienced much less greening. In contrast, a widespread 'browning' indicated by the decreasing maxima of NDVI (hereafter NDVImax) was found, in particular, in the most biologically productive forest ecosystems (Beck and Goetz 2011). The browning seems to be a relatively recent and northward advancing tendency, as earlier studies found the NDVImax increased in 1981–1999 in the boreal forest zone (Zhou et al 2001). The reports of extensive tree growth decline or mortality in northern mid- and high- latitudes (Bunn et al 2007) highlight the large spatial–temporal variability in tree growth responses to climate change in northern areas. It indicates that the response to warming varies substantially between different tree species (Goetz et al 2011). Elevated temperatures enhance growth in deciduous species more than in evergreen, and the main boreal forest may respond with a loss of evergreen trees and a shift toward deciduous trees (Way and Oren 2010). At the same time climate warming induces species- and forest-type-shifts and northward migration. Siberian forests may collapse in south areas and become greener in the north (Bonan 2008). The projection is that ecosystems unique to larch forests occurring on permafrost will move farther north, and Siberian pine, spruce, and fir will expand northward into traditional larch habitat (Tchebakova et al 2009). It is very probable that this process will be driven substantially by changes in permafrost regimes. Changing species distribution and biodiversity will have important biological, ecological and social consequences. Tree species play a key role, providing habitat, food or mutualisms with many animals, fungi, microorganisms, and other plants in addition to other ecosystem services and resources for human use. Forests also contain around three quarters of the Earth's terrestrial biomass and thus are tightly linked with atmospheric carbon budgets (Aitken et al 2008).

This study utilizes NDVI data products from the Moderate Resolution Imaging Spectroradiometer (MODIS) onboard of the Earth Observing System-Terra platform satellite. Several previously published studies by Frey and Smith (2007), Epstein et al (2012), Macias-Fauria et al (2012) and others demonstrated that MODIS NDVI can be successfully used in vegetation productivity assessment for Siberian biomes. We used MODIS NDVI 16-day composites with 250 m spatial resolution (MOD13Q1). This product is widely used in phenological and vegetation dynamics studies in the Arctic and boreal zone (e.g., Blok et al 2011). In this paper, we update, detail and expand upon previous satellite studies of high northern latitude vegetation productivity. Relevant previous studies have had either a broader scope (e.g., Zhou et al 2001, Bunn and Goetz 2006, Beck and Goetz 2011, Barichivich et al 2014) or else covered different regions of northern Eurasia (e.g., Elsakov and Teljatnikov 2013) including arctic tundra sub-regions of West Siberia (e.g., Walker et al 2009, Frost et al 2014). Here we focus on the entire NWS spanning four bioclimatic zones: tundra, forest-tundra, northern taiga and middle taiga. Two novel aspects should be mentioned: (1) whereas previous studies analyzed coarse-resolution data, which is likely to exaggerate the extent and magnitude of NDVI trends, we use moderate-resolution (250 m) data; and (2) moderate-resolution data give an opportunity to trace the changes within the same bioclimatic zone, and help to reveal differences between forest types and even down to the species level. The objectives are to: (1) identify the recent (2000–2014) spatial–temporal patterns of variability and trends of NDVImax across NWS in detail, (2) examine the variability of trends in NDVImax for different bioclimatic zones, and (3) examine within-group variability of trends in NDVImax. We assess the distribution of different forest types and species within the same bioclimatic zone in relation to the observed increases or decreases of productivity. In particular, we compare and contrast trends between main forest forming species in NWS: evergreen coniferous forest, with the dark (Picea abie, Picea obovata) or light (Pinus sylvestris) evergreen and evergreen-majority mixed forests, and deciduous needle-leaf forest dominated by larch (Larix sibirica) species.

2.1. Study area and its vegetation

West Siberia is the territory east of the Ural Mountains but west of the Yenisei River. The West Siberian landscape is nearly flat, with large winding rivers that flow thousands of kilometers to the Arctic marginal seas. The Western Siberian plain is the most bogged region of the world (in total 10⁶ km²); in some parts, up to 70%–80% is covered by mires (Kirpotin et al 2009). Forest covers only 36% of northern West Siberia (hereafter NWS) (Sedykh 1997). The length of growing season in the northern latitude is short: 50–60 days in the tundra zone and 60–90 days in taiga zone. June, July and August are the months with highest biomass level, with the peak of growing activity in July. Vegetation zones across NWS are delimited in terms of phytogeographic zone or bioclimatic (latitudinal) zone, according to Sedelnikov et al (2011). The four major bioclimatic zones within NWS, from north to south are: tundra, forest-tundra, northern taiga and middle taiga (figure 1).

Zoom In Zoom Out Reset image size

In the tundra zones of northern NWS, vegetation develops under limiting conditions of a short growing period, low air and soil temperatures, and high humidity. Tundra expands on the continuous permafrost region and keeps it frozen most times of the year. The permafrost active layer here is approximately 25–100 cm (Ukrainceva et al (2011)). The main vegetation types are lichens and bryophyte and also graminoids, forbs, and shrubs. Larix sibirica together with Siberian spruce (Picea obovata) forms the polar limit of tree-stand distribution (see figure 1) over the territory across NWS between the Ob and Yenisei Rivers. Trees cover only 2% of the tundra zone and comprise only 0.5% of the entire area of NWS.

The forest-tundra is a transition zone or ecotone that expands on sporadic permafrost. It has low and open vegetation canopy and forest covers slightly more than 30%, comprising 2.6% of the total NWS area. Evergreen forests here are often swampy, with abundant dwarf shrubs and green mosses dominated by spruce (Picea obovata) and pine (Pinus sylvestris). Better-drained localities are occupied by larch forest (Larix sibirica).

In the northern taiga, the deciduous needle-leaf forest dominated by Larix sibirica and light evergreen coniferous dominated by Pinus sylvestris are the major forest type. Forests cover 52% of northern taiga, which is 16% of the total NWS area. In the middle taiga, evergreen light forest (Pinus sylvestris) is common on the north while in the south deciduous broadleaf birch (Betula) forests, which mainly represent secondary succession, are widely distributed. Forests cover 49% of the middle taiga and it is 17% of the total NWS.

In NWS, forests develop generally in proximity to rivers where drainage is better and, because of the exceptional swamping, they do not occupy large areas but stretch for great distances along rivers. Forests also occur on elevated drier sites in wetlands. However, these forests are sparse and tree growth is impeded (Shahgedanova 2003).

2.2. Data and methods

2.2.1. NDVI data

In this study we examine changes in vegetation productivity in NWS (2000–14) using variations in NDVI, which is a well-established proxy for gross photosynthesis at different spatial scales (Goetz et al 2005) and an index of vegetation greening, density and development. NDVI is defined as a normalized ratio of reflectance factors in the near infrared (NIR) and red spectral radiation bands

$$\begin{eqnarray}&&{\rm{N}}{\rm{D}}{\rm{V}}{\rm{I}}=\,\displaystyle \frac{{\rho }_{{\rm{N}}{\rm{I}}{\rm{R}}}-{\rho }_{{\rm{R}}{\rm{E}}{\rm{D}}}}{{\rho }_{{\rm{N}}{\rm{I}}{\rm{R}}}+{\rho }_{{\rm{R}}{\rm{E}}{\rm{D}}}}\end{eqnarray} \tag{ 1 }$$

where ${\rho }_{{\rm{N}}{\rm{I}}{\rm{R}}}$ and ${\rho }_{{\rm{R}}{\rm{E}}{\rm{D}}}$ are the surface bidirectional reflectance factors for their respective bands. NDVI is based on the contrast between red and NIR reflectance of vegetation, as chlorophyll is a strong absorber of red light, whereas the internal structure of leaves reflects highly in the NIR (Gamon et al 1995). The greater the difference between the reflectance in the red and NIR portions of the spectrum, the more chlorophyll in the vegetation canopy. Vegetation yields positive NDVI values and approaches +1 with increasing plant chlorophyll content or green biomass. Being a synthesizing dimensionless indicator, NDVI has strong correlations with the productivity and biomass of vegetation for various ecosystem types (Gamon et al 1995, Raynolds et al 2008, Epstein et al 2012).

MODIS data for the 2000–2014 summers (June–July–August, JJA) across NWS were collected and processed for this study. Five tiles to cover the entire area of interest (total 35 tiles per summer) were downloaded and imported into the ArcGIS geographic information system (GIS). Images were mosaic and re-projected from original Sinusoidal (SIN) to the universal Transverse Mercator projection (UTM Zone 42N, WGS84 ellipsoid). NDVI data were quality-filtered by the MODIS reliability data provided together with the MOD13Q1 product, to retain only data of the highest quality, excluding snow/ice- and cloud covered pixels (Solano et al 2010).

The data gaps in the raster mosaic pixels were then filled with information using the nearest neighbor statistical interpolation from the surrounding pixels with data. The percentage of excluded pixels is variable, ranging from 10%–30%. (Images with a higher percentage of excluded pixels were eliminated from the analysis, usually one or two images per summer.) Because the spatial patterns of missing pixel values are not random but naturally clustered (e.g., a significant Moran's I test statistic value >0.8 (Moran 1948)), we assessed the potential problems when interpolating over large-connected patches of missing values. We examined the spatial autocorrelation structure of the observed values of NDVI, finding that autocorrelation remains positive and significant to spatial lag 10–20, e.g., ∼0.5 at lag₁₀ in the tundra zone where missing values are most commonly found. This is larger than the typical size of the areas of excluded pixels and the 9 × 9 kernel operator that we used in the interpolation. Moreover, we are focused more on the ecosystem response rather than a single pixel and consider that each pixel-value is a part of surrounding ecosystem values. In our case we believe that interpolation of data is necessary, because by interpolation, each pixel contributes to the general ecosystem response value for each date rather than using no data. This is important for producing the final NDVImax composite.

The NDVI represents a measure of canopy 'greenness' where values below 0.2 are generally non-vegetated surfaces and green vegetation canopies are generally greater than 0.3 (Gamon et al 1995). A 0.3–1 NDVI threshold was used to exclude water, bare soil and other non-vegetated pixel from the analysis.

In accordance with most remote sensing studies of arctic and sub-arctic vegetation (e.g., Raynolds et al 2008, Walker et al 2009, Beck and Goetz 2011, Blok et al 2011), we use annual maximum NDVI (NDVImax). NDVImax values characterize the maximum development and represent the peak greenness achieved by vegetation during the growing season. Summarizing NDVI into NDVImax composites eliminates any seasonal variation in NDVI and reduces the errors in beginning of phenological phases between different vegetation zones (Reidel et al 2005). An NDVImax map for each year (summer, JJA) was compiled by selecting the maximum NDVI value from each 16-day composite for each pixel. This results in a 250 m resolution, 15-year dataset of NDVImax and generated an up-to-date 15-year mean NDVImax map throughout NWS.

2.2.2. Trend analysis

Statistical methods were applied to the time series of NDVImax and derived parameters, in order to provide metrics and to identify and test for changes. Time series characteristics of interest are mean, variability and trend. Temporal linear trends in yearly summer NDVImax over NWS for the period 2000–2014 were estimated using Ordinary Least Squares (OLS) regression for each pixel stack, with year as the independent variable and NDVI as the response variable. The purpose is to do a pixel-wise trend analysis and extract only significant trends over a certain period of time for a selected region for rejecting the null hypotheses. To identify pixels with statistically significant trends, we masked out all pixels with an estimate p-value > 0.05. The slope of the trend line is the annual rate of 'greening' or 'browning' for pixel stacks with significant trends (p < 0.05).

OLS trend analysis remains the most frequently used in environmental studies in general, and in arctic climate change studies in particular. Our OLS trend analysis can be easily compared with previous vegetation trend studies and with the trend analysis of other climate-relevant quantities. The limitation of OLS estimates is that linear trend model is sensitive to unusual data values or outliers, particularly if they are fitted to small data samples (n = 15 in this study). There are a number of robust linear regression methods (e.g., Yu et al 2014b), which are less sensitive to the outliers. The difference becomes particularly significant when the sampled data are suspected to contain heteroscedasticity.

To reduce the above-mentioned uncertainties, we performed a comparison of the OLS regression versus the non-parametric Theil–Sen regression (TSR) method. We found only small systematic OLS-TSR differences for the largest negative and positive trends—the TSR values are smaller than the OLS trend values, as expected. Furthermore, we tested the significance of the difference between trends using the variance statistics for the OLS trends (e.g., Santer et al 2000). We found that the expected and revealed data heteroscedasticity was surprisingly insignificant. As a result, the widely used OLS trends and the robust TSR trends here differ only statistically insignificantly at the 95% confidence level.

To define the similarity/dissimilarities of NDVImax changes between each bioclimatic zone and each forest type, we calculated correlation matrices. Pearson's correlation matrix was used to describe the correlations among variables. These variables were mean NDVImax extracted every year for the area of each zone and each forest type. The correlation coefficients (r) are between −1 and 1 (0: no correlation and ±1: high correlation). A positive (or negative) value indicates that the two correlated parameters increase (or decrease) simultaneously.

We performed a more detailed examination of different land cover and forest type interannual trends in productivity and variation of summer maximum greenness. To distinguish how productivity responses differ between forest dominant species we calculated the proportion of areas browning and greening for each species. The number of pixels that exhibited positive and negative trends was determined for landcover (forest versus non-forest treeless areas) and forest types (forest type refers to the dominant tree species in the overstory of a given site) across biomes in NWS.

Forest dominant species information (figure 1) was generated from Proba-V-TerraNorte digital map of forest cover for Russia with 345 m spatial resolution that utilises the same forest types as Bartalev et al (2014). To avoid discrepancies between data, the forest map was resampled to 250 m resolution using nearest neighbor interpolation methods, which is typically used for discrete data, such as a land-use classification, since it will not change the values of the cells. The maximum spatial error is one-half the cell size.

We calculated the area fractions covered by different forest types from the total area of NWS in each bioclimatic zone (table 1) from Bartalev et al (2014) forest map. The treeline was digitized using same data by delimiting the edge of tree vegetation growth (figure 1)

Table 1.  Proportion of different forest types in different bioclimatic zones in NWS, calculated here using digital dataset from Bartalev et al (2014). Description of forest types by dominant (at least 80%) or subdominant species (at least 20%–60%) of the forest canopy (Bartalev et al 2014): EG—evergreen dark needleleaf forest consisting of spruce (Picea), fir (Abies), and Siberian pine (Pinus sibirica); EL—evergreen light needleleaf forest consisting of pine (P. sylvestris); BrD—deciduous broadleaf forest consisting of birch (Betula), aspen (Populus tremula), oak (Quercus), linden (Tilia), ash (Fraxinus), maple (Acer), and some other deciduous broadleaf tree species; MBM—mixed broadleaf majority forest consisting of deciduous broadleaf species and evergreen needleleaf species; M—mixed forest—proportions of the evergreen needleleaf and the deciduous broadleaf tree species are approximately equal; MNM—mixed needleleaf majority forest consisting of the evergreen needleleaf tree species and the deciduous broadleaf tree species; DN—deciduous needleleaf forest consisting of larch (Larix); UF—unforested areas—area not covered with trees.

	EG	EL	BrD	MNM	M	MBM	DN	UF
Tundra	0.03%	0.00%	0.00%	0.00%	0.00%	0.00%	0.48%	26.32%
Forest-Tundra	0.46%	0.00%	0.00%	0.00%	0.00%	0.00%	2.20%	5.60%
Northern Taiga	3.48%	4.11%	0.94%	0.75%	0.14%	0.06%	6.67%	14.91%
Middle Taiga	1.90%	8.23%	2.96%	2.09%	0.54%	0.14%	0.68%	17.32%
Total per class	5.86%	12.35%	3.90%	2.84%	0.67%	0.20%	10.02%	64.15%

Note: total sums to 100%.

3.1. Patterns and trends in NDVImax within biomes

Figure 2 provide the region-wide overview of the NDVI pattern and its statistically significant trends for the past 15 years. The mean NDVImax (figure 2(a)) in the NWS generally decays from the southwest to the northeast of the territory. The tundra naturally has the lowest NDVImax values. The central, swamped part of the NWS is characterized by a much lower NDVImax, and the largest NDVImax found for the most productive vegetation along the Ob River and between the Ob River and Ural mountains. Typically, the NDVImax is significantly higher on river terraces with better-drained and therefore warmer soils.

Zoom In Zoom Out Reset image size

Only 18% of the total NWS areas display significant trends, with 8.4% being an increase and 9.6% being a decrease in productivity. The map of statistically significant (p < 0.05) trends in NDVImax from 2001–2014 (figure 2(b)) confirms continuing greening in tundra and forest-tundra biomes. Indeed, more than half of all positive changes in NWS are localized in the tundra zones (table 2). However, this greening is highly fragmented and overlay with patches of browning and to the large degree could be associated with different local disturbances. The most significant areas of greening are found in Taz and southern Gydan peninsulas. In general, almost all greening lies between 65° and 70°N latitude (figure 2(b)). The main negative trends in the tundra ecotones are in the northern part of the Gydan peninsula and eastern part of Yamal peninsula are found along the river, in the coastal zone and in the floodplain landscapes.

Table 2.  Proportion (percentage) of statistically significant trends of NDVImax in terms of the NWS total area, zone area and total trend, all stratified by bioclimatic zone.

	% of the total NWS area1		% of the zone area		% of the trend
	Greening	Browning	Greening	Browning	Greening	Browning
Tundra	18.4%	2.3%	89.0%	11.0%	42.8%	4.0%
Forest-Tundra	7.9%	1.1%	87.6%	12.5%	18.4%	2.0%
Northern Taiga	10.0%	17.7%	36.0%	64.0%	23.2%	31.0%
Middle Taiga	6.7%	36.0%	15.7%	84.3%	15.7%	63.1%

¹NWS areas with statistically significant trend. The values in the right columns (% of the trend) are plotted in figure 2(a). 5.

In contrast, the taiga zones of NWS demonstrate widespread negative trends of NDVImax. The negative trends are prevailing here and concentrating in the southwest, south and southeast parts of the NWS (figure 2(b), table 2).

The highest NDVImax in Tundra on record was in 2012–14 (figure 3). For the transition forest-tundra zone and two taiga sub-zones, a dip in NDVImax occured in 2011, and recovery occurred in 2012. In 2014 the NDVImax declined in tundra zone, however in other zones continued to rise. The correlation is higher between neighboring zones (table 3). Forest tundra and the two taiga sub-zones show better correlation with each other than with the tundra zone. The temporal variability in annual NDVImax for the tundra is different than in other zones, suggesting that the causes of changes in vegetation are different here.

Table 3.  Correlation matrix (r) between mean annual NDVImax in each bioclimatic zone in NWS from 2000–14, n = 15.

	Forest-Tundra	Northern Taiga	Middle Taiga
Tundra	0.43	0.09	−0.90
Forest-Tundra		0.68	0.19
Northern Taiga			0.76

3.2. Trends in NDVImax within forest types

Both positive and negative trends occurred in forested and non-forested land of NWS, but declines are dominating for forested areas. The majority of forest types exhibit more negative than positive trends. We observe two main tendencies: increasing share of areas with positive trends from south to the north for each forest type and non-forested areas and vice versa: increasing percent of areas with negative trends from north to the south for each forest type and non-forested areas.

Zoom In Zoom Out Reset image size

A significant portion of the area with declines is observed in evergreen and mixed forest types (table 4, figure 4). Negative trends indicate decreased productivity in dark coniferous dominated by Picea. In line with this observation we observe expansion of browning in light coniferous dominated by P. sylvestris. Dominant pattern of negative trend was detected in the southwest part of NWS associated with decline in light evergreen forest, while in the south and southeast part it is mostly related to dark evergreen, deciduous broadleaf and mixed forests. Declining spruce and pine productivity appears to be the major driver of the downward trend in the NWS taiga zone (table 4). Nevertheless, decreasing productivity is observed in areas dominated by deciduous species as well. Moreover, the negative trend that prevailed in taiga including a drop in greenness for non-forested areas as well.

Table 4.  Proportion of positive and negative trends in each zone within each forest type.

		EG	EL	BrD	MNM	M	MBM	DN	UF
Negative	Tundra	0.07%						0.46%	3.47%
	Forest-Tundra	1.17%						2.10%	0.90%
	Northern Taiga	38.27%	23.53%	11.24%	24.20%	24.71%	17.51%	21.98%	13.25%
	Middle Taiga	50.43%	59.59%	56.71%	73.42%	71.04%	70.59%	8.81%	26.22%
	Total	89.94%	83.12%	67.96%	97.62%	95.75%	88.10%	33.35%	43.85%
Positive	Tundra	0.04%						3.67%	31.91%
	Forest-Tundra	0.62%						16.02%	10.05%
	Northern Taiga	8.04%	6.12%	4.43%	0.46%	0.18%	2.95%	44.15%	8.14%
	Middle Taiga	1.37%	10.76%	27.61%	1.92%	4.07%	8.95%	2.81%	6.05%
	Total	10.06%	16.88%	32.04%	2.38%	4.25%	11.90%	66.65%	56.15%

Zoom In Zoom Out Reset image size

The only greening observed within the background of browning is for forested areas dominated by Larix. Moreover, Larix stands are found to increase in NDVI in tundra and northern taiga zone (table 4). Forest-tundra and northern taiga are larch-dominated zones and highest larch greening is seen there (table 1). On the southern extent of its bioclimatic range (the middle taiga zone) however, larch show prevail browning. The most interesting observation is that among variant type of forest stands growing close to each other the Larix show positive changes when at the same time the others are showing negative trend.

Changes in NDVImax were heterogeneous among different forest types (figure 5), although some common features are evident. Similar interannual variations were seen among the six forest types (ED, EL, Brd, MNM, M, MBM). Linear regression trends for these types over the analysis period were significantly negative (p < 0.05). For these types, the peak vegetation greenness was relatively low in 2001 and then increased with the highest peak in 2003 and a drop again in 2004. In 2009, a sharp decline was observed for each forest type and in non-forested areas. More notably, larch-dominated forest (DN) is only the forest type with a significant positive linear trend and a negative correlation to other species (see table 5), e.g., DN NDVImax increased in 2004 when other types shows drop in greenness. The strongest positive NDVImax anomalies are observed for DN in 2007, and after a slight decline in 2012 there was an ascending trend. The linear trend for UF over the period analyzed was slight but significantly positive (p < 0.05), with an increase since 2012 after some years of decline.

Zoom In Zoom Out Reset image size

Table 5.  Correlation matrix (r) between mean annual NDVImax for each forest class 2000–14, n = 15.

	EL	BrD	MNM	M	MBM	DN	UF
ED	0.94	0.81	0.90	0.84	0.74	0.07	0.47
EL		0.89	0.88	0.86	0.75	−0.00	0.46
BrD			0.75	0.76	0.81	−0.10	0.35
MNM				0.96	0.77	−0.15	0.20
M					0.78	−0.28	0.21
MBM						−0.37	0.23
DN							0.15

This study has shown contrasting trends between forest and treeless areas, and notably also within-biome differences down to the species level. Contrasts in NDVI tendency may depend on different factors controlling phytomass productivity. Despite the recent summer cooling (Bhatt et al 2013) the positive NDVImax trends in tundra could be still linked to warmer temperature. The dip in 'greenness' in 2009 was circum-Arctic and according to Walker et al (2011) was a response to generally cooler summer temperatures across the Arctic due to elevated atmospheric aerosols, including volcanic dust. Temporary dip in 'greenness' in 2009, followed by a circumpolar recovery in 2010, when the mean NDVI for North America and the Northern Hemisphere overall was the greatest on record. In the latest report tundra greenness, derived from remote sensing data, has been declining consistently for the past 2–4 years throughout the Arctic. Elsakov and Teljatnikov (2013) showed a dip in NDVImax for Gydan in 2006 and West Yamal in 2007. Our record (figure 3) of interannual variability and trends for the tundra zone of NWS shows a slow decline started in 2006, with recovery occurring in 2011. The highest NDVImax on our record was in 2012–14 contrary to the general declining trend of greening (Epstein et al 2015). For the transition forest-tundra zone and two taiga sub-zones, the dip in NDVImax happened in 2011, and recovery occurred in 2012 when summertime surface air temperatures in the NWS were 5 to 15 °C higher than those observed during the previous three decades (Pokrovsky et al 2014).

The bioclimatic (latitudinal) NDVI variations mostly related to changes in temperature, but at the same time latitudinal variations in NDVI mostly related to change in precipitations and very much depend on site factor (elevation, hydrology, soils, plant density, etc). For instance, the temperature stress on Larix sibirica is pronounced only in the southern part of its areal range, whereas accelerated growth and its general northward proliferation along river terraces was reported in the northern part of the area (Tchebakova et al 2009). Gridded climatic data associated with tree ring-width around the circumpolar boreal forest (Lloyd and Bunn 2007) show inverse growth responses to temperature. Browning was mainly concentrated in dark coniferous (Picea) species and occurred more frequently in the warmer parts of species ranges, and the Larix species have a higher-than expected frequency of greening. Declining spruce productivity appears to be the major driver of the decreasing trend in Alaskan boreal forest (Beck and Goetz 2011). Wilking and Juday (2005) found that longitudinal variation in tree growth responses to recent warming depend on site type, tree density and precipitation pattern. Similar studies from four sites along a 30 km north–south transect in Alaska show diverse growth responses between and within study sites (Driscoll et al 2006, D'Arrigo et al 2009). Wilking and Juday (2005) and Beck and Goetz (2011) found that low density stands show higher sustainability either on latitudinal or longitudinal variations; for example tundra and forest-tundra characterized by low tree density show increase in productivity whereas taiga with denser forest showing decrease. Larch forests usually form low-density stands due to shade intolerance. Trends in productivity are thus heavily dependent on cover type and the underlying vegetation density (Beck and Goetz 2011).

In general, most recent boreal forest decline appears connected with increased drought conditions (moisture stress) (Bunn et al 2007, Kharuk et al 2013). In northern West Siberia, however, climatic warming actually results in increases in precipitation, rainfall frequency, accumulated snow depth and a reduction in temperature range (Frey and Smith 2003, Bulygina et al 2011). These factors cause bog development processes to become more active on the flat and poorly-drained surface of the NW Siberian plain (Sedykh 1997, Crawford et al 2003, Moskalenko 2013). A reversal of succession from forest to bog results in a 26% decrease of all above-ground biomass (Moskalenko 2013).

Kirpotin et al 2009 describes contrasting processes occur in the southern and northern parts of the West Siberian Plain. In the taiga zone bogs are expanding and these progressive swamping leads to forest death. However, in the tundra zone the situation is totally different. Bogs there are reducing their area and becoming actively colonized by shrubs and trees (Crawford et al 2003, Kirpotin et al 2009). Similar effect has been discovered during dendroclimatic studies in mid- and northern Sweden: trees growing on bogs in the forest zone are more sensitive to humidity variations than trees in the north (Linderholm et al 2002). The predominance of relatively large deadwood and young pines at ridges testifies the beginning of pine generation change in the tree stand in forest-mire ecosystem of middle taiga in West Siberia (Makhatkov 2010). In Eastern Siberia forest decline is caused by high soil water conditions in a permafrost region due to high amount in precipitations (Iwasaki et al 2010).

The observed contrasting trends indicate that along with climatic changes, site factors (soils, permafrost, hydrology, etc) and plant factors (structure, health, phenology, species composition, etc) determine vegetation vitality, growth and productivity in NWS. Greening is highly fragmented and overlay with patches of browning and to the large degree could be associated with different local disturbances.

The most significant areas of greening are found in this study in Taz and southern Gydan peninsulas. This finding is in good agreement with previously reported plot-scale studies using Landsat images (Frost et al 2014). In general, almost all greening lies between 65° and 70°N latitude (figure 2(b)). The main negative changes observed in this study in the tundra ecotones are in the northern part of the Gydan and eastern part of Yamal peninsulas, also reported by Frost et al (2014). Moreover, vegetation in the area experiences considerable pressure from various disturbance factors. The territory of NWS is an extensive area of oil and gas exploitation and production. Contrasting NDVI trends in the same zone are probably due to a combination and variety of these factors. For example, NDVI patterns on the Yamal, Taz and Gydan are correlated with the temperature and also strongly related to difference in landscape factors, erosion, reindeer yearly grazing and anthropogenic landscape transformation (Walker et al 2009, Forbes et al 2009, 2010, Frost et al 2014, Yu et al 2015, Esau et al 2016). We did not have enough auxiliary data to distinguish the nature of the trend, but this must be done to be able to characterize the vegetation changes due to climate or anthropogenic influence on positive or negative changes.

Consistent trends in the NDVI in the arctic and sub-arctic have been reported regardless of the dataset being analyzed (Goetz et al 2005, Bunn et al 2007, Beck and Goetz 2011, Bhatt et al 2013, Elsakov and Teljatnikov 2013, Kharuk et al 2013, Epstein et al 2015 and etc). However, in other regions, the use of different datasets has led to conflicting findings, potentially due to differences in the processing and corrections applied to the satellite data (Alcaraz-Segura et al 2010). In our study, the trend is generally consistent with previous studies (e.g., Raynolds et al 2008, Walker et al 2009, Beck and Goetz 2011, Blok et al 2011). The robustness of the results is supported by several aspects: (1) we used NDVI, an establish metric to study vegetation change (Walker et al 2009); MODIS NDVI series are widely used and well agreed with trends generated from AVHRR series (Goetz et al 2005); (2) we used higher resolution data versus lower resolution, which can cause misleading trends (Alcaraz-Segura et al 2010, Ju and Masek 2016); (3) our results are in agreement with previous high-resolution studies using Landsat (Frost et al 2014); (4) we used time-tested data processing and correction methods (Verbyla 2008; Bjerke et al 2014); and (5) we used standard OLS trend analysis, which is the most frequently used in environmental studies in general and in arctic climate change studies in particular (Bhatt et al 2013).

Several factors influence the stability of the satellite-derived NDVI, e.g., cloud contamination, atmospheric variability, and bidirectional reflectance (Gutman 1991). The changes in NDVI caused by these factors are seen as undesirable noise in vegetation studies. Thus, compositing methods have been developed to eliminate these effects. Three major composite methods: bidirectional reflectance distribution function composite (BRDF-C), constrained-view angle-maximum value composite (CV-MVC) and maximum value composite (MVC) are applied to generate MOD13Q1 product (Huete et al 1999, van Leeuwen et al 1999, Huete et al 2002). It ensures the high accuracy for MODIS NDVI dataset (Yu et al 2014a). General Accuracy Statement of MOD13 products are well discussed at https://landval.gsfc.nasa.gov/ProductStatus.php?ProductID=MOD13.

It should be noted that Wang et al (2012) reported on sensor degradation having an impact on trend detection in North American boreal and tundra zone NDVI with Collection 5 data from MODIS. The main impacts of gradual blue band (Band 3, 470 nm) degradation on simulated surface reflectance was most pronounced at near-nadir view angles, leading to a small decline (0.001–0.004 yr⁻¹; 5% overall between 2002 and 2010) in NDVI under a range of simulated aerosol conditions and high-latitude surface types. This degradation has been reduced in the MOD13Q1 product by the above-mentioned compositing methods.

In this study, the summer maximum NDVI (NDVImax) values from MODIS data have been aggregated and analyzed for the entire northern West Siberia (NWS). The study provides a detailed analysis of 15 years (2001–14) of NDVI variability and trends in four different bioclimatic zones, for forest and non-forest vegetation and different forest types in NWS. The main conclusions are:

• 1.  
  At the regional level, our analysis shows a general zonal trend in the relative proportion of vegetation responding positively or negatively. The general trend for the NWS is an increase in tundra and decrease in boreal forest zone, which appears to be comparable to the previous vegetation response study in the Arctic (e.g., Raynolds et al 2008, Walker et al 2009, Beck and Goetz 2011, Blok et al 2011). About 18% of the region displays statically significant trend, with just over half negative and the reminder positive trend. Browning of forest and non-forest vegetation has a latitudinal gradient and occurred more frequently in the warmer part of the species ranges; in particular the middle taiga sub-zone shows significant decline in both forest and non-forest vegetation.
• 2.  
  There are contrasting trends between bioclimatic zones. However there is substantial spatial heterogeneity in the NDVI trends, with patches of negative and positive anomalies amongst a background of insignificant change.
• 3.  
  There are also contrasting trends for different species within the same bioclimatic zone. In particular, most negative trends in the taiga appear to be related to a decline in evergreen coniferous forest, with the dark (Picea abie, Picea obovata) or light (Pinus sylvestris) evergreen and evergreen-majority mixed forests showing highest decline in productivity (76% of all negative trends for forested areas) in every zone. In contrast, deciduous needle-leaf forest dominated by larch (Larix sibirica) shows a significant increase in productivity (54% of all positive trends for forested areas), even while neighboring different species show decreasing productivity.

These contrasting results underscore the complexity of the patterns of variability and trends in vegetation productivity. The results show the importance of understanding climatic effects in combination with other factors that are causing the vegetation change. More detailed studies of the effect–response of vegetation change in the region are necessary because it have direct impact on local ecosystem and human well-being here. This suggests the need for spatially and thematically detailed studies to better understand the response of different northern forest types and species to climate warming and other disturbances; higher resolution data, such as Landsat (Ju and Masek 2016) could be used to landcover change detection. Furthermore, the results might be useful to coordinate the human practice in the region.

This study was supported by the Belmont Forum through the Norwegian Research Council grant HIARC: Anthropogenic Heat Islands in the Arctic: Windows to the Future of the Regional Climates, Ecosystems, and Societies (no. 247268). We thank the reviewers for his/her thorough reviews and highly appreciate the comments and suggestions, which significantly contributed to improving the quality of the publication. The authors thank Martin Miles for technical discussions and editing.