Varying climate response across the tundra, forest-tundra and boreal forest biomes in northern West Siberia

Satellite studies using the normalized difference vegetation index (NDVI) have revealed changes in northern Eurasian vegetation productivity in recent decades, including greening in tundra and browning in the boreal forests. However, apparent NDVI changes and relationships to climate depend on the temporal and spatial sampling and the biome and forest-land cover type studied. Here we perform a consistent analysis of NDVI and climate across four bioclimatic zones (tundra, forest-tundra, northern and middle taiga) in northern West Siberia (NWS), further stratified into eight forest-land cover types. We utilize NDVI data from the Moderate Resolution Imaging Spectroradiometer and climate reanalysis data from 2000 to 2016, a period including the record warm anomaly in 2016 (+2 °C–5 °C June–July surface air temperature (SAT) across NWS). Statistically significant (α = 0.05) correlations were found for two bivariate relationships at the biome level: between NDVImax and June–July surface air temperature (SAT)(r ∼ +0.79), and between middle taiga NDVImax and July precipitation (r ∼ +0.48). No significant statistical relationships were found for the northern taiga and forest-tundra biomes. However, within these biomes we found that deciduous needle-leaf (larch) NDVImax is significantly correlated with July temperature (r ∼ +0.48). Qualitatively, spatial composites of NDVI and climate variables were effective for revealing insights and patterns of these relationships at the sub-regional scale. The spatial heterogeneity of NDVI patterns indicates divergent reactions of specific types of vegetation, as well as local effects that are clearly important on the background of a regional climate response.


Introduction
Satellite-based studies of the normalized difference vegetation index (NDVI) data have revealed changes in vegetation productivity in northern Eurasia and in subarctic and arctic North America in recent decades. A widespread 'greening' indicated by the increasing growing-season maxima of NDVI has been reported in the tundra (e.g. Walker et al 2009, Epstein et al 2012, Bhatt et al 2017, while some recent 'browning' in the boreal forests has been detected (e.g. Beck and Goetz 2011, Elsakov and Teljatnikov 2013, Buermann et al 2014. However the observed trends are highly sensitive to the temporal and spatial sampling and have large spatial heterogeneity . The causes of these changes are not well understood but have been attributed to several factors including regional climate warming. The climate response of northern Eurasian vegetation does not appear to be universal but varies between different periods (Buermann et al 2014), different biomes-e.g. tundra versus middle taiga-and can vary substantially even within a given biome (e.g. Kharuk et al 2013, Hellman et al 2016. Comparable results have been found in pan-Arctic and North American studies (e.g. Bunn et al 2007, Olthof and Latifovic 2007, Epstein 2016, Bhatt et al 2017. Bunn and Goetz (2006) highlight that vegetation changes are highly spatially variable, and that certain species and landscape types are likely more sensitive to changes than others, thus indicating that different vegetation species respond differently to climatic changes (Tsebakova et al 2010, Andreu-Hayles et al 2011. Identifying spatial patterns provides important insights into the controlling factors of vegetation and the potential responses to climate and environmental change, as well as separating climate effects versus disturbance effects (e.g. Sulla-Menashe et al 2018).
Studies of the tundra response to climate in northern Eurasia have been based on local plot-level experiments (Hollister et al 2005, Walker et  Studies have strongly suggested that tundra growth positively reflects temperature change, although Walker et al (2009) found only a weak relationship for tundra on the Yamal peninsula in northern West Siberia (NWS) and Blok et al (2011) found no consistent temperature relationship for tundra in northeast Siberia.
Boreal forest-climate response studies in northern Eurasia have been based on field measurements of tree-ring width (e.g. Lloyd andBunn 2007, Hellman et al 2016) and satellite studies using NDVI (e.g. Goetz et al 2011, Buermann et al 2014. The adverse response of the boreal forest under recent warming in western central Eurasia found by Buermann et al (2014) appears consistent with the findings from recent studies elsewhere that suggest that this biome is becoming more vulnerable to warming-related factors including temperature induced drought.
There have been essentially no studies of northern Eurasia that have consistently examined NDVI-climate relationships spanning from the tundra to the middle taiga-and across the forest-tundra and northern taiga transitional biomes-except for an overview chapter by Goetz et al (2011). It is therefore challenging to compare and contrast results from the above-mentioned studies of tundra to those of the boreal forests, because of differences in temporal and spatial sampling, data and methods-thus there is some ambiguity assessing the climate response. Further, most previous satellite NDVI studies in northern Eurasia that explored ecosystem responses to climate change have relied on trend analysis or changes between time intervals. Some notable exceptions are, e.g. Olthof and Latifovic (2007), Bauermann et al (2014) and Hellman et al (2016). Of particular relevance here is the study of Buermann et al (2014), who effectively examined interannual covariability in boreal forest NDVI and temperature and precipitation with focus on western central Eurasia (50-65°N, 50-75°E).
Here we examine the interannual variability in NDVI and the relationships to temperature and precipitation across four major bioclimatic zones (tundra, forest-tundra, middle taiga and northern taiga) in NWS. We perform a consistent analysis-same temporal and spatial sampling, datasets and methodsacross these bioclimatic zones and further stratified into eight forest-land cover types. This provides a means to: (a) identify the detailed spatial-temporal variability in NDVImax in recent years including the record warm summer 2016, (b) identify the particular climate variables and month(s) from May through July that influence NDVImax, and (c) explore the varying climate response of vegetation between and within different biomes and forest-land cover types.
2. Study area, data and methods 2.1. Study area and its vegetation NWS lies east of the Urals Mountains and west of the Yenisei River, spanning approximately 55-72°N, 54-80°E (figure 1). The distribution of zonal vegetation within NWS follows general climatic patterns and forms four major bioclimatic zones: tundra and forest-tundra in the north, and northern taiga and middle taiga in the south. NWS has relatively flat terrain with predominantly latitudinal zonation. Different climatic factors promote or limit plant growth and development in every bioclimatic zone in NWS. While average air temperature is clearly important for vegetation across the region, precipitation also plays an important role, as is elucidated here in our study of the varying climate response of the different biomes.
Tundra vegetation develops under limiting conditions of a short growing period, low temperatures (air and soil), and high humidity due to low soil penetration and evaporation. Tundra expands on the continuous permafrost region and keeps it frozen most of the year. The main vegetation types are lichens, bryophyte, graminoids, forbs, and shrubs. Trees cover only 2% of the tundra zone and comprise only 0.5% of the entire forested area of NWS. Larix sibirica together with Siberian spruce (Picea obovata) forms the polar limit of tree-stand distribution (see figure 1).
Forest-tundra is a transitional zone 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 forested areas in NWS. 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). The tundra and forest-tundra are the coldest zones of NWS. The mean temperature of the warmest month (July) is 5°C-10°C. The length of growing season in the northern latitude is short: 50-60 d in the tundra zone and 60-90 d in taiga zone. June, July and August are the months with highest biomass level, with the peak of growing activity in July.
The boreal forests of NWS are comprised of the northern and middle taiga. In the northern taiga, the deciduous needle-leaf (DN) forest dominated by Larix sibirica (larch) and light evergreen (EL) coniferous dominated by Pinus sylvestris are the major forest types. 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.
NDVI in NWS shows highest values around the 60°N zone, and gradually decreases northward toward the low temperature polar region (see Miles and Esau 2016, their figure 2(a)).    Miles and Esau 2016). The NDVI data and methods used here are essentially the same as in Esau et al (2016) and Miles and Esau (2016) and are summarized only briefly here -see those references for details. It is acknowledged that for global vegetation studies, using the Collection 5 (C5) for long-term trend analysis may give different results than the newer collection C6 (e.g. Zhang et al 2017). However the C5 data are considered commensurate for our purposes of studying patterns of yearto-year variability on a regional basis.

NDVI data and methods
MODIS NDVI data for the 2000-2016 summers (June-July-August, JJA) across NWS were collected and processed for this study. We used MODIS NDVI 16 d 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). Five tiles to cover the entire area of interest (total 35 tiles per summer) were downloaded and imported into the ArcGIS geographic information system. 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% to 30%. See Miles and Esau (2016) for further details.
Here, vegetation productivity was characterized by annual maximum NDVI (hereafter NDVImax), in accordance with Esau et al (2016) and Miles and Esau (2016) and several previous satellite 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. NDVImax values characterize the maximum development and represent the peak greenness achieved by vegetation during the growing season and eliminates any seasonal variation in NDVI between different vegetation zones. We selected the maximum NDVI value for each year from 16 d composites for each pixel and compiled a baseline dataset of NDVImax, from which we generated: (1) A map of 17 year mean NDVImax throughout NWS, providing the mean state of the spatial distribution of vegetation productivity, (2) NDVImax maps for each year, providing information on temporal and spatial variability of vegetation productivity, and (3) NDVImax values for each year stratified by biome and forest-land cover type, aggregated into mean values for each category and used as input to statistical analysis. Forest dominant species information was generated from the Proba-V-TerraNorte digital map of forest cover for Russia with 345 m spatial resolution and utilizes the same forest types (Bartalev et al 2014) (figure 1). To avoid data discrepancies, the forest map was resampled to 250 m resolution using nearest neighbor interpolation methods, which is appropriate for discrete data such as a land-use classes, because 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 from the forest map. The tree-line was digitized using same data by delimiting the edge of tree vegetation growth. Further details are available in Bartalev et al (2014).

NDVI variability and trends
Interannual variability of NDVImax aggregated for each biome and forest-land cover type was identified through generating a set of time series of NDVImax for each class. From these we evaluated combined years with highest or lowest mean NDVImax value for a specific biome zone compared with mean NDVImax for the entire study period. The spatial composites were based on the particular years of the most positive and negative anomalies, denoted by green and brown dots in the time series in figure 3(a). Spatial maps of NDVImax anomalies were produced as the difference between the average summer NDVImax over the specified years and the average summer NDVImax for the entire study period (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016). The resulting anomaly composites are useful to characterize the state of the vegetation relative to normal (mean) and in comparison to composites of climate anomalies.
Trends in NDVImax for the period 2000-2016 were estimated using ordinary least squares (OLSs) linear regression for each pixel stack, thus yielding the spatial distribution of NDVImax trends across NWS. To identify pixels with statistically significant trends, we masked out all pixels with an estimated p-value>0.05. Details on statistical significance and robustness are provided in Miles and Esau (2016). Trends were also aggregated into biome and forestland cover type, in order to identify the distribution of positive versus negative trends within each classfigures 2 and S1(c), (d) is available online at stacks.iop. org/ERL/14/075008/mmedia.

Climate data
Climate data from the NCEP/NCAR reanalysis dataset (Kalnay et al 1996, updated) were extracted and plotted using the monthly seasonal composites functionalities provided by the National Ocean and Atmospheric Administration. Environmental Science Research Laboratory (ESRL) Physical Sciences Division (PSD). For each month May to August in each year 2000-2016, we extracted monthly-resolution data on surface air temperature (SAT) and precipitation averaged for the grid cells approximating each of the four biome regions. We further divided regions into the grid cells for the western, central and eastern parts of the forest-tundra, northern taiga and middle taiga regions, and for the western and eastern parts of the tundra region (i.e. Yamal and Gydan, respectively). These climate time series were used in univariate and bivariate statistical analysis and composite analysis.

Statistical and composite analysis
First, bivariate statistics were used to test first-order relationships between NDVImax and climate, by calculating Pearson's correlation coefficients r for NDVImax-SAT and NDVImax-precipitation for each biome. For the climate variables, different months and combinations of months from May through July were tested to see which had the highest correlation with NDVImax for each biome and land cover-forest type. Previous research has demonstrated that time-lags are important when studying the relationships between climate and vegetation characteristics (see e.g. Lin et al 2014, Wu et al 2015, Seddon et al 2016. Statistical significance at α=0.05 for the small sample size n=17 (degrees of freedom=16) required an r> |0.5|. For those bivariate relationships found to be statistically significant, we also calculated the regression slopes for NDVImax as a function of temperature or precipitation as an estimate of the general climate response.
Second, as a supplement to the 'black box' statistical approach of aggregated values, we performed composite analysis for those bivariate relationships found to be statistically significant. Composite analysis is recognized as a straightforward and effective tool to identify conditions observed during specific states of the climate, i.e. mean, and positive and negative anomalies. Spatial composite maps of anomalies in SAT and precipitation covering the entire study area (corner coordinates 54°N-73°N, 54°E-80°E) were generated using the monthly seasonal composites functionalities on the ESRL PSD website. We examined SAT and precipitation for composited years of large positive and negative anomalies of NDVImax for each biome. The composited maps were assessed qualitiatively through comparing the sign, magnitude and spatial patterns of anomalies of NDVImax and the climate variables.

NDVI variability and trends for different biomes and forest-land cover types
The spatial distribution of significant positive and negative trends in NDVImax stratified by major forest-land cover type is shown in figure 2. Greening is most prevalent in the north, unforested (UF) tundra and deciduous needle-leaf (DN, larch) in the foresttundra and adjacent areas. Browning is most prevalent in the middle taiga, especially in the mixed forest classes (M and MBM) and the dark evergreen (ED) in the middle and northern taiga. While greening and browning generally predominate in the north and south, respectively, there is a large degree of spatial heterogeneity ('patchiness') and differences depending on forest type, e.g. patches of browning of deciduous needle-leaf evident in the eastern northern taiga, and larch greening even in a background of browning (or no change) of other forest-land cover types.
Trends aside, it is the variability that is of interest for identifying the relationships of boreal vegetation and climate. The interannual variability of NDVI for the different biomes and forest-land cover classes is shown in figure 3. For each biome, the years of large positive and negative anomalies of NDVI-green and brown dots, respectively, in figure 3(a)-provide the basis for composite analyses to examine spatially the vegetation response to variations in temperature and precipitation.
It is noteworthy that NDVImax values in 2015 and 2016 were high in the taiga ( figure 3(a)) for a variety of forest types, e.g. broadleaf deciduous (BR), mixed broadleaf majority (MBM), light evergreen (EL), as well as deciduous needle-leaf (DN) and unforested areas (UF) in the forest-tundra and tundra regions respectively ( figure 3(b)). These changes are a result of increases in NDVImax in two years 2015 and 2016, which were anomalously warm.
The climate conditions in 2016 were exceptional, with the Arctic average SAT for the year ending September 2016 being by far the highest on record (Richter-Menge et al 2016), with a +2.0°C anomaly (relative to the 1981-2010 baseline) over land north of 60°N. The largest temperature anomaly in the Arctic happened to occur over NWS during summer 2016. Figure 4(a) shows that in June-July SATs over our study area were 2°C-5°C above average, with highest values over the tundra and forest-tundra regions. These regions had negative anomalies in precipitation, especially in the eastern part, whereas the forested regions of NWS had strong positive precipitation anomalies ( figure 4(b)). There is are pronounced positive anomalies in NDVImax spanning nearly the entire NWS, with a few smaller areas with negative anomalies, e.g. northern Yamal and Gydan and easternmost part of the northern taiga region ( figure 4(c)).

NDVI and climate relationships: biome and land cover-forest types
Relationships between biome-averaged NDVImax and climate variables were found to be statistically significant (α=0.05) for just two bivariate pairs: tundra NDVImax-SAT (r∼+0.79, p=0.02) and middle taiga NDVImax-precipitation (r∼+0.49, p=0.05), with insignificant results for the foresttundra and northern taiga transitional biomes. The lack of significant correlations for other biomes does not necessarily mean the lack of a climate response, but presumably reflect limitations of bulk averaging, given spatial heterogeneity of NDVImax arising from different forest-land cover types within each region, local effects, as well as the small sample size (n=17). Notably, within the forest-tundra and northern taiga biomes, we found a significant relationship between the deciduous needle-leaf (larch) NDVImax and SAT (section 3.2.3).
The particular month(s) having the strongest correlations were also identified. June-July had the highest correlation for tundra NDVImax and SAT. Forested areas had varying results. July was found to be a better predictor than June-July for middle taiga NDVImax and precipitation and for deciduous needle-leaf NDVImax and SAT. However, for broadleaf deciduous and light evergreen forests, correlations with precipitation were highest for June-July and May-June-July respectively, albeit these correlations were below the significance level.

Tundra and temperature
Tundra NDVImax increased from 0.63 to 0.65 (+0.013 per decade, an increase of ∼3.1% per decade) during the period, while the June-July (JJ) SAT increased from 6.0°C to 8.7°C (+1.6°C per decade) estimated by the OLS linear trend. The statistical relationship between tundra NDVImax and SAT for 2000-2016 is shown in figure 5(a). The correlation coefficient r between tundra NDVImax-SAT is +0.79 (p=0.02). The linear relationship is even stronger (r=+0.86) and the slope (response) is steeper if excluding the 2016 extreme temperature anomaly ( figure 5(b)). The climate conditions in summer 2016 shown above in figures 4(a), (b) indicate a ∼4°C-5°C positive SAT anomaly in the tundra region, together with negative anomalies in precipitation, especially in the eastern part (Gydan). This resulted in lower NDVImax than would be expected from the 4°C temperature anomaly, thereby suggesting an effect of the negative anomaly of precipitation, and as supported by independent observations of dry vegetation and fire conditions in the region during summer 2016.
The patterns of SAT and precipitation anomalies during summers of tundra greening and browning are shown in figure 6. Tundra greening and browning patterns are generally associated with positive and negative anomalies in SAT, respectively. Here the composite of the greenest years (figures 6(a)-(c)) shows large positive temperature anomalies of 2°C-3.5°C together with moderately above-average precipitation (0.2-0.7 mm d -1 or up to 4.2 cm accumulated over JJ). Note that the tundra greening observed in these years is negatively correlated with widespread browning observed in the northern taiga and middle taiga. The brownest tundra years (figures 6(d), (e)) correspond to moderately negative temperature anomalies (0.2°C to −0.8°C)  above-average precipitation (0-0.5 mm d -1 , or 0-3 cm for JJA).

Middle taiga and precipitation
The correlation between middle taiga NDVImax and July precipitation for 2000-2016 is moderate (r=+0.48, p=0.05), and the scatterplot shows a generally linear relationship ( figure 7). The relationships between middle taiga NDVI and climate variables are examined further in the spatial composite analysis ( figure 8).
Greening in the middle taiga (figures 8(a)-(c)) is associated with moderate precipitation anomalies across the region, +0.6 to 1 mm d −1 or +1.8 to 3 cm for July during these years. The areas of observed greening (west and central middle taiga) are in qualitative agreement with the precipitation anomaly composite. The middle taiga forest types had increased NDVI in 2015 and 2016 (figure 2(a)), both years with warm anomalies but with positive anomalies in precipitation in the region, in contrast to the negative precipitation anomalies in the tundra region in 2016 (see figures 4(a), (b)) suggesting that strongly positive SAT anomalies (e.g. 2016) do not stress the vegetation if there is above normal precipitation in the forest regions, (i.e. adequate moisture).
The browning composite for the middle taiga (2010, 2011 and 2012) is shown in figures 8(d)-(f) suggests an association with reduced precipitation (-1 to -1.8 mm d −1 , i.e. 3-5.4 cm for July) in the central and west middle taiga, coinciding with above average SAT. 1.0°C-1.4°C. Even at sub-regional scales, the spatial patterns of precipitation and NDVImax anomalies match each other for both the browning and greening composites.

Deciduous needle-leaf forest and temperature
Although neither the northern taiga nor forest-tundra biomes were found to have a general relationship to climate, within these two regions, one class of vegetation had a statistically significant relationship. The correlation coefficient r between deciduous needle-leaf NDVImax-SAT in July is +0.53 (p=0.03) (figure 9). Not only was July the month with the highest correlation, but no other months or combination of months were significant (e.g. June-July). The strength of the correlation was essentially identical for SATs calculated for the northern taiga and foresttundra regions, and is qualitatively consistent with field-based studies in central Siberia (Kharuk et al 2015).

Discussion
This study of a 17 year record (2000-2016) of maximum NDVI (NDVImax) calculated from MODIS data has identified contrasting variability and trends in vegetation between tundra and forest areas, and notably also within-biome differences. Contrasting results highlight the large spatial-temporal variability in vegetation growth responses to climate change in northern areas. The response is determined by the limiting factors of plant growth inherent for particular biome.
We find the strongest general relationship is a positive response of tundra vegetation productivity to SAT, as found in the biome-averaged correlation (r∼+0.79) and composite analysis. Tundra greening and browning patterns appear to be associated with positive and negative anomalies in SAT, respectively, as expected from the bulk-averaged correlation and previous analyses of tundra vegetation reaction to changes in temperature (e.g. Temperature is an important factor that determines the plants growth in tundra. However, the soil moisture controls and strengthens the direction of plant reaction (Bjorkman et al 2018). The wetness of the soil is controlled by the microtopography of the land much more than by annual precipitation because permafrost prevents rain and meltwater from moving into the ground (Grosse et al 2016). In addition, evaporation is slowed by low temperatures. In the case of 2016, anomalies in increased SAT must have increased evaporation and in combination with lower precipitation resulted in evaporation exceeding precipitation and soil moisture. This is consistent with the Ackerman et al (2018) statement that 'shrubs exhibited diminishing marginal growth gains in response to increasing temperatures, indicative of alternative growth limiting mechanisms in particularly warm years, such as temperature-induced moisture limitation.' The other reasons can affect the NDVI during warmer periods is that increasing summer warmth could lead to melting of ground ice and/or thawing of permafrost, the development of thermokarst, and the ponding of surface water, which would decrease the landscape NDVI signal (Reichle et al 2018). In 2016, the depth of seasonal thawing of permafrost was at least 10 cm more than in the past (www.arctic.noaa.gov/Report-Card/Report-Card-2017/ArtMID/7798/ArticleID/694/Terrestrial-Permafrost,earthchronicles.com/natural-catastrophe/ heat-waves-melting-the-permafrost-on-the-yamalpen insula.html). NDVI declines in northern West and Central Siberia can also be the result of disturbances related to warming, such as fire (e.g. Goetz et al 2005, Cuevas-Gonzalez et al 2009, Jin et al 2012, Liebman et al 2015. Warmer and drier summers make vegetation more combustible, and higher temperatures bring about more thunderstorms and lightning strikes that can spark wildfires, such as observed in the tundra-forest in NWS in 2016. On the other hand, some field studies showed no agreement and suggest that interannual variability in tundra productivity is generally constrained and not strongly linked to summer temperatures (e.g. Hudson andHenry 2010, Bjorkman et al 2018), although others showed positive relationships between summer temperatures and tundra plant growth (Callaghan et al 1989, van der Wal and Stien 2014, cited from Reichle et al 2018. Here the composite analysis of the spatial variability of NDVImax during tundra greening and browning (figure 6) shows a degree of agreement between SAT and NDVI anomaly patterns in Yamal, especially during the strong  greening in 2012 and 2013, when precipitation and soil moisture were adequate for enhanced growth.
The middle taiga NDVImax shows an apparent positive response to precipitation, in both the greening and browning composites. Some insight and illustration may be from the composites of green tundra under high temperatures in 2012 and 2013 (figures 6(a), (c)), which were years of boreal forest browning under high temperatures, which in combination with the precipitation deficits gave clear browning (figures 8(d)-(f)). In contrast, the 2016 high temperatures were associated with relatively brown tundra and relatively green forest areas. The difference is that precipitation patterns were different across NWS, essentially mirror images of each other (figures 6(a), (b) vis-à-vis figures 8(a), (b)). This is evident in inferred moisture deficits in the forest in 2012 and 2013 and a moisture surplus in 2016 in the forest and vice versa for the tundra region. In the boreal forest biomes, we find no indications of direct temperature stress that caused browning; i.e. positive SAT anomalies can be associated with a positive NDVI response, if precipitation is also above average, e.g. 2015 and 2016. Therefore, in order to predict changes in the boreal forest productivity it is important for future scenarios to have summer precipitation properly modeled.

Conclusions
In this study, the recent variability in summer maximum NDVI values from MODIS data have been consistently analyzed across four biomes for NWS. Knowledge on the varying response of arctic and subarctic vegetation types to climate has been nuanced here, with the following conclusions: 1. NDVImax in NWS has high interannual variability across and within different biomes and forest-land cover types, which has a large influence on apparent trends in decadal-scale time series.
2. Interannual variability in NDVImax appears to reflect the short-term response of vegetation to temperature and/or precipitation during the growing season. The climate response itself is varying between biomes in NWS, e.g. tundra NDVImax is most closely related to June-July temperature, whereas middle taiga NDVImax is most closely related to July precipitation.
3. There is large spatial heterogeneity of NDVI patterns within biomes such that forest type differences are evident. In particular, deciduous needle-leaf NDVImax is positively correlated with July temperatures in the northern taiga and forest-tundra biomes, whereas other forest types in the taiga do not have a clear response temperature.
The results here represent a first-order characterization of the short-term climate response of vegetation in different biomes in NWS. While some climate response is evident from this analysis, the unexplained variability is far larger than that explained. There is the need for further refined analysis, including multivariate analysis, lagged response effects (e.g. snow cover preceding the growing season) and more specific spatial analysis of NDVI and climate variables. However, this remains challenging because of the lack of climate and environmental data of the same accuracy and spatial precision as the NDVI data. Further, local effects-both natural and anthropogenic-are clearly important in explaining more of the spatial and temporal variability than achieved here.