Detecting long-term changes to vegetation in northern Canada using the Landsat satellite image archive

Analysis of coarse resolution (∼1 km) satellite imagery has provided evidence of vegetation changes in arctic regions since the mid-1980s that may be attributable to climate warming. Here we investigate finer-scale changes to northern vegetation over the same period using stacks of 30 m resolution Landsat TM and ETM + satellite images. Linear trends in the normalized difference vegetation index (NDVI) and tasseled cap indices are derived for four widely spaced national parks in northern Canada. The trends are related to predicted changes in fractional shrub and other vegetation covers using regression tree classifiers trained with plot measurements and high resolution imagery. We find a consistent pattern of greening (6.1–25.5% of areas increasing) and predicted increases in vascular vegetation in all four parks that is associated with positive temperature trends. Coarse resolution (3 km) NDVI trends were not detected in two of the parks that had less intense greening. A range of independent studies and observations corroborate many of the major changes observed.


Introduction
The Arctic has experienced marked temperature increases in recent decades that have been larger than the global average (McBean et al 2005). This has been linked to widespread changes in permafrost condition, sea ice and snow cover, glacier mass balance, and tundra vegetation growth (Walsh et al 2005, Richter-Menge andOverland 2010). These impacts are predicted to intensify, as climate scenario modeling projects further arctic temperature increases of 5-7 • C by the end of this century (Kattsov and Kallen 2005).
Satellite remote sensing represents a promising, efficient means of documenting changes occurring over vast and remote arctic regions. Indeed, much of the evidence for recent changes in the region comes from this technology (Stow et al 2004). One such change is the 'greening' of tundra vegetation documented by studies employing the normalized difference vegetation index (NDVI) calculated from the ∼25 year archive of coarse resolution NOAA advanced very high resolution radiometer (AVHRR) satellite imagery (Myneni et al 1997, Goetz et al 2005, Jia et al 2009. The NDVI, providing a surrogate for the amount of green foliage, has demonstrated broad circumpolar increases that have been attributed to longer growing seasons and/or higher temperatures. Increased growth of graminoids and arctic shrubs, such as willow ( Salix spp.) and birch (Betula spp.), has also been recently documented at local scales based on repeat photography (Tape et al 2006) and in situ measurements (Gould et al 2009, Hudson and Henry 2009, Hill and Henry 2011.
Although the AVHRR sensor has the advantage of providing daily, continental-scale coverage, its large heterogeneous pixels permit only identifying regional-scale changes. The Landsat TM and ETM+ satellite sensors provide an image archive going back to 1984. These sensors cover a 185 km wide swath, have a 16 day revisit, and provide six 30 m resolution reflectance channels. This freely accessible archive provides an opportunity to study decadal vegetation changes at finer spatial scales compared to AVHRR. However, one challenge in employing Landsat for this purpose in the Arctic is its limited observation frequency in an environment where vegetation has high inter-annual variability (Epstein et al 2004) and a short growing season. Dense time-series 'stacks' of Landsat images have recently been used in temperate environments for detecting both long-term vegetation trends (Röder et al 2008) and land cover transitions (Huang et al 2010, Kennedy et al 2010. The main advantage of the stack approach is that the greater number of observations can allow long-term changes to be detected with greater sensitivity and reliability by comparison to conventional two-date change detection. This letter presents a Landsat-based, vegetation change study over four, widely spaced national parks in northern Canada. Linear trends in several Landsat spectral indices are analyzed from stacks of 16-41 growing season images spanning periods of 17-25 years. Trends are summarized by land cover type and also related to changes in the fractional cover of functional vegetation types using regression tree models (Olthof and Fraser 2007). The results are assessed by comparing them to available independent studies and observations.

Study areas
Study areas were focused on four Canadian national parks: Ivvavik, Sirmilik, Torngat Mountains and Wapusk (figure 1). This work was performed as part of a collaboration with Parks Canada Agency and the Canadian Space Agency that since 2004 has been developing remote-sensing-based methods to track and report on changes in the ecological integrity of Canada's national parks (Fraser et al 2009, Woodley 2010. The four study parks, covering five terrestrial ecozones (table 1) (Wiken 1986), were selected to represent a wide range of northern ecological and physical conditions. All 17 classes from the Circumpolar Arctic Vegetation Map (Walker et al 2005) are represented by the study areas, with the exception of prostrate/hemiprostrate dwarf shrub tundra. The limit of tree growth occurs within three of the parks (Ivvavik, Wapusk and Torngat) and all four contain shrub tundra vegetation communities. Wapusk encompasses a plain within the Hudson James lowlands and has extensive wetland ecosystems. The other three parks include both mountainous and flat, coastal regions. A detailed description of the vegetation and physical landscape within each park is provided by Parks Canada Agency (2011).
All 78 images were Level 1T, precision-geocoded and terrain-corrected products. These were supplemented with 23 images from the Canada Centre for Remote Sensing Landsat archive. We exploited the overlap between adjacent orbital paths to maximize observation frequency over the study areas (figure 1). Between 16 and 41 images from the nominal July-August growing season were analyzed over each park for collection periods spanning 17-25 yr (table 1).
Images that deviated strongly from peak annual vegetation phenology were removed to ensure that the selection of images did not introduce artificial trends. Deviation was calculated for green targets in each park by subtracting average, acquisitiondate NDVI from peak NDVI derived from 10 day AVHRR  (Fraser et al 2012). Note that none of the four stacks of Landsat images analyzed showed a significant ( p < 0.05) linear trend for year of Landsat acquisition versus deviation from peak NDVI.
Satellite digital numbers were converted to top-ofatmosphere (TOA) reflectance using USGS calibration coefficients designed to provide compatibility between TM and ETM+ imagery (Chander et al 2004). Clouds and shadow were manually masked by digitizing on-screen against a clear-sky reference Landsat image. Missing data in ETM+ images following the 2003 scan line corrector failure were also masked. The NDVI and brightness, greenness, and wetness tasseled cap (TC) indices were computed from the TOA channels. TC indices are linear transformations of the six optical Landsat channels that can be interpreted in terms of physical image characteristics (Crist and Cicone 1984).

Change detection approach
Change detection involved analyzing each pixel's unique timeseries of values in the stack of masked Landsat images (Fraser et al 2012). Linear trends in the four indices were quantified using Robust (TheilSen) regression analysis. Significant ( p < 0.05) trends were then related to predicted changes in fractional vegetation and land covers using regression tree modeling.
Regression tree classifiers were created for Sirmilik and Torngat using a finer resolution (2.4-10 m) baseline image and field plot measurements providing per cent composition. The fine resolution, fractional products were then aggregated within coincident 30 m Landsat pixels to train a second set of tree classifiers using the trend-based Landsat values. Since suitable plot measurements were not available for Ivvavik and Wapusk, other available reference data were used to create 'hard' classifications of fine resolution imagery (table 1). These were then aggregated to 30 m to train the Landsat classifiers.
The regression tree models were applied using trendbased values corresponding to the first and last Landsat image dates. If no significant trend was detected, the average Landsat index value was inputted. A simple differencing of predicted first-and last-date fractions was used to estimate fractional vegetation cover change.

Spectral trends
Long-term spectral index trends derived from 30 m Landsat and 3 km AVHRR imagery  are shown in table 2. The table shows the percentage of the study areas over which each index is did not change, decreased (negative slope), or increased (positive slope) at a 95% confidence level. It also shows the average slope for each index, where non-significant slopes are assigned a value of zero.
AVHRR-NDVI change statistics were similar for Sirmilik and directionally consistent for Torngat Mountains by comparison to Landsat-NDVI (table 2). We also found that the distribution of AVHRR-NDVI trends were similar to greenness trends from Landsat for these two study areas (figures 2(a)-(d)), except that the greening valleys in Torngat were poorly discerned at 3 km. The Landsat-NDVI trends in Ivvavik and Wapusk were not detected at 3 km, likely due to the weaker greenness trends observed in these study areas (table 3).
Landsat TC greenness (TCG), representing the contrast between the near-infrared and other five Landsat channels, consistently showed more areas increasing (6.1-25.5%) than decreasing (0.3-4.1%). The NDVI, which measures the difference between the near-infrared and red bands that is characteristic of green foliage, showed similar results to TCG for all study areas, both statistically and spatially. In three of the parks (Sirmilik, Torngat and Ivvavik), areas of increasing greenness corresponded to relatively moist, lower elevations (figures 2(b), (d) and (f)). Greening in Wapusk was more highly concentrated in the northern portion of the park lying above the treeline, with some patches in the forested zone associated with fire regeneration (figure 2(e)). The only major Table 2. Summary of Landsat-based changes in spectral indices and estimated vegetation and land cover fractions. For each index or fraction, the table shows the percentage of the analysis area that is non-changing (%NC), decreasing (%Dec), and increasing (%Inc). The row labeled 'Avg' represents average change for each index or land cover fraction over the study area, where no-change pixels are assigned a value of zero.

Change in spectral indices
Change in cover fractions  (Brook and Kenkel 2002) was cross-walked to the same legend for consistency.
The largest increases in greenness occurred in the higher biomass graminoid and shrub land covers (1-6, 11) (table 3). This suggests that increased growth of photosynthetically active vegetation is occurring most strongly in areas with already favorable growing conditions. More sparsely vegetated classes (7-10, 12) generally showed smaller changes in TCG. Overall, Torngat and Sirmilik had more strongly positive slopes by comparison to Ivvavik and Sirmilik.
The brightness (TCB) index overall showed the largest areas of significant change. Three of the study parks (Ivvavik, Sirmilik and Wapusk) had larger percentage decreasing trends (13.9-17.5%), while Torngat had a larger percentage increasing trend (5.3%). TCB represents brightness in all six optical Landsat channels and may be influenced by many surface and vegetation properties. Areas of declining TCB in Ivvavik and Sirmilik were associated with predicted increases in vascular vegetation cover (section 3.2). This would be consistent with increased growth of shrub vegetation and resulting canopy shadowing, especially if it were occurring over a lighter soil or rock background.
TC Wetness (TCW) contrasts the two shortwave-infrared channels with the other channels and is sensitive to moisture content of soil and vegetation and to vegetation density and structure. TCW trends were more difficult to interpret and did not show any consistent patterns among the study parks: they decreased mainly on south-facing mountain slopes in Ivvavik, were highly variable in Wapusk's wet plain, and increased in many of Sirmilik's and Torngat's areas of greening.

Fractional land cover trends
Regression tree models of fractional pixel land cover were created by scaling-up plot and/or high resolution estimates of land cover to 30 m and relating this to the trend-based Landsat indices. Applying the models to the trend values for the first and last Landsat dates provided an estimate of how the spectral trends translated to changes in land and vegetation cover composition.
The proportion of each study area over which fractions did not change, decreased, and increased, as well as the average, absolute change for three key cover types within each park is shown in table 2. Models were assessed by comparing predicted fractions against 'real' up-scaled fractions from a 20% hold-out testing sample. The correlation coefficient averaged among all models created for each park ranged from 0.73 to 0.92 and average absolute error ranged from 2.3 to 4.4% (table 2). Regression models derived for Wapusk using available field data were poor (r 2 = 0.29-0.44) and therefore not applied to the Landsat trends.
All three modeled study areas had an overall increase in estimated shrub cover fraction (vascular for Sirmilik) and a decrease in bare land cover fraction during the study periods. Shrub increases were generally predicted to occur in areas of increasing TCG and/or decreasing TCB. Figure 3 shows fractional changes for the Torngat study area, which were the most striking among the parks. Here, the major trends were replacement of bare cover by shrub in the broad valleys and also by herbaceous vegetation on the lower hillslopes. The more sparsely vegetated, higher elevations were estimated to have relatively stable fractional composition. In these areas, environmental factors other than temperature, such as soil moisture, wind exposure and nutrient flux, may be more limiting to vegetation growth (Walker 2000).
Similar to the greenness indices, we found that vascular vegetation fractions increased mainly in areas already supporting productive vegetation communities. Areas with unfavorable soil and moisture conditions, especially upland sites, were generally stable. What then, could be causing this apparent increase in greenness and vascular vegetation growth across the study areas during the last 17-25 years? The simplest explanation is that this is a consequence of warmer temperatures.
All four study areas showed a trend toward milder winters during the study periods. Figure 4 plots an annual temperature index based on cumulative temperatures below 0 • C for each year derived from NCEP North American Regional Reanalysis (NARR) daily mean data. Note that the slopes are larger for Torngat and Sirmilik, which also had the stronger greening trends. An annual summer warmth index based on cumulative daily temperature above 0 • C was similarly positive for Torngat ( p = 0.003), Sirmilik ( p = 0.09) and Wapusk ( p = 0.06) (summer NARR data is missing over Ivvavik). Epstein et al (2008) documented a strong spatial relationship (r 2 = 0.84) between total aboveground biomass and summer temperatures across a 1800 km north-south transect through North American arctic tundra. Several mechanisms have been proposed for how increasing temperatures may stimulate growth of arctic vascular vegetation. As a biome limited by temperature, a longer and/or warmer growing season would be expected to directly increase the growth rate of woody shrubs (Walker 1987, Forbes et al 2010. Warming may also have an indirect effect on vegetation by increasing ground temperatures and nitrogen nutrient mineralization by soil microbes (Chapin et al 1995). This process may be positively reinforced by shrub growth as this leads to additional trapping of insulating snow (Sturm et al 2005). Herbivores may have also influenced trends in vegetation abundance in some of the study areas. The southern coastal plain of Sirmilik's Bylot Island hosts a large colony of Greater Snow Geese that consume a significant proportion of the graminoid vegetation biomass near their nesting sites (Gauthier et al 1995). The Bylot goose population peaked in 1993 at 150 000 and then declined to 106 000 in 2003 and 94 000 in 2008, likely resulting in a reduction in overall grazing pressure. Similarly, the George River caribou herd that migrates into Torngat during the summer has dropped from almost 800 000 animals in 1993 to less than 75 000 in 2010. Grazing by this herd within its summer range causes reductions in biomass of dwarf birch (Manseau et al 1996) and therefore the population decline may be promoting shrub regeneration.

Corroboration of results
A higher confidence should be assigned to the physical, spectral trend results by comparison to the modeled, fractional changes. The Landsat image processing methods, including calibration and phenology screening, were unlikely to lead to systematic biases in the index trends. Also note that the trend patterns were spatially coherent and not random or scene-wide, as would be expected from sources of noise or bias.
A major challenge for validating the satellite-based evidence for arctic vegetation growth is the paucity of repeat ground-based measurements. The number of long-term, plot-based studies measuring vegetation change, although increasing (Gould et al 2009, Hill andHenry 2011), is still very small across the tundra biome. For each of the study areas, we attempted to locate published studies and reports related to vegetation and other surface changes that could corroborate (or contradict) the Landsat trend and fractional change results.
The Landsat greening trends observed over coastal areas of Ivvavik and nearby Herschel Island are consistent with field observations recorded by Kennedy (2008). A sampling of dominant vegetation types in 1985 and 1999 on Herschel and in 1988-9 and 2001 on the coastal plain suggested increases in the cover of lupine, polargrass and netted willow. These observations also matched the predicted increasing shrub and herbaceous fractions in the sampled areas. We also noted that long-term coastline erosion and retrogressive slumping documented on Herschel (Lantuit and Pollard 2008) could be clearly identified as areas with predicted decreasing shrub and herbaceous fractions and increasing bare and water fractions.
The only systematic vegetation monitoring within Sirmilik is being performed on Bylot Island's lowland plain. Annual plot measurements of wetland graminoid vegetation have been conducted since 1990 and indicate an increased aboveground live biomass of 84% between 1990 and 2008 (Cadieux et al 2008). Our analysis showed extensive greening trends ( figure 2(b)) and vegetation fraction increases over this same region. The reduction in glacier extent on Bylot described by Dowdeswell et al (2007) was associated with the only extensive negative greenness trends in the study areas ( figure 2(b)). These areas also had predicted decreasing ice cover and increasing water cover fractions.
In Wapusk, the large patches of positive and negative greenness trends in the southern, treed portion could be matched to fire disturbance records from the Canadian National Fire Database. An area of declining greenness in the northwest limit of the park was related to severe damage of coastal wetland vegetation caused by snow geese (Abraham et al 2005). The other greenness changes in the park, which mainly corresponded to increases, could not be assessed using available ground-based information.
Finally, for Torngat, traditional knowledge provides some support for the most dramatic greening observed among the study areas. Inuit Elders have observed that the landscape is becoming greener and that there are more and larger woody shrubs, such as willow in the park (Parks Canada Agency 2008). This is corroborated by a pair of replicate photographs taken by Parks Canada employees in 1991 and 2011 that shows increased density and height of birch and willow shrub (figure 5). A series of aerial photographs from 1964 and 2003 around George River 140 km to the west similarly shows that erect shrubs (mainly dwarf birch Betula glandulosa) have substantially increased in density over more than half of available surfaces (Tremblay et al 2011).

Conclusions
This study demonstrated that the Landsat TM/ETM+ satellite image archive can be effective for identifying trends since 1984 in the greenness of arctic tundra vegetation. The 30 m resolution trends provide the ability to more closely investigate regional-level productivity changes that can complement studies based on 1 km AVHRR imagery. For example, results from two of our study areas showed that Landsat imagery can identify weaker or more diffuse greenness changes that may be undetectable using 1 km imagery. The Landsat spectral trends can also be converted to an estimate of fractional vegetation and other land covers by up-scaling plot and/or fine resolution land cover maps to train a regression tree classifier.