Trends and transitions in the growing-season normalized difference vegetation index (NDVI) from the Moderate Resolution Imaging Spectroradiometer (MODIS) satellite sensor at 250-m resolution were analyzed for the period from 2000 to 2018 to understand recent patterns of vegetation change in ecosystems of the Arctic National Wildlife Refuge (ANWR) in Alaska. Statistical analysis of changes in the NDVI time series was conducted using the breaks for additive seasonal and trend method (BFAST). This structural change analysis indicated that NDVI breakpoints and negative 18-yr trends in vegetation greenness over the years since 2000 could be explained in large part by the impacts of severe wildfires. At least one NDVI breakpoint was detected in around 20% of the MODIS pixels within both the Porcupine River and Coleen River basins of the study area. The vast majority of vegetation cover in the ANWR Brooks Range and coastal plain ecoregions was detected with no (positive or negative) growing-season NDVI trends since the year 2000. Results suggested that most negative NDVI anomalies in the 18-yr MODIS record have been associated with early spring thawing and elevated levels of surface moisture in low-elevation drainages of the northern ANWR ecoregions.
The Arctic National Wildlife Refuge (ANWR) was established by the Alaska National Interest Lands Conservation Act of 1980 and covers 19 million acres (77 000 km2) in northeast Alaska. Proponents of development in the ANWR view its 1.6 million acre (6475 km2) coastal plain as a promising onshore oil reserve (Comay et al. 2018). Nonetheless, wildlife habitats in the ANWR are vulnerable to long-lasting effects from any disturbance, in part because short growing seasons in the Arctic provide limited time for species to recover.
Caribou (Rangifer tarandus) are the most numerous large animals in the ANWR. The Porcupine herd (named after the Porcupine River) numbers approximately 150 000 animals. The herd’s Arctic Village–South Brooks Range spring migration route within Alaska and the ANWR crosses the East Fork of the Chandalar River, the Sheenjek River, and upper Coleen River, and follows the Firth River into Canada where it joins the Old Crow migration route. Caribou herds will use areas in the northern foothills of the Brooks Range during summer at elevations of 300–600 m, with vegetation communities consisting of graminoid meadows, dwarf shrub, and alpine tundra communities (Nicholson et al. 2016).
Caribou consume a variety of vegetation including fungi, lichens, woody browse, graminoids, and forbs (Thompson et al. 2015; Denryter et al. 2017). Barboza et al. (2018) reported that female caribou and their calves in northern Alaska select a mixture of graminoids, browse, and forbs to achieve adequate dietary concentrations of digestible energy and nitrogen requirements. Graminoids were the most abundant forage, accounting for 77% of the digestible nitrogen and 74% of the digestible energy in forage biomass.
Localized patterns of vegetation growth and decline in the ANWR may be affecting the survival and reproduction of the caribou herds, which are a crucial subsistence resource for Alaska native communities in the region. Because of changing growing-season lengths, caribou migration may be increasingly out of sync with the timing of key food resources and forage availability. Parturition (birthing of calves) precedes or coincides with the onset of the Arctic growing season and is timed to maximize the period of peak nutrient availability to mother–offspring pairs (Loe et al. 2005).
The trophic mismatch hypothesis for large herbivores postulates that earlier green-up will shift peak nutrient availability away from peak nutritional demand leading to lower productivity. Recent work by Gustine et al. (2017) showed that climate-related effects on forage in the summer and autumn ANWR migration ranges corresponded more closely with the demands of female caribou and their offspring to gain mass for the next reproductive cycle. Therefore, these researchers implied that the window of time to examine the mismatch hypothesis in Arctic ungulates is not at parturition, but rather in late summer–autumn, at which time the effects of small changes in forage quality are amplified by forage abundance, peak forage intake, and resultant mass gains in mother–offspring pairs.
In a study of land-cover change using high-resolution aerial imagery for the ANWR, Jorgenson et al. (2018) estimated that 18% of the area had been altered over the past 50 years, mainly by wildfire and postfire succession, shrub and tree increase in the absence of fire, river erosion and deposition, and ice-wedge degradation. Changes in tundra communities tended to be related to landscape wetting, mainly from increased wet troughs caused by ice-wedge degradation. Boreal forest cover showed changes associated with landscape drying and decreases in lake area. Prior to this, Pattison and Welker (2014) found that decreasing snow cover was detrimental to tundra plants in northern Alaska, particularly for various species of shrubs and grasses, including the diamond-leaf willow (Salix pulchra) and tussock cottongrass (Eriophorum vaginatum). As snow depth declined, so too did general productivity of these species throughout the growing season.
To evaluate satellite imagery for vegetation greening trends in northern Alaska, Pattison et al. (2015) reported on changes in species composition for coastal plain field plots in the ANWR from 1984 to 2009. They linked these changes to the trends in the normalized greenness vegetation index (NDVI) at both fine and coarse scales. Field plot data implied few changes in plant community composition and no detectable increases in total vegetative cover. Deciduous shrub cover did not show large increases. Field-measured NDVI was positively correlated to deciduous and evergreen shrub composition, suggesting that these functional groups had a strong influence on NDVI values, but the field data showed that NDVI did not change over time on these tundra types. Although Landsat (30 m) NDVI was consistent with field-measured NDVI, NDVI at 1-km resolution from the National Oceanic and Atmospheric Administration Advanced Very High Resolution Radiometer (AVHRR) satellite record was not. AVHRR NDVI values showed increases that were in neither the field-measured nor Landsat NDVI records. This result suggested that AVHRR may be demonstrating increasing trends in NDVI that are not occurring on the ground in some Arctic tundra ecosystems.
To date, most of the large-scale studies of vegetation greening or browning in Alaska have not included comprehensive structural breaks analysis, designed to simultaneously detect all major disturbances that can alter greening trend statistics and the conclusions about gradual change in vegetation cover density and tundra or shrubland health. Gradual change analysis applied to a time series is designed to test for changes in the coefficients of a regression model, and generally assumes that there is just a single change under the alternative or that the timing and the type of change are known (Zeileis et al. 2002). A structural break can occur when a time series abruptly changes at a point in time. Detection of multiple breaks or disturbances in a time series of NDVI can occur in wilderness areas as a result of periodic wildfires, insect outbreaks, and/or from repeated cycles of extreme weather events.
The objective of this study was to detect both abrupt and gradual changes in vegetation cover throughout the ANWR since the year 2000 using the 250-m-resolution regional MODIS NDVI record for structural change analysis. The principal question posed in this analysis of the highest-spatial-resolution MODIS NDVI available, and the longest time series yet tested, was, Has the green vegetation cover changed over large areas in the ANWR since the year 2000 and particularly along major caribou migration routes? Statistical analysis of changes in the MODIS NDVI time series was conducted using the breaks for additive seasonal and trend method (BFAST; Verbesselt et al. 2010a,b; de Jong et al. 2012).
2. Study area
The ANWR covers nearly 80 000 km2 from the Beaufort Sea coast in the north, across the Brooks Range to the boreal forest and tributaries of the Yukon River in the south (Figure 1). The mean annual temperature is below freezing, and all parts of the ANWR are underlain by continuous permafrost, except for larger river valleys in the far southern portion (Jorgenson et al. 2018). Tundra vegetation predominates from the Brooks Range to the coastal plain, composed mainly of dwarf shrubs, sedges, and mosses (USFWS 2015). Boreal forest vegetation is found south of the Brooks Range. White and black spruce forests predominate in the northern portion of the subarctic lowlands, whereas the uplands support moist sedge–shrub tundra.
The ANWR comprises five ecoregions in a roughly north-to-south direction: the Beaufort Sea coastal plain, the Brooks Range Foothills, the Brooks Range, the Davidson Mountains, the Yukon–Old Crow basin, and the North Ogilvie Mountains (Nowacki et al. 2002). The coastal plain is composed of extensive networks of ice-wedge polygons, peat ridges, frost boils, and pingos (ice-cored hills) (USFWS 2015). Topography throughout the Brooks Range is rugged, with a history of glaciation and differential erosion of tilted, folded, and faulted rock layers. River valleys are wide and flat floored, cut by glaciers and then filled with alluvium (Figure 2).
Despite its remote and rugged topography, the ANWR has not been impervious to substantial changes over the past 50 years. For instance, Jorgenson et al. (2015) reported that the mean annual temperature at the Kuparuk weather station on Beaufort Sea coastal plain increased by 2.5°C between 1984 and 2009. Rapid increases in May air temperature have been associated with earlier snowmelt dates in northern tundra zones, which advanced by about 10 days between 1941 and 2004 (Hinzman et al. 2005). Gustine et al. (2017) also determined that day of spring ground thaw (≥0°C) occurred 8 days earlier on average and the length of the vegetation growing season was 11 days longer in 2013 than in the 1970s. Fast snowmelt transition has been tied to the presence of anticyclones in arctic tundra, which can in turn affect summer soil moisture across those latitudes with an excess in latent heat return to the atmosphere in midsummer (Fochesatto et al. 2015).
3.1. Fire boundary data from Landsat
Digital maps of burn severity classes at 30-m spatial resolution were obtained from the Monitoring Trends in Burn Severity (MTBS; www.mtbs.gov) project, which has consistently mapped fires greater than 1000 acres across the United States from 1984 to the present (Eidenshink et al. 2007). The MTBS project is conducted through a partnership between the U.S. Geological Survey (USGS) National Center for Earth Resources Observation and Science (EROS) and the USDA Forest Service. Landsat data have been analyzed through a standardized and consistent methodology by the MTBS project. The normalized burn ratio (NBR) index was calculated by MTBS using approximately 1-yr prefire and postfire images from the near-infrared (NIR) and shortwave infrared (SWIR) bands of the Landsat sensors, with reflectance values scaled to between 0 and 10000 NBR units:
Pre- and postfire NBR images were next differenced for each Landsat scene pair to generate the dNBR severity classes.
3.2. MODIS vegetation index time series
NASA’s Moderate Resolution Imaging Spectroradiometer (MODIS) satellite sensors Terra and Aqua have been used to generate a 250-m-resolution NDVI (MOD13) global product on 16-day intervals since the year 2000 (Huete et al. 2002; Didan et al. 2016; Shao et al. 2016). The MODIS Collection 6 NDVI dataset provides consistent spatial and temporal profiles of vegetation canopy greenness according to the equation
where NIR is the reflectance in the spectral band of wavelengths from 0.7 to 1.0 μm and Red is the reflectance from 0.6 to 0.7 μm, with values scaled to between 0 and 10000 NDVI units to preserve decimal places in integer file storage. Low values of NDVI (near 0) indicate barren land cover, whereas high values of NDVI (above 8000) indicate dense canopy greenness cover (Potter 2018).
The MOD13 250-m vegetation indices (VIs) have been retrieved from daily, atmosphere-corrected, bidirectional surface reflectance. These MOD13 datasets were downloaded from the files available online (https://modis.gsfc.nasa.gov/data/dataprod/mod13.php; NASA LP DAAC 2007) for time series analysis across the study area. The VIs were computed from MODIS-specific compositing methods based on product quality assurance metrics to remove all low-quality pixels from the final NDVI value reported. Low-quality snow-covered, cloudy, and water pixels were identified and excluded using other MODIS atmospheric data masks. From the remaining good-quality growing-season (May–October) NDVI values, a constrained view-angle approach (closest to nadir) then selected the optimal pixel value to represent each 16-day compositing period.
3.3. Elevation and land-cover map layers
Digital elevation (in vertical meters) for Alaska was derived from USGS (2016) mapping at 300-m ground resolution. Vegetation cover was mapped for the state at 30-m ground resolution from the 2011 National Land Cover Dataset (NLCD) of Alaska (Homer et al. 2015; Selkowitz and Stehman 2011). The overall thematic accuracy for the Alaska NLCD was 76% at Level II (12 classes evaluated).
3.4. Statistical analysis methods
The BFAST methodology was applied to MODIS NDVI monthly time series. BFAST was developed by Verbesselt et al. (2010a,b) for detecting and characterizing abrupt changes within a time series, while also adjusting for regular seasonal cycles. A harmonic seasonal model is first applied in BFAST to account for regular seasonal phenological variations. BFAST next computes the ordinary least squares moving sum (OLS-MOSUM) by considering the moving sums of the residuals after the harmonic seasonal model has been removed from the time series data values. MOSUM tests for structural change using the null hypothesis that all regression coefficients are equal; that is, every observed value can be expressed as a linear function with the same slope (Zeileis et al. 2002). If the null hypothesis is true, the values can be modeled by one line with that slope and the sum of residuals will have a zero mean. MOSUM compares moving sums of residuals to test the likelihood of the regression coefficient for a certain time period based on a user’s input stating the minimum time between potential “breakpoints.” A rejection of the null hypothesis indicates that the regression coefficient changes at that point in time.
The MOSUM uses a default p value of 0.05, meaning that the probability of it detecting a structural change when none has occurred is less than 5%. If MOSUM does not detect some structural change with a confidence level of 95%, it returns a “no breakpoints” result. If MOSUM detects some structural change with a confidence level of 95%, it then processes the time series through a second test, which is used to determine where the breakpoints are located in time. The optimal number of breakpoints is estimated by minimizing the Bayesian information criterion (BIC) and the monthly timing of a breakpoint is estimated by minimizing the residual sum of squares of this regression. The output of this function is a 95% confidence interval for each breakpoint (expressed as two date numbers that define a range).
For BFAST time series analysis, MOD13 NDVI data values (2000–18) were subsampled to include only the growing-season values, during the low snow-cover period from 1 May to 1 October, leaving about 10 observations per year. If a “no data” value was present in the growing-season MOD13 record, then the NDVI from the previous 16-day period was substituted.
4.1. Breakpoint frequency and locations
NDVI time series analysis detected at least one breakpoint at 19% and 22% of the 250-m-resolution MODIS pixels within the Coleen River and the Porcupine River basins, respectively, in the southern portion of the ANWR (Figure 3). The majority of the locations with more than two breakpoints over the entire MODIS NDVI time period were clustered mainly within boundaries of the Boulder Creek Fire (2004), the Rabbit Mountain Fire (2005), the Coleen Mountain Fire (2005), the Coleen River Fires (2006 and 2007), the Rock Slough Fire (2009), and the Sheenjek Fire (2009).
Far fewer pixel locations were detected with one or more NDVI breakpoints within all the other river basins of the ANWR, from the Brooks Range north to the coastal plain, and particularly within the Porcupine caribou herd’s spring migration route (Figure 3). The distribution of BFAST outputs showed that less than 1% of the 250-m MODIS pixels in these seven combined northern ANWR river drainages were detected with one or more NDVI breakpoints. The spatial patterns of breakpoint locations followed closely along the main river valleys of the Chandalar, Sheenjek, Camden Bay, and Firth Rivers, generally below 1000-m in elevation.
4.2. BFAST output sample results
Structural change in the NDVI time series from BFAST analysis at four selected locations (Figure 4; locations shown in Figure 3) were plotted to illustrate breakpoint detection within large wildfire areas that burned since the year 2000 (according to MTBS fire boundary records). All four results detected an NDVI breakpoint at the MTBS-reported year of wildfire ignition. Both BFAST results for the Stringo Lake Fire of 2002 (in the Chandalar River basin) and for the Rabbit Mountain Fire of 2005 (in the Coleen River basin) showed the postfire slope of the NDVI (Tt) fitted “trend” component to be strongly positive until the end of the time series in early 2018. In the case of the Stringo Lakes Fire location, the postfire trend recovered about 2000 NDVI units after the fire event, to approximately prefire NDVI levels. Two additional burned area examples from within the Porcupine Flats and Graphic Lakes subbasins displayed BFAST results for locations that were detected with two breakpoints, and showed that the NDVI trend (Tt) component declined slightly from the prefire mean NDVI level over the entire 18-yr time series.
The outputs for BFAST’s “noise” (et) component (Figure 4, bottom portion of each panel) identified the dates and magnitudes of the largest negative residuals (LNR) from the deseasonalized and detrended NDVI record. These LNR dates for all four of the selected locations were consistently clustered within the years 2012 and 2015, and commonly in early growing-season months. The magnitudes of these LNRs were consistently between −3000 and −4000 NDVI units. In all these cases, the dates of LNR did not correspond to the MTBS-recorded dates of large wildfires, nor to the breakpoints detected in the BFAST analysis.
Within major ANWR river basins along the Porcupine caribou herd’s spring migration route, most breakpoint dates were detected in the years 2003, 2006, and 2017 (Figure 5). In contrast, the majority of LNR dates for locations detected with at least one breakpoint were identified during the years 2001 and 2005 (Figures 5 and 6). The overall 18-yr trend in NDVI for locations detected with at least one breakpoint was skewed strongly to negative (Tt) slope values, with a majority having negative trends between −500 and −1500 NDVI units over the time series (Figure 5).
4.3. Greening and browning trends
For the major ANWR river basins along the Porcupine caribou herd’s spring migration route, BFAST could detect positive (greening) NDVI trends between the years 2000 and 2018 at less than 7% of the total MODIS 250-m pixel area coverage (Figure 7). Most of these greening NDVI trend locations were detected along the East Fork of the Chandalar River valley and its tributaries below 1000-m elevation.
Conversely, the vast majority of MODIS time series data in the Brooks Range and coastal plain ecoregions showed no detectable change in NDVI over the 18-yr study period (Figure 7). BFAST results at four selected locations with no detectable change in NDVI were plotted (Figure 8) to show that NDVI mean annual values (determined from the Tt component) at elevations greater than 1400 m in the Brooks Range were all lower than 1000 NDVI units. Owing to a relatively low NDVI seasonal amplitude of less than 500 units (determined from the St seasonal component), there was no structural change detected for such sparsely vegetated locations of these northern ANWR ecoregions.
Among the small number of locations with no breakpoints detected but with some structural change in the time series, NDVI trends (determined from the Tt component) along the Porcupine caribou herd’s spring migration route (Ryder et al. 2007) were skewed toward positive slope values (Figure 9). The distribution of NDVI trends showed that the majority of these slope values were between +200 and +500 NDVI units over the 18-yr time series. The timing of LNR values from BFAST results for pixels with no NDVI breakpoints detected along the Porcupine caribou herd’s spring migration route showed that the years 2001, 2008, 2011, and 2013 had the highest numbers of single negative NDVI anomalies (Figure 9).
To address the seasonal timing of abrupt negative NDVI anomalies and potential shifts in the timing of forage availability for large grazing herbivores in tundra-dominated ecoregions of the ANWR, it was calculated that the majority (39%) of LNR values detected over the entire 18-yr MODIS time series occurred during the month of May, whereas 27% of LNR values detected during the month of September. The months of June, July, and August each had less than 15% of LNR values detected from the full MODIS NDVI record. There was a slight southern-to-northern lag effect in the early growing-season timing of LNR values, with the majority of negative NDVI anomalies (34%) detected in June rather than May (18%) in the northern Kongakutr River basin.
This is the first study to use 18 years of observations and the highest-resolution MODIS NDVI (at 250-m resolution) to generate detailed map results of vegetation change over the ANWR. A region-wide application of BFAST revealed that NDVI breakpoints and negative 18-yr trends in vegetation greenness over the years 2000–18 could be explained in large part by the impacts of large wildfires in the southern boreal forest ecoregions of the Refuge. The vast majority of MODIS time series data in the Brooks Range and coastal plain ecoregions of the ANWR showed no detectable trend in NDVI over the 18-yr study period. Vegetation change patterns that are occurring at scales of several meters on the ground may, however, not be detectable in MODIS NDVI data at 250-m resolution.
The MODIS BFAST results reported here were consistent with those of Jorgenson et al. (2018), who used high-resolution aerial imagery sampled across the ANWR to estimate that 18% of the area had been altered over the past 50 years, mainly by wildfire and postfire succession, shrub and tree increase in the absence of fire, river erosion and deposition, and ice-wedge degradation. At least one MODIS 250-m NDVI breakpoint was likewise detected in around 20% of the MODIS pixels within both the Porcupine River and Coleen River basins of the ANWR, southern ecoregions that made up the vast majority of change in green vegetation cover over the study area.
These MODIS BFAST results were also consistent with those of Pattison et al. (2015), who reported on changes in species composition and NDVI for field plots located on the coastal plain in the ANWR from 1984 to 2009 and found few changes in plant community composition and no detectable increases in total vegetative cover. Positive greening trends were detected in the MODIS NDVI data since the year 2000 in less than 7% of the total MODIS 250-m pixel area coverage of the tundra-dominated ecoregions of the ANWR.
It is worth noting that Pattison et al. (2015) have postulated that coarser AVHRR satellite data (at 1-km to several-kilometer resolution) may be demonstrating increasing trends in NDVI that are not occurring on the ground in tundra ecosystems of northeastern Alaska. The positive (but unvalidated) trend in AVHRR NDVI observed at long-term coastal plain plots by Pattison et al. (2015) from 1988 to 2006 was similar to the positive trends in NDVI seen in other studies in the Arctic using the AVHRR data (Goetz et al. 2005; Verbyla 2008; Beck and Goetz 2011).
Forkel et al. (2013) used BFAST with the third-generation Global Inventory Monitoring and Modeling System (GIMMS) NDVI (NDVI3g) dataset AVHRR data for change analysis across the state of Alaska. In comparison to a previous version of the GIMMS dataset, this NDVI time series was improved for applications in high-latitude regions through calibrations to stable targets in these regions (Kaufmann et al. 2000). The dataset covered the period from July 1981 to December 2011 with 2-weekly temporal resolution and 8-km spatial resolution. Greening trends from the 1980s were detected mostly in northern Arctic regions, although BFAST detected breakpoints in the Alaskan tundra around the year 2000, followed by browning trends. As in the present BFAST study of MODIS 250-m NDVI, these authors reported that most breakpoints in NDVI time series coincided with large wildfire events.
The effects of high burn severity wildfires since the year 2000 on MODIS greenness changes in all wetlands of Alaska were documented by Potter (2018) using BFAST time series analysis. This study showed that burned wetlands generally recovered their green cover seasonal profiles to relatively stable prefire levels in less than 10 years. Negative changes in NDVI the 18-yr MODIS time series were more commonly detected in burned wetland areas after 2013. In years prior to 2013, the NDVI change tended to be positive at burned wetland elevations lower than 400 m, whereas higher-elevation burned wetland locations showed relatively few positive greenness recovery changes by the year 2017.
The magnitude and timing of LNR values from BFAST analysis of satellite NDVI can provide a new and useful metric of abrupt declines in ecosystem green cover that commonly recover rapidly from whatever agent of disturbance was present at the time of the LNR. The majority of LNRs of NDVI for low-elevation tundra cover in the ANWR were detected during the growing seasons of 2008, 2011, and 2013. Gustine et al. (2017) reported that the growing seasons of 2011–13 were warmer than in 1977 on the Arctic coastal plain and in the Brooks Range. This study found that the day of spring ground thaw (≥0°C) occurred 8 days earlier on average along the Arctic coastal plain and the length of the vegetation growing season was 11 days longer in 2013 than in the 1970s.
This BFAST analysis revealed an 18-yr baseline observation of a majority of LNR values for tundra ecoregions detected during the month of May, and the second highest frequency of LNR values detected during the month of September. One can hypothesize from this evidence that a substantial number of abrupt vegetation browning events have been associated with early spring thawing and elevated levels of surface moisture in low-elevation drainages of the ANWR. Further research is necessary to judge whether these types of anomalies are linked to changing growing-season lengths in the Arctic tundra and to caribou migration being increasingly out of sync with the timing of forage availability.
The majority of vegetation cover in the ANWR extending from the Brooks Range to the coastal plain ecoregions was detected with no (positive or negative) growing-season NDVI trends since the year 2000. Most of the detected breakpoints and strongly negative growing-season NDVI trends were confined to the area within the boundaries of large wildfires in the Coleen River and Porcupine River basins. Using the highest-resolution MODIS data at 250-m resolution and an 18-yr time series to generate detailed map results over this vast wildlife refuge, new insights and metrics of change can be derived from BFAST statistical outputs.
The author thanks two anonymous reviewers for comments on an earlier version of the manuscript.