The browning of Alaska’s boreal forest

We used twelve Landsat scenes from the 1980s–2009 and regional 2000–2009 MODIS data to examine the long-term trend in the normalized difference vegetation index (NDVI) within unburned areas of the Alaskan boreal forest. Our analysis shows that there has been a declining trend in NDVI in this region, with the strongest “browning trend” occurring in eastern Alaska where the climate during the growing season is relatively dry and warm. Possible reasons for the "browning trend" are decreased vegetation due to temperature-induced drought stress and increased infestations of insect pests.


Introduction
The greatest climate warming over the past 50 years in North America has occurred in Alaska and northwest Canada [1].This warming has led to substantial physical and biological changes including record sea-ice retreat [2], permafrost thawing [3], record summer warmth [4], and an increase in growing season length [5].In the arctic tundra, there has been an increase in vegetation productivity as indexed by the Normalized Difference Vegetation Index (NDVI) [6], consistent with field-based historic photography [7], and experimental warming studies [8].In contrast to the increasing NDVI in arctic tundra, there has been a declining trend in NDVI in boreal Alaska documented at several spatial scales [9][10][11].These boreal results were all based on the NDVI data produced as part of the National Aeronautics and Space Administration Global Inventory, Monitoring and Modeling Studies (GIMMS) project.However, this observed -browning‖ trend may be due to a bias towards negative trends in the OPEN ACCESS GIMMS NDVI data set.For example, a recent study [12] detected positive NDVI trends in re-vegetating burned areas in central boreal Canada using a 1-km AVHRR NDVI dataset, yet the GIMMS NDVI did not increase in most of these burned areas.The objective of this study is to determine whether the browning trends in Alaska's boreal region exists based on other NDVI sources (Terra MODIS and Landsat TM/ETM+ NDVI).
In this study, we analyzed the long term (mid-1980s to present) trend using NDVI computed from historic Landsat sensor data by developing a time series of NDVI for 12 different scene footprints located within boreal Alaska.Within each scene footprint, we also use MODIS data to document the NDVI trend since 2000.MODIS NDVI data is obtainable as composite data available every 16 days, while the Landsat sensor data is typically much more limited in Alaska.In this study, we acquired at least six years of Landsat images, with at least one image from the 1980s as a time series within each scene footprint.We also assess the potential of a browning bias in the GIMMS NDVI by comparing GIMMS and MODIS NDVI from a large regional perspective.We used polygons of burned areas since 1980 to exclude wildfire burns in our analysis.

Study Region
Most of Alaska's boreal forest lies in an intermountane region bounded to the north and south by large mountain ranges resulting in a west to east climatic gradient (Figure 1).The region is characterized by isolated mountains, large areas of hilly uplands, meandering rivers with broad floodplains and extensive wetland regions.Black spruce forests dominate cold sites, such as wetlands and north-facing slopes, often underlain by permafrost.White spruce occurs on warmer sites such as active floodplains and south-facing slopes.Deciduous forest (aspen, birch, and balsam poplar) also occur on warmer sites.The vegetation mosaic of this region is primarily controlled by disturbance legacies (primarily wildfire) and topographic control of microclimate.

GIMMS NDVI Data
The GIMMS 1981-2008 NDVI dataset was obtained from the GIMMS team [14], and also available from http://www.landcover.org/data/gimms/.The data are maximum NDVI values composited over a 15-day period, using Advanced Very High Resolution Radiometer (AVHRR) Global Area Coverage 1B data (4km resolution at nadir) and projected to a North American Albers Equal Area Projection with 64 km 2 pixel size.The data had been corrected for sensor degradation, inter-sensor differences, solar zenith angle, and viewing angle effects due to satellite drift.No atmospheric correction has been applied with GIMMS data, except for major volcanic stratospheric aerosol periods.Cloud screening was based on an AVHRR channel 5 thermal mask of 0 °C [14].

MODIS NDVI Data
The MODIS 1-km 2 NDVI product, MOD13A2, version5 [16] was used in this study, downloaded from https://wist.echo.nasa.gov/~wist/api/imswelcome/.This product is based on the Terra MODIS level 2 daily surface reflectance product, which provides red and near-infrared surface reflectance corrected for the effect of atmospheric gases, thin cirrus clouds and aerosols.This 16-day composite product was downloaded in a global sinusoidal projection and then reprojected to the Alaska Albers Equal Area projection at 1-km 2 pixel size.Only pixels with reliability value of zero were used in this study.To minimize the effect of unvegetated pixels, only pixels with an MODIS NDVI of at least 0.4 were used in this study.Forested pixels typically had maximum NDVI values above 0.8 for broadleaf forest and above 0.7 for spruce woodland.We computed the linear trend in annual maximum NDVI from 2000 through 2009 for each 1-km pixel in the MODIS NDVI time series.Using linear regression, we excluded all pixel trends that were not statistically significant (p > 0.05), and excluded any pixels that had burned since 1980.With the remaining pixels, we produced a raster of significant regression slopes for 1-km pixels within the Alaskan boreal region.We also computed linear regressions for each Landsat scene footprint, based on 2000 through 2009 MODIS maximum NDVI values.
The day of year for each NDVI pixel in the composite was retained to determine the time period during the growing season when boreal NDVI was near its maximum.This time period was used to define the annual peak growing season period that was then used to restrict our analysis to this period (Figure 2).

Landsat TM/ETM Data
The 1985-2009 Landsat NDVI data were compiled from Landsat 4 and 5 Thematic Mapper (TM) and Landsat-7 Enhanced Thematic Mapper Plus (ETM+) data downloaded from the United States Geological Survey (USGS) website from (http://glovis.usgs.gov/).The twelve scene footprints were chosen to span across the east to west climatic gradient (Figure 3) within the boreal zone of Alaska.Based on MODIS NDVI data, June 19 through August 8 represented the period of peak NDVI (Figure 2); thus Landsat TM/ETM+ scenes were chosen only if they were during this period and had less than 30% cloud cover.Landsat has a 16-day repeat cycle, and due to cloud cover, a sparse time series of TM/ETM+ was available for most areas.All TM/ETM+ data were spatially referenced to WGS 1984 UTM.A total of 106 Landsat TM/ETM+ images were processed (Table 1).For each Landsat WRS-2 scene, band 3 and band 4 digital numbers were converted to red and near-infrared at-sensor spectral reflectance [17] using sun angles and gain and bias values downloaded with the TM/ETM+ scenes.We then used a radiometric normalization method [18,19] based on dark and bright unvegetated surfaces to account for possible variation of atmospheric conditions between scenes.A master image, usually the one with the least cloud cover, was chosen from each time series of images (Table 1) and all other images within a given scene location were adjusted to the master images using radiometric rectification based on at least 30 dark (deep water bodies) and 30 bright (rock outcroppings or urban areas) unvegetated pixels.Linear regression was then used to adjust the red and near infrared spectral reflectances to the master image.Each regression had an R 2 of at least 0.95 with a sample size of at least 60.
After radiometric normalization, the adjusted red and near-infrared reflectance values were used to compute NDVI for each image.We eliminated cool cloud-contaminated pixels based on the radiant temperature threshold computed from thermal band 6 and we eliminated other unvegetated pixels based on a NDVI threshold specific to each scene to mask cloud shadows, water bodies and rock outcroppings.We excluded any areas that were burned since 1980 (http://agdc.usgs.gov/data/blm/fire/index.html ), since recently burned areas would have a sharp decline in NDVI, while older burned areas would have an increase in NDVI due to revegetation, and our primary focus was on the regional NDVI trend in unburned forests.With these areas eliminated, the mean NDVI was calculated for the entire image in the time series.Using simple linear regression, mean NDVI trend between 1985 and 2009 was then determined for each scene location, based on the slope and statistical significance of the regression line.For each scene we also computed the mean NDVI at each time period for those pixels that had reliable NDVI values (cloud-free, vegetated) throughout the entire time series.

Comparing GIMMS and MODIS NDVI Trends
Criticism of the GIMMS NDVI as a data set biased to -browning trends‖ in the boreal region [12], was based on a comparison of the coarse resolution GIMMS NDVI (64 km 2 pixels) with a finer resolution NDVI dataset (1 km 2 pixels) within relatively small burned areas.Thus scale may have been a factor in this reported -browning bias‖.We felt that a fairer comparison would be to examine mean GIMMS NDVI versus MODIS NDVI from a larger regional perspective, with regions larger than 15 million hectares.Therefore we compared the mean annual maximum NDVI within the eastern and western boreal region (Figure 1) based on GIMMS and MODIS NDVI from 2000 through 2008.Only pixels with annual maximum NDVI values above 0.4 were used in this analysis to minimize the effect of variations in unvegetated pixels.

GIMMS and MODIS NDVI Trends
There was a strong correlation (Pearson's r = 0.91, p < 0.001) between GIMMS and MODIS annual maximum NDVI values within the eastern boreal region.However, there was a weak correlation (Pearson's r = −0.12,p = 0.750) from the western boreal region (Figure 4).One possible explanation for this pattern might be more cloud-contamination of the GIMMS data in the western boreal region.For example, in 2001 and 2007 the GIMMS NDVI decreased, while the MODIS NDVI increased from the western boreal region.During the late July 2001 and 2007 composite periods, most of the pixels within the western boreal region were flagged as unreliable in the MODIS product, likely due to cloud contamination (Figure 5), while the GIMMS product flagged most of these pixels as good quality, resulting in reduced mean NDVI within the region.Because of the potential of cloud-contamination artificially causing a -browning trend‖ based on GIMMS NDVI data, we restricted the remaining analysis to MODIS and Landsat TM/ETM+ NDVI values.

Longer Term Trends Based on Landsat TM/ETM+ NDVI
The scenes with the strongest -browning trend‖ or decline in NDVI were from the eastern boreal zone where summer growing season is the hottest and driest within boreal Alaska (Table 2; Figure 3), with no significant trends in western boreal Alaska based on MODIS and TM/ETM+ NDVI.The number of pixels of each Landsat TM/ETM+ image used to calculate NDVI could be as low as one third of the scene (Table 1).The varying amounts of pixels used was due to high cloud contaminations, smoke from active forest fires, and omitting previous fire areas.This shifting of usable pixels was likely a source of substantial variation in the TM/ETM+ NDVI time series.To address this issue we evaluated the pixels that were vegetated and used in every scene in the time series.This left 0.43 to 9.5% of the scene available to be analyzed.As seen in Table 2, there was little change between using all available pixels and using only reliable pixels from all time periods even though many fewer pixels were used.The scenes in eastern boreal Alaska that had a strong significant negative trend when looking at all pixels kept similar trends when assessing only the common pixels.The regional pattern of pixel-level trends based on linear regressions of MODIS 2000-2009 NDVI also reflects the west to east climatic gradient, with pixels having the strongest declining trend in NDVI occurring in the warmest, driest region of boreal Alaska (Figure 6).

Potential Causes of NDVI Browning Trends
Although wildfire occurs throughout boreal Alaska, we excluded all areas within burn perimeters since 1980, and therefore the NDVI trends are not due to wildfires.There are several likely reasons for the decreasing NDVI trends in boreal Alaska.Based on tree ring and stable isotope data [4], summer warming in eastern boreal Alaska since the mid-1970s is unprecedented over the past 100 years and 2004 was the warmest summer in the past 200 years.With warmer summers and longer growing seasons, insect life cycles can be shortened and survival increased, leading to major new insect outbreaks especially in eastern Alaska.For example, since 1989 there has been a dramatic increase in the area of boreal forest infested especially with the aspen leaf miner, the willow leaf blotch miner and the spruce budworm (Figure 7).Examples of growing pest infestations in boreal Alaska.Data acquired from http://agdc.usgs.gov/.Aspen leaf miner (A) infests primarily quaking aspen and balsam poplar while willow leaf blotch miner (B) infests all species of willow excluding feltleaf willow.Spruce budworm (C) infests primarily white spruce, a species also experiencing widespread regional growth decline based on tree-ring studies.
The appearance of willow leaf blotch miner and spruce budworm infestations in boreal Alaska is new and did not occur prior to 1990.The prolonged multi-year infestation of aspen leaf miner is also new, with previous outbreaks crashing after one or two years.
The aspen leaf miner (Phyllocnistis populiella) hardly existed in Alaska in 2000 but in 2007 it occupied over 300,000 hectares.The aspen leaf miner can be found feeding on the epidermal cells on both the top and bottom of aspen leaves and gives the leaves a silver appearance.They do not enter the mesophyll cells and thus do not significantly impact NIR reflectance.However, damage done to the leaf epidermis can damage stomates and affect leaf water.The trees are able to survive when the aspen leaf miner is present, but stunted growth has been found and highly damaged leaves abscise about four weeks earlier than unharmed leaves [20].
Willow leaf blotch miner (Micrurapteryx salicifoliella) has also expanded significantly since it was first found in Alaska in 1991 [21].Willow are broadleaf shrubs that would likely contribute substantially to the NDVI response in black spruce stands and other wetland types.Eggs are cemented to the underside of the leaf in late May and hatch into the leaf in early June.They feed on the mesophyll which would likely cause a decline in NIR reflectance and NDVI by mid-summer.
Spruce budworm (Choristoneura fumiferana and Choristoneura orae) outbreaks have been very distructive across boreal Alaska.Eggs are deposited a few inches from where the new bud will grow and hatch over a period of a few weeks so that some will hatch when the buds emerge [22].All ages of trees are suseptible to spruce budworm but the more mature spruce stands have a better chance of survival [23].Population fluxuations have been observed from 5 to 40 years [22] which is more reflected in the Figure 7 where the spruce budworm has been present in the late 1980's while the aspen leaf miner and willow leaf blotch miner were not.
Summer water deficits in eastern boreal Alaska have increased in recent years due to increased evapotranspiration, and decreased regional water balance associated with longer, hotter summers [4,24].This has likely lead to decreased tree growth due to drought stress.Based on tree-ring studies, there has been a widespread negative correlation between tree growth and summer temperature in both white spruce [25] and black spruce [26].Growth of boreal deciduous species is thought to be less sensitive to the higher summer temperatures, but more sensitive to drought since the relatively high transpiration rates in deciduous stands hastens the onset of soil water deficits in early summer [27].

Conclusions
Based on both Landsat TM/ETM+ and MODIS NDVI, there has been a browning trend in eastern boreal Alaska.The strongest browning trend occurred in the warmest, driest region of boreal Alaska and at the lowest elevation zones within that region.There are many potential reasons for the observed browning trend including decreased vegetation production due to temperature-induced drought stress, reduction of phytomass due to increased insect infestations, and increased tree and shrub mortality due both to drought stress and insect infestations.There were no significant NDVI trends in western Alaska.This was likely due to the more maritime climate of this boreal region leading to less impact of a warming climate and associated drought-stress and increased insect infestations relative to eastern boreal Alaska.
Although Landsat and MODIS NDVI gave similar results, each had unique advantages.Landsat TM/ETM+ with a higher spatial resolution allowed for the elimination of cloud and cloud shadows, as well as small unvegetated areas.With TM data available since the 1980s, longer-term trends are possible to document using Landsat TM/ETM+ data.Landsat scenes typically had large amounts of cloud and cloud shadow contamination while the MODIS NDVI product was composited over a 15 day-period and available for several composite periods during each growing season since 2000.The GIMMS NDVI data were not significantly correlated with MODIS NDVI data for the eastern boreal region, likely due to cloud contamination substantially influencing GIMMS NDVI values.

Figure 2 .
Figure 2. Mean MODIS NDVI by day of year of NDVI observation for two Landsat scenes.Based on these data, we defined our annual peak growing season study period as from day 170 (June 19) to day 220 (August 8).Day to day variation in mean NDVI is likely due to daily variation in proportion of conifer versus broadleaf pixels.Interannual variation is likely due to climatic conditions, with 2007 being an earlier spring start of growing season.(A) P75R15 (WRS Path75 Row15) covers a western part of the region.(B) P70R14 (WRS Path70Row14) covers an eastern part of the Alaskan boreal region.

Figure 3 .
Figure 3. Location of twelve Landsat scenes used in this study.

Figure 4 .
Figure 4. Mean annual maximum NDVI from MODIS and GIMMS NDVI products within 2 large boreal regions (eastern boreal region = 34.5 million ha, western boreal region = 19 million ha).The lower NDVI and poor interannual correlation of the GIMMS Western Boreal NDVI is likely due to clouds reducing the GIMMS NDVI in this region.

Figure 5 .
Figure 5. MODIS and GIMMS mean NDVI by composite periods of 2001.The period of maximum NDVI is typically from late June to early August in boreal Alaska.Declines in mean NDVI during this period is likely due to clouds reducing the mean NDVI value.

Figure 6 .
Figure 6.Pixel-level regression slopes based on 2000-2009 MODIS annual maximum NDVI values.Only regression slopes with p-values of <0.05 are displayed in this figure.Areas burned since 1980 were excluded and are portrayed as black polygons in this figure.

Figure 7 .
Figure 7.Examples of growing pest infestations in boreal Alaska.Data acquired from http://agdc.usgs.gov/.Aspen leaf miner (A) infests primarily quaking aspen and balsam poplar while willow leaf blotch miner (B) infests all species of willow excluding feltleaf willow.Spruce budworm (C) infests primarily white spruce, a species also experiencing widespread regional growth decline based on tree-ring studies.

Table 1 .
Landsat scenes used in this study.Asterisk * indicates the image used as master for radiometric normalization.

Table 2 .
Interannual linear trends of NDVI in boreal Alaska.(A) Linear regressions for each of the twelve Landsat scene footprints, based on TM/ETM+ images from the 1980s to 2009.(B) Linear regressions based on pixels that had reliable NDVI values in every time period.(C) Linear regressions of the twelve Landsat scene footprints, based on annual maximum MODIS NDVI from 2000 to 2009.