The response of Arctic vegetation to the summer climate: relation between shrub cover, NDVI, surface albedo and temperature

Recently observed Arctic greening trends from normalized difference vegetation index (NDVI) data suggest that shrub growth is increasing in response to increasing summer temperature. An increase in shrub cover is expected to decrease summer albedo and thus positively feed back to climate warming. However, it is unknown how albedo and NDVI are affected by shrub cover and inter-annual variations in the summer climate. Here, we examine the relationship between deciduous shrub fractional cover, NDVI and albedo using field data collected at a tundra site in NE Siberia. Field data showed that NDVI increased and albedo decreased with increasing deciduous shrub cover. We then selected four Arctic tundra study areas and compiled annual growing season maximum NDVI and minimum albedo maps from MODIS satellite data (2000–10) and related these satellite products to tundra vegetation types (shrub, graminoid, barren and wetland tundra) and regional summer temperature. We observed that maximum NDVI was greatest in shrub tundra and that inter-annual variation was negatively related to summer minimum albedo but showed no consistent relationship with summer temperature. Shrub tundra showed higher albedo than wetland and barren tundra in all four study areas. These results suggest that a northwards shift of shrub tundra might not lead to a decrease in summer minimum albedo during the snow-free season when replacing wetland tundra. A fully integrative study is however needed to link results from satellite data with in situ observations across the Arctic to test the effect of increasing shrub cover on summer albedo in different tundra vegetation types.


Introduction
Climate change is expected to have a great impact on the Arctic tundra ecosystem (Kattsov et al 2004). Changes in temperature, precipitation and growing season length are likely to have great effects on future vegetation development in the Arctic, resulting in ecosystem feedbacks that may lead to further climatic changes (Wookey et al 2009). For example, changes in tundra vegetation could influence the local water and energy balance by altering the surface albedo and evapotranspiration (Eugster et al 2000).
As evidence for increasing Arctic vegetation growth, greening trends have been observed in multiple high-latitude regions during the last decades, based on the satellite-inferred Normalized Difference Vegetation Index (NDVI) (Jia et al 2009, Bunn et al 2007. NDVI records are a good measure for plant productivity (Tucker and Sellers 1986) and spatially correspond with variations in summer warmth in the Arctic region (Raynolds et al 2008). Changes in NDVI may indicate a response of Arctic tundra vegetation productivity and/or vegetation composition to increases in temperature. Most of the high-latitude areas with a significant positive trend in NDVI are restricted to tundra areas, in contrast to many boreal areas that show browning over the last decades (Bunn and Goetz 2006.
Within Arctic vegetation types, deciduous shrub species such as dwarf birch (Betula nana) are expected to respond most to changes in temperature (Euskirchen et al 2009. Photographic and satellite data suggest that shrubs are at least partly responsible for the recent greening trends (Tape et al 2006, Jia et al 2004. Recently, dendrochronological studies showed that radial growth of several willow (Salix spp.) shrub species in Siberia is positively related to local peak-summer NDVI and summer temperature (Forbes et al 2010, Blok et al 2011b, providing further support for a linkage between increased shrub growth, tundra greening and increasing air temperature.
Along with an increase in NDVI, shrub expansion could lead to a higher absorption of shortwave solar radiation, thereby decreasing the summer surface albedo and potentially leading to additional local atmospheric heating , Euskirchen et al 2009. Indeed, a negative correlation between shrub cover and summer albedo has been observed in Alaskan shrub tundra (Sturm et al 2005, Beringer et al 2005. However, it remains unclear how spatial and temporal differences in summer albedo are related to summer climate and vegetation composition at a larger spatial scale. Our hypotheses were that summer albedo is negatively related to summer NDVI and that summer NDVI increases and summer albedo decreases with increasing summer temperature and shrub fractional cover. To test these hypotheses, we (i) examined the relationship between deciduous shrub cover, NDVI and albedo using data from a field campaign in a low-Arctic tundra site in NE Siberia and (ii) evaluated recent (2000-10) trends and correlations in satellite-derived peak growing season NDVI and minimum albedo with climate station data of four low-Arctic regions in Yakutia (E-Siberia), Nunavut (Canada), Yamal (W-Siberia) and the North Slope of Alaska.

Materials and methods
2.1. Field measurements 2.1.1. Site description. Field measurements were made in a low-Arctic tundra site in the Kytalyk nature reserve (70 • 49 N, 147 • 28 E) in the Indigirka lowlands in northeast Siberia, Russia. The vegetation at the research site consists of a mixture of graminoids, forbs, mosses and shrubs and is classified as a mixture of vegetation classes G4 (tussocksedge, dwarf-shrub, moss tundra) and S2 (low-shrub tundra) on the Circumpolar Arctic Vegetation Map (CAVM) (Walker et al 2005). Measurements of NDVI, albedo, leaf area index (LAI) and plant species cover were made in 20 experimental 10 m-diameter plots during summer 2008 and 2009. Half of these plots were located in Betula nana-dominated patches in a former thermokarst lake bed. The B. nana-dominated patches in this former lake bed area are alternated by wet sedge depressions. The other half of the plots were located on top of a Pleistocene river terrace, hereafter referred to as ridge plots. Vegetation in the ridge plots consisted of a mixture of Eriophorum vaginatum tussocks, evergreen and deciduous dwarf shrubs. Vegetation height was approximately 10 cm in ridge location and 15 cm in the lakebed plots. Plots were selected pairwise on the basis of similarity in vegetation cover. In one plot of each plot pair (5 lakebed and 5 ridge plots), all Betula nana was removed during summer 2007 (Blok et al 2010). Some regrowth of B. nana had occurred in the following years, albeit significant differences in B. nana fractional cover remained during summer 2008 and 2009, when measurements presented in this paper were made. Two plot types could thus be distinguished: plots with dense (control) and with thin (removal) B. nana fractional cover.

2.1.2.
Spectral reflectance measurements and NDVI calculation. Vegetation spectral reflectance measurements were made during July 2008 using a field spectrometer (ASD Field Spec Classic FR, Analytical Spectral Devices Inc., Boulder, CO). The reflectance measurements represent hemispherical-directional reflectance factors (Schaepman-Strub et al 2006), covering the spectral range of 350-2500 nm. Reflectance was measured from a distance of 1 m above the vegetation (from nadir) in all 20 plots, replicated 20 times per plot. The mean reflectance was calculated for each plot by averaging the 20 replicates. Wavelengths affected by atmospheric water vapor (350-399 nm, 1361-1409 nm, 1801-1959 nm, 2401-2500 nm) were excluded from further analysis. The mean reflectance per plot was spectrally resampled to the MODerate resolution Imaging Spectroradiometer (MODIS) satellite response function using ENVI v4.5 software (ITT Visual Information Solutions, Boulder, CO).
NDVI for each plot was calculated using resampled red (band 1, 620-670 nm) and near infrared (band 2, 841-876 nm) field spectrometer reflectance according to equation (1): in which R NIR = near-infrared reflectance and R red = red reflectance. NDVI values were compared against plant species fractional cover in the plots as recorded by point intercept.

2.1.3.
Fractional cover, albedo and leaf area index measurements. Species fractional cover was measured in all plots using the point intercept method. A description of the method and plant fractional cover data for summer 2008 and 2009 are given in Blok et al (2010) and Blok et al (2011a), respectively.
Albedo in the 10 lakebed plots was measured during July 2009. Net shortwave radiation (SWnet) was measured alternatingly in all plots (Blok et al 2011a). Incoming shortwave radiation (SWin) was recorded continuously by a CM7b shortwave albedometer (Kipp and Zonen B.V., Delft, The Netherlands) mounted on a meteorological tower at 4.5 m height (Van der Molen et al 2007). Midday average albedo (13:00-16:00 h) was calculated from 30 min average SWnet (measured at the experimental plots) and SWin data (measured at the meteorological tower) using equation (2): Leaf area index (LAI) of the vascular vegetation in plots was measured at approximately two cm above the moss layer at 20 random points in each plot on a cloudless day on 31 July 2009 using a SunScan canopy analysis system (SS1, Delta-T Devices, Cambridge, UK), which gives LAI values based on Wood's SunScan equations (SunScan user manual v.2.0).

Statistical analysis of field measurements.
Pearson's correlation coefficients (r ) were calculated between NDVI, albedo, LAI and Betula nana fractional cover, as measured in the experimental plots. Multiple linear regression analysis was used to assess how much of the variation in NDVI, albedo and LAI in the plots is explained by fractional cover of the dominant plant species and plant functional types. A stepwise model input sequence was used, whereby only significant plant species or functional types were added to the model. We used a mixed model analysis to determine differences in surface albedo between control and removal plots in the lakebed. Because surface albedo was measured pairwise in one control and one removal plot at a time, factor 'pair' was included in the model as a random factor to take into account potential differences in surface albedo between plot pairs caused by weather conditions. All data were checked for homogeneity of variance and normal distribution and natural-log transformed where necessary to achieve normal distribution and homogeneity of variance. Note that in graphs the non-transformed raw data are shown for better interpretation but in statistical analysis transformed data were used when necessary to meet assumptions of normality. SPSS v.17.0 was used for all statistical analysis.

Satellite and meteorological data
2.2.1. NDVI and albedo data. We analyzed the relationship of annual growing season maximum NDVI, minimum albedo and meteorological station data in a selected 200 km × 200 km study area in Yakutia, NE Siberia centered at our field site (70 • 49 N, 147 • 28 E), along with three other low-Arctic tundra areas that can be characterized as plain landscapes (CAVM Team 2003). These similar-sized study areas are located on the North Slope of Alaska (center at 70 • 0 N, 154 • 0 W), on the Yamal peninsula, NW-Siberia (center at 69 • 3 N, 70 • 0 E) and in the Queen Maud Gulf Migratory Bird Sanctuary, Nunavut, Canada (center at 67 • 0 N, 101 • 3 W) (figure 1). Hereafter, the study areas are referred to as YAK, ALA, YAM and NUN, respectively. All four areas are located within bioclimatic subzones D (summer warmth index (SWI) 12-20 • C) and E (SWI 20-35 • C) and represent CAVM classes barren tundra, graminoid tundra, shrub tundra and wetlands (figure 1). Substrate soil chemistry is slightly acidic to circumneutral (pH 5.5-7) in mostly all parts of the four areas (Walker et al 2005).
MODIS albedo, NDVI, and land cover data were obtained from the Oak Ridge National Laboratory (ORNL) Distributed Active Archive Center (DAAC) for biogeochemical dynamics using the MODIS subsetting tool (http://daac.ornl.gov/ MODIS/modis.shtml), accessed January, February and July 2011. Sinusoidal projected data were acquired for the four low-Arctic tundra areas (figure 1) for the period 2000-10. We used the MODIS NDVI product (MOD13Q1, 250 m spatial resolution, 16-days composites produced every 16 days) to compile annual maximum NDVI (maxNDVI) maps by selecting the maximum NDVI from the 16-days composites for each pixel. (Holben 1986). NDVI data were quality-filtered by the MODIS subsetting tool, excluding snow-and cloudcovered pixels.
We used the MODIS shortwave white sky albedo (WSA) product (Collection 5 MCD43A3, 500 m spatial resolution, 16-days composites produced every 8 days), quality-filtered by the subsetting tool, thus excluding magnitude bidirectional reflectance distribution function (BRDF) inversion data (i.e. backup algorithm data), restricting the analysis to full BRDF inversion data (highest quality). We constructed annual peak growing season minimum WSA (minWSA) maps by selecting the minimum WSA value of 5 possible albedo images produced between day of year (DOY) 193 and 225 (12th July-13th August) for each pixel. This period corresponded with the period when maxNDVI was reached at the peak of the growing season.
We used the MODIS land water mask product (MOD44W, 250 m spatial resolution) to restrict further analysis to mostly vegetated pixels. The MOD44W product was resampled to match the 500 m spatial resolution of the MCD43A3 albedo product using the nearest neighbor method in ENVI v.4.5. Annual peak growing season maxNDVI and minWSA values were spatially averaged per vegetation class and bioclimatic subzone by projecting the CAVM vegetation and bioclimatic descriptions (http://www.geobotany.uaf.edu/cavm/) onto the maxNDVI and minWSA layers. Only maxNDVI and minWSA data were analyzed that contained a representative number of pixels per vegetation class and bioclimatic subzone. This resulted in the exclusion of barren tundra in ALA and wetland tundra in NUN for further analysis. Temporal average (2000-10) values of maxNDVI and minWSA were calculated per CAVM bioclimatic subzone and vegetation class.

Meteorological station data.
Raw daily temperature data were obtained for all study areas from regional meteorological stations. These data were used instead of gridded satellite-derived temperature to ensure that completely independent data sets were compared against each other. We combined two independent meteorological station records for areas without a centrally located station to approximate the average regional climate. Only data with a full record (no missing days) were used in this study. As a result of this quality criterion, meteorological stations located close to a research area but lacking a full data record were excluded from analysis. Instead, data from the second-closest station that could provide a full temperature data record were applied (table 1). Summer warmth index (SWI) was calculated for each region as the annual sum of mean monthly temperatures above 0 • C.

Statistical analysis of satellite and meteorological
data. Satellite-derived annual maxNDVI and minWSA were analyzed for trends over the available MODIS data period, 2000-10. Regression slopes in maxNDVI, minWSA and monthly air temperature were calculated per study area, Arctic vegetation class, and bioclimatic subzone. Linear mixed models were used to analyze how much of the inter-annual variation in maxNDVI and minWSA could be explained by variations in summer climate and how inter-annual variations in maxNDVI and minWSA were related to one-another. Hereto, data from all regions were combined with region specified as random factor. Pearson correlation coefficients r were calculated between annual peak growing season maxNDVI, minWSA, and meteorological temperature data to analyze relationships between these variables.

Field data: relation of shrub fractional cover, NDVI and albedo
NDVI was positively related to Betula nana fractional cover (r 2 = 0.69, p < 0.001, n = 20 plots), with highest NDVI values being measured in lakebed control plots and lowest NDVI values in lakebed removal plots (figure 2(a)). Differences in NDVI between control and removal plots were greater in the lakebed than in the ridge location because of higher initial differences in B. nana cover before removal. Together, B. nana, graminoid and moss fractional cover explained 82.7% of the variation in NDVI, as shown by multiple regression analysis (data not shown). NDVI values saturated at B. nana fractional cover values above 40%. Control and removal plots in the lakebed area differed significantly in surface albedo ( p < 0.05, n = 5 plot pairs, figure 2(b)), with higher albedo values being measured in plots with a thin B. nana canopy. Overall, we however did not observe a significant relationship between albedo and NDVI in the field ( p > 0.05, data not shown). B. nana fractional cover explained 46.5% of the variation in albedo between plots. Fractional cover of species other than B. nana did not significantly improve the regression model explaining variation in albedo (not shown). LAI was positively related to B. nana fractional cover (r 2 = 0.52, p < 0.001, n = 20 plots, figure 2(c)), but showed no significant relationship with albedo.

NDVI and albedo relations with tundra vegetation class and bioclimatic subzone.
On average, annual maxNDVI values were found to be greatest in areas classified as shrub tundra, closely followed by graminoid tundra. Lowest maxNDVI values were observed in barren tundra (figure 3). In all areas we observed a large difference in annual maxNDVI between vegetation classes, except for Yamal (YAM), where no clear differences were found. Neither did maxNDVI differ much between bioclimatic subzones for the YAM area, in contrast to the two other Arctic regions with both bioclimatic subzones, Yakutia (YAK) and Alaska (ALA). There, we did see large differences in maxNDVI between bioclimatic subzones, with southern tundra areas showing consistently higher maxNDVI values than northern tundra areas. Annual growing season minWSA generally followed a similar pattern in differences between Arctic vegetation classes and bioclimatic subzones as maxNDVI, with the exception of graminoid tundra showing considerably higher minWSA values than the other Arctic vegetation classes for YAK and NUN. Lowest minWSA values were observed in barren and wetland tundra areas (figure 3), which is most likely due to the relative great proportion of water and bare soil in these areas. Unexpectedly, highest minWSA values were observed in the warmest bioclimatic subzones, which also had the highest maxNDVI. A significant negative relationship between inter-annual variation in regional average growing season maxNDVI and minWSA ( p < 0.01) was observed when data from all  regions were analyzed together using a mixed model analysis. However, when data were analyzed separately per region, a significant negative relationship between maxNDVI and minWSA was observed only for YAM ( figure 4, table 2).

3.2.2.
Correlations with temperature. When regional summer temperature and NDVI data were analyzed together, mixed model results showed that variations in July temperature were not significantly related with maxNDVI ( p > 0.05). However, when data were analyzed separately per region, large differences in correspondence of maxNDVI and minWSA with monthly air temperature data were observed, with some regions showing a high correspondence with temperature and other regions showing no correspondence (table 2). For NUN, July temperature was most important to annual maxNDVI (r = 0.79, p < 0.01, n = 11 years) whereas for YAK mean May temperature was most important to annual maxNDVI (r = 0.72, p < 0.05, n = 11 years). No relationship between maxNDVI and mean monthly temperature was observed for YAM and ALA. Summer warmth index did not correlate significantly with maxNDVI for any region. For interannual variations in growing season minWSA there were no significant relations with summer temperature variables, except for YAM showing a significant positive relationship with August temperature (r = 0.72, p < 0.05, n = 11 years).
We did not detect any significant temporal trends in regional average maxNDVI during the period 2000-10, neither did we detect trends in regional average minWSA for any of the four Arctic regions (data not shown). In contrast, we did observe an increase in SWI, with positive trends observed for ALA (0.91 SWI • C yr −1 ) and YAK (0.67 SWI • C yr −1 ) (figure 5).

NDVI
Previously, NDVI was found to be a good measure for shrub photosynthetic biomass in Alaskan tundra (Kushida et al 2009). Likewise, we observed a clear positive correlation between B. nana leaf fractional cover and NDVI in our experimental plots. However, our results also indicated that NDVI seemed to saturate at B. nana fractional cover values above 40%, corresponding to LAI values of 0.8-1.0 m 2 m −2 . This observation is in agreement with the saturation level of NDVI with LAI found by Van Wijk and Williams (2005). Both studies therefore seem to support that NDVI may not be useful for detecting future shrub expansion in tundra areas that are already high in green shrub biomass and LAI.  Summer NDVI has been shown to spatially correspond with SWI for pan-Arctic vegetation based on satellite data of 1993 and 1995 (Raynolds et al 2008). Jia et al (2006) found a high positive temporal-spatial correlation of NDVI with SWI in two transects of the Alaskan North Slope from 1991 to 2000. We indeed observed greater maxNDVI values at lower Arctic latitudes, suggesting that peak growing season biomass is positively related to summer temperature. Furthermore, NUN showed a strong correspondence between maxNDVI and July temperature. For YAK, maxNDVI was mostly affected by early growing season temperature, suggesting that an increase in growing season length may also influence peak biomass. In contrast, we did not observe any significant temporal correspondence between maxNDVI and SWI inter-annual variation for any of the four study areas for the period 2000-10. Despite positive trends in summer temperature for YAK and ALA during the last decade, no trends in maxNDVI were detected, suggesting that currently vegetation growth does not linearly respond to increasing air temperature. Of all regions studied, YAM showed least inter-annual variation in maxNDVI between vegetation classes and bioclimatic subzones. This may be a result of intensive reindeer grazing and/or other disturbance factors such as landscape erosion, resulting from increased permafrost thaw, affecting tundra vegetation and NDVI . Furthermore, YAM showed little inter-annual variation in summer temperature during the last decade, which may contribute to the low correspondence of maxNDVI with temperature in this region.
Several studies suggest that observed greening trends across the Arctic are related to an expansion of deciduous shrubs (Tape et al 2006, Sturm et al 2001, Jia et al 2009, Stow et al 2007, 2004, Jia et al 2003. Our analysis of satellite data reveals that shrub and graminoid tundra show highest maxNDVI values of all vegetation classes. This corresponds with a recent study by , which showed highest summer NDVI in shrub and tussock tundra for the pan-Arctic region. These results, combined with earlier observations of strong correlations between annual radial shrub growth and summer NDVI during recent decades (Forbes et al 2010, Blok et al 2011b, support the hypothesis that tundra greening is related to increasing deciduous shrub growth but does not exclude the possibility that this greening may also be related to an increase in graminoid biomass .

Albedo
An increase in deciduous shrub cover in the Arctic as a response to increasing air temperatures may lead to a decrease in summer surface albedo and cause atmospheric heating, thereby potentially triggering a positive feedback loop leading to a further increase in shrub growth . Taking NDVI as a proxy for green biomass and growth, we indeed observed a significant negative relationship between inter-annual variations in peak growing season maxNDVI and minWSA, which supports the hypothesis that a temperaturedriven increase in green biomass may lead to a reduction in summer albedo.
Unexpectedly, values of annual minWSA were on average higher in bioclimatic subzone E than in the cooler bioclimatic subzone D, suggesting that summer albedo does not linearly decrease with increasing summer temperature along a high to low-Arctic tundra gradient. Albedo may be influenced more by spatial extent and radiative properties of nonphotosynthetic landscape elements (e.g. water bodies, plant litter, bare ground) than biomass: toward the south, shrub fractional cover generally increases but fractional cover of water bodies decreases, which may explain the higher albedo values observed in southern tundra with relative high shrub fractional cover. Even though we masked out water bodies from the albedo maps using the highest-resolution land water mask available, sub-pixel sized water bodies may for a large part contribute to the overall albedo signal in wetland tundra and cause lower albedo values than observed in shrub tundra. In contrast to this observation, a study that examined differences in albedo between vegetation and climate zones across the Arctic during the snowmelt season, showed that subzone D had higher albedo than subzone E, presumably because of shrubs protruding through the snow layer in the more southern climate zone E .
We observed that peak growing season minimum albedo values are generally not lower in shrub tundra than in other tundra vegetation zones.
For three out of four Arctic areas investigated, lowest albedo values occurred in wetland tundra. Even though we masked out water bodies at a spatial resolution of 250 m, sub-pixel scale in undated areas and soil moisture may have been larger in wetland tundra than in shrub tundra, resulting in a greater absorption of solar radiation by the relatively dark water and soil surface and causing lower albedo values. From the vegetation classification map it is evident that in more northern areas wetlands tend to dominate. For the low-Arctic shrub tundra zones, our field data however suggest that a decrease in summer albedo can be expected with increasing shrub growth. Given the high temperature sensitivity of vegetation in the low-Arctic zone , these areas may experience the most rapid vegetation transition toward shrub dominance with increasing temperature, triggering a positive warming feedback through decreasing albedo. These results suggest that a northwards shift of shrub tundra might not lead to a decrease in summer minimum albedo during the snow-free season when replacing wetland tundra.

Conclusions and outlook.
From field data, we show that in mixed graminoid/shrub tundra vegetation, NDVI rises and albedo declines with increasing deciduous shrub cover. In support of this observation, MODIS satellite data (period 2000-10) from four large regions spread across the North American and Eurasian Arctic unequivocally show highest NDVI in areas classified as shrub and graminoid tundra. For albedo, wetland and barren tundra showed lower values than graminoid and shrub tundra. This suggests that an expansion of shrubs in tundra areas with climate warming might not in all cases lead to a lowering of summer albedo. However, not just changes in vegetation and regional climate may affect local summer NDVI and albedo, but also differences in soil pH, substrate chemistry, elevation and water availability (Lafleur and Humphreys 2008). Changes in vegetation composition, such as shrub encroachment, are influenced by climate. They may have a significant feedback to the climate system, but are at the same time an indicator for changing growing conditions. To evaluate the full climate feedback of shrub encroachment in the Arctic through changes in albedo and energy partitioning, integrative studies are needed. These necessarily include regional predictions of temperature, precipitation and cloud cover, changes in vegetation composition and biomass, as well as hydrological and permafrost conditions and related landscape changes.