Changes in vegetation cover and snowmelt timing in the Noatak National Preserve of Northwestern Alaska estimated from MODIS and Landsat satellite image analysis

ABSTRACT 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 Noatak National Preserve in northwestern 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 at least one NDVI abrupt change breakpoint was detected in 25% of the MODIS pixels covering the NOAT. All of the large wildfires mapped by Landsat burn severity classifications within the NOAT since the year 2000 coincided with multiple negative NDVI breakpoints. Results showed that six large drainage basins, all within the Eli River and Kugururok River systems of the western NOAT, had significant correlations between early spring snowmelt and annual area burned. Later-than-average snowmelt dates were strongly associated with a high number of abrupt shifts in greenness cover detected by BFAST. Negative (browning) NDVI trends in the de-seasonalized residuals were detected as significant (p < 0.05) at 15% of the total MODIS 250-m pixel coverage of the NOAT.


Introduction
The arctic tundra region of northwestern Alaska, from the Pacific coast across the interior Brooks Range Mountains, includes several national parks that represent important protected areas for climate change research. This region also stands out as the largest in Alaska where the length of the snow cover season has decreased significantly, with complete snowmelt typically coming two weeks earlier in 2015 than it did in 2001 (Cherry, Zhu, & Kirchner, 2017). Moreover, long-term (1982Moreover, long-term ( -2016 analysis of satellite imagery has shown that tundra vegetation cover has been in decline (i.e. "browning") over extensive areas of western Alaska (Epstein et al., 2016).
When snow cover is present on the tundra, approximately 80-95% of incoming solar radiation is reflected directly back into the atmosphere (Armstrong & Brown, 2008), which favors stable cold air masses and atmospheric temperature inversions. When snow cover is not present, the relatively lower albedo of exposed plants and soils can lead to elevated sensible heating of the atmosphere (Chapin et al., 2005). In addition, early melt results in soils losing the insulating properties offered by snowpacks, and early vegetation growth becomes vulnerable to frost and freeze damage from spring-time low temperature extremes and to a longer season of insect outbreak impacts (Bjerke et al., 2014(Bjerke et al., , 2017Epstein et al., 2016).
The Noatak National Preserve (NOAT) covers 80% of the watershed of the Noatak River, one of the largest drainages in northwestern Alaska. The NOAT also encompasses strong climate gradients, from maritime influences in the western region of the Preserve near the Chukchi Sea, to the strongly continental climate found in the upper reaches of the Noatak River. For the time period 1930 to 2016, Sousanes and Hill (2017) reported that the date of the last freeze in spring in the NOAT is trending earlier and the first freeze in fall is trending later. More recently, Cherry et al. (2017) used MODIS snow cover data for the NOAT to document that spring snowmelt dates in the years 2004 and 2010 were 10-20 days earlier across the Preserve than the mean snowmelt date since the year 2000.
Recent changes in arctic tundra ecosystems may be related to and landscape wetting or drying and snow cover dynamics. For instance, in a study of land cover change using high-resolution aerial imagery for the Arctic National Wildlife Refuge (ANWR) in northeastern Alaska, Jorgenson, Jorgenson, Boldenow, and Orndah (2018) estimated that 18% of the area had been notably altered over the past 50 years, mainly by wildfire and post-fire succession, shrub and tree increase in the absence of fire, river erosion and deposition, and ice-wedge degradation. Boreal forest cover further south showed changes associated with landscape drying and decreases in lake area. 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 Mountains. 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. Pattison and Welker (2014) reported 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. Bokhorst, Bjerke, Street, Callaghan, and Phoenix (2011) reported that extreme (+2 to +10°C), short-lived, winter warming in the arctic can rapidly melt the snowpack and subsequently expose plants to much colder temperatures due to the loss of insulating snow. Such single freezing events have been shown to reduce plant reproduction and increase shoot mortality. Root growth was reduced by more than 25% and gross primary productivity was reduced by more than 50% in the summer growing season following repeated experimental snow loss events.
At a (Alaska) statewide scale, Potter, Li, and Crabtree (2013) mapped tundra greening and browning from 2000 to 2010 from monthly MODIS data, and correlated trends in vegetation dynamics with climate index variations. Correlation analysis showed a notable association between the largest positive slopes in tundra greening with annual warming trends of greater than +2 growing degree days (GDD) per year. Ecoregions that revealed the highest density of negative slope (tundra browning) points were the northeastern Arctic Coastal Plain and the Brooks Range Foothills, the latter extending into the central portions of the NOAT.
Nevertheless, to date, few of the large-scale studies of vegetation greening or browning in Alaska have included comprehensive structural breaks analysis, designed to simultaneously detect all major disturbances that can alter greening trend statistics and the conclusions about gradual change (greening and browing) in vegetation cover density and tundra productivity. 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, Leisch, Hornik, & Kleiber, 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 NOAT 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 and disturbance frequency changed over large areas in the Noatak National Preserve since the year 2000 and particularly in relation to variation in the timing (early versus late season) of annual snowmelt?" Statistical analysis of changes in the MODIS NDVI time series was conducted using the "Breaks for Additive Seasonal and Trend" method (Verbesselt, Hyndman, Newnham, & Culvenor, 2010a;Verbesselt, Hyndman, Zeileis, & Culvenor, 2010b;Potter, 2018). Landsat satellite images at 30-m resolution, acquired during summer (July-August) months, were used to validate interannual MODIS NDVI changes. This is the first study to use the entire 18-yr MODIS NDVI record to capture changes in vegetation cover over the NOAT and the western Brooks Range Mountains, and the first to This is the first study to use the entire 18-yr MODIS NDVI record to capture changes in vegetation cover over the NOAT and the western Brooks Range Mountains, and the first to relate inter-annual variations in snowmelt timing to the frequency and severity of disturbances of all kinds in this region.

Study area
The NOAT drains 33,670 km 2 of the western Brooks Range just above the Arctic Circle ( Figure 1) and contains a range of lowland and alpine tundra, shrub tundra, and treeline habitats. Vegetation cover in the Preserve is predominantly Eriophorum tussock-shrub tundra. Shrub communities consist mainly of Salix spp., Betula glandulosa, and Alnus crispa. Boreal forests of spruce (Picea marina and P. glauca) and birch (Betula neoalaskana) are sometimes present at low elevations (Swanson, 2017).
The Noatak River flows across the Preserve from east to west. The northern boundary of the NOAT follows the crest of the DeLong Mountains and the southern border crosses the Baird Mountains, forming the western edge of the Brooks Range. Elevation ranges from around sea level at the mouth of the Noatak River to about 1750 m in the Schwatka Mountains.
Annual precipitation in the NOAT averages between 40-70 cm, with the majority falling as rain during the summer and fall months of July to September. The mean annual temperature at the coastal site of Kotzebue (at the southwestern margin of the NOAT) is −5.1°C, ranging from a July average of 13°C to a January average of −19°C (Sousanes & Hill, 2017). The months of June to September generally have average temperatures above freezing. Temperatures are cooler at the higher elevation Brooks Range sites of the NOAT in the summer than at the low elevation location of Kotzebue. Over the available 70-year weather station record , mean annual temperature in the NOAT has increased significantly at +2.1°C, and particularly in the winter months, by a significant increase of +4.3°C (Sousanes & Hill, 2017).
Wildfires have burned more than 900 km 2 of tundra in the NOAT since the year 2000 (Eidenshink et al., 2007), with the vast majority of area burned in the years 2005, 2010, and 2012 and started during the months of June or July. The largest single fires recorded were the Uvgoon Creek Fire of 2012 and the Kaluktavik River and Sidik Lake fires of 2010.

MODIS vegetation index time series
NASA's MODIS (Moderate Resolution Imaging Spectroradiometer) 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 (Didan, Munoz, Solano, & Huete, 2016;Shao, Lunetta, Wheeler, Iiames, & Campbell, 2016). The MODIS Collection 6 NDVI data set provides consistent spatial and temporal profiles of vegetation canopy greenness according to the equation: where NIR is the reflectance 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 10,000 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.
The MOD13 250-m vegetation indices (VIs) have been retrieved from daily, atmosphere-corrected, bidirectional surface reflectance (LP-DACC, 2007). These MOD13 data sets were downloaded from the files available at modis.gsfc.nasa.gov/data/dataprod/mod13.php 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. Snow-covered, cloud and water pixels were identified and excluded using other MODIS atmospheric data masks. From the remaining good-quality growing season (May to October) NDVI values, a constrained view-angle approach (closest to nadir) then selected the optimal pixel value to represent each 16-day compositing period.

MODIS snowmelt timing products
Snowmelt timing maps (STMs) for North America were obtained from the database developed by O'Leary, Hall, Medler, Matthews, and Flower (2017), based on the MODIS 8-day composite snow-cover product (MOD10A2) and covering the time period from 2001 to 2015. The MOD10A2 product is derived from the Normalized Difference Snow Index (NDSI) using a 50% snow/no-snow classification threshold. Snowmelt timing (for the first no-snow Julian day of the year) was defined as a snow-free reading following two consecutive snow-present readings for a given 500-m STM pixel. STM products have been validated using in-situ observations and comparisons with SNOTEL meteorological stations from across the western United States (Natural Resources Conservation Service (NRCS), 2016). O' Leary, Hall, Medler, and Flower (2018) reported that differences between MODIS STMs and the available SNOTEL stations had a mean error of 4.3 days, which was less than the 8-day temporal resolution of the MOD10A2 product and therefore indicated that their STM algorithm captures on-the-ground snowmelt timing accurately.
Following the approach of O'Leary, Bloom, Smith, Zemp, and Medler (2016) for comparisons of snow-free dates across drainage basins to disturbance areas and numbers, the dates of STM were transformed into relative values by Z-score calculation (Aho, 2014) of number of standard deviations from the mean over the years 2001 to 2015. A negative Z-score indicated a year with an earlier snowmelt relative to the multi-year mean.
Landsat NDVI and fire boundary products NDVI images from Landsat sensors was selected from the United States Geological Survey (USGS) Earth Explorer data portal (available online at earthexplorer.usgs.gov) for the years 2004 to 2018. The most cloud-free Landsat scenes over the NOAT were acquired between July 1 and August 31 each year. All images used in this study were geometrically registered (UTM Zone 2 -CHECK) using terrain correction algorithms (Level 1T) applied by the USGS EROS Data Center and cloud-filtered using the pixel quality layer derived for each scene and date.
For the Landsat 4-5 Thematic Mapper (TM) images acquired between 1985 and 2011, 30-m resolution surface reflectance data were generated from the Landsat Ecosystem Disturbance Adaptive Processing System (Masek et al., 2006). Moderate Resolution Imaging Spectroradiometer (MODIS) atmospheric correction routines were applied to Level-1 TM data products. Water vapor, ozone, geopotential height, aerosol optical thickness, and digital elevation are input with Landsat data to the Second Simulation of a Satellite Signal in the Solar Spectrum (6S) radiative transfer models to generate top of atmosphere (TOA) reflectance, surface reflectance, brightness temperature, and masks for clouds, cloud shadows, adjacent clouds, land, snow, ice, and water. Landsat 8 (after 2012) surface reflectance products were generated from the L8SR algorithm, a method that uses the scene center for the sun angle calculation and then hard-codes the view zenith angle to 0. The solar zenith and view zenith angles are used for calculations as part of the atmospheric correction.
Digital maps of burn severity classes (low, moderate, and high) 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 one-year pre-fire and post-fire images from the near infra-red (NIR) and shortwave infra-red (SWIR) bands of the Landsat sensors, with reflectance values scaled to between 0 and 10,000 NBR units.
Pre-and post-fire NBR images were next differenced for each Landsat scene pair to generate the dNBR severity classes. All available, relatively cloud-free (< 30%), Landsat NDVI images from surface reflectance processing were obtained for the mid-summer months (July and August) from the years 2000 to 2018 from the USGS Earth Explorer image service (online at earthexplorer. usgs.gov). Path/rows 81/12, 80/12, and 79/12 covered most of the NOAT study area. NDVI valuesat 30-m resolution in selected areas of the NOAT, where cloud-free Landsat images overlapped, were compared before and during years of extreme early or extreme late snow melt dates, as derived from the MODIS products from of O'Leary et al. (2017).

Elevation, land cover, and drainage basin map layers
Digital elevation (in vertical meters) for the NOAT was derived from USGS 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 & Stehman, 2011; available at www. mrlc.gov/nlcd11_leg.php). The overall thematic accuracy for the Alaska NLCD was 76% at Level II (12 classes evaluated). Drainage basin boundaries and names at the HU12 level for the study area were obtained from the USGS hydrologic unit (HU) database (Seaber, Kapinos, & Knapp, 1987).

Statistical analysis methods
The BFAST (Breaks for Additive Seasonal and Trend) methodology was applied to MODIS NDVI monthly time series. BFAST was developed by Verbesselt et al. (2010a), (2010b)) 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 i.e. 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 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 to 2018) were subsampled to include only the growing season values, during the low snow cover period of May 1 to October 1, leaving about 10 observations per year.
The 18-year trends in de-seasonalized NDVI values from BFAST outputs at each 250-m MODIS pixel were calculated using the Sen's slope test for linear rate of change (Sen, 1968). The Sen's test is a method for robustly fitting a correlation line to sample points in the linear regression plane by choosing the median of the slopes of all lines through pairs of points. The test is relatively insensitive to outliers and can be significantly more accurate than non-robust simple linear regression (least squares) for skewed data sets.

Breakpoint frequency and locations
NDVI time series analysis with BFAST detected at least one breakpoint between 2000 and 2018 at approximately 25% of the 250-m resolution MODIS pixels within the NOAT (Figure 2). Out of all pixel locations where at least one breakpoint was detected, about 39% if these locations were also detected with two breakpoints and 10% of these locations were detected with three breakpoints. More than half (56%) of all NDVI breakpoints were detected as abrupt negative shifts in green vegetation cover. It should be noted however that about 13% of all pixel locations in the NOAT did not have a sufficient percentage of cloud-free NDVI values in the time series to carry out a valid structural change analysis.
All of the large wildfires mapped by the MTBS within the NOAT since the year 2000 coincided with multiple negative NDVI breakpoints ( Figure 2). The fires with the highest frequency of negative breakpoints during the mid-summer period (July and August) were the Uvgoon Creek Fire (2012) Figure 5). The distribution of NDVI trends showed that the majority of positive slope values were between +20 to +80 NDVI units per year over the entire time series. Generally speaking, greening (or browning) trends that were significant at p < 0.05 had a change of greater than +50 NDVI units per year. The greatest concentration of NDVI greening locations were detected in the western NOAT, within the Eli River, Deadlock Mountain, and Kivivik Creek-Noatak River valleys, commonly below 200-m elevation. Other drainage basins with extensive and significant greening trends were detected in the Aliktongnak Lake and Kugururok River valleys of the western NOAT.
Locations that showed a relatively high positive NDVI anomalies during the first few years of the 18-year time series, followed by a notably lower variance in anomalies later in the time series, were commonly detected by the Sen's test with a negative (browning) slope. This type of non-significant browning trend more accurately reflects a decrease in interannual variability in NDVI over about the past five years of the time series (2013)(2014)(2015)(2016)(2017)(2018). At these locations, there has been a stabilization of the growing season NDVI phenology, after earlier years with repeated departures in the positive greening direction. The histogram distribution of highest positive residual (HPR) occurrences for all MODIS 250-pixels in the NOAT in the 18-year time series confirmed that the years 2000 and 2003 had the highest numbers of extreme positive NDVI anomalies (Figure 6(a)).
Negative (browning) NDVI trends in the deseasonalized residuals over the years 2000 and 2018 were detected as significant (p < 0.05) at 15% of the total MODIS 250-m pixel coverage of the NOAT ( Figure 5). The distribution of NDVI trends showed that the majority of negative slope values were between −30 to −80 NDVI units per year over the entire time series. The greatest concentration of these browning NDVI locations were detected within the Nimiuktuk, Noatuk, Makpik Creek, Akiknnaak Peaks-Noatak River, and the Rough Mountain Creek drainages of the northern and eastern NOAT, mainly above 400-m elevation.

Lowest negative residuals
The timing of LNR occurrences from the BFAST results (for pixels with no NDVI breakpoints detected) showed that the year 2001 had the highest numbers of extreme negative NDVI anomalies, followed by 2005, 2007 and 2009, each with about 60% fewer LNR occurrences than 2001 (Figure 6(b)). The majority (55%) of LNR occurrences detected over the entire 18-yr MODIS time series occurred during the month of May, and another 32% of LNR occurrences were detected during the month of June. The months of July, August, and September each had between 8% and 17% of LNR occurrences detected from the full MODIS NDVI record.

STM correlations with burned areas and breakpoint totals
Based on the date of snow-free cover each year derived from MODIS STM values and averaged over the entire NOAT (Figure 7 As a single-year example, the relative (Z-score) timing of snow-free cover shown in Figure 8 for 2005 revealed high variability across the NOAT, with the earliest snow-free dates in drainage basins of the central Noatak River basin, particularly the Kugururok River and Sapun Creek drainages, and in the Noatuk and Middle Anisak River drainages farther to east. The latest snow-free dates for this year were detected in drainage basins of the southern Agashashok River and Culter River sections of the NOAT.
Correlations over the years 2001 to 2015 of the average STM snow-free date with moderate burn severity (MBS) and high burn severity (HBS) area totals each year from the MTBS Landsat record of large fires showed that six drainage basins, all within the Eli River and Kugururok River systems of the western NOAT (Figure 9), had significant (p < 0.1) linear regression coefficients of between −0.4 and −0.5 (r values). In those drainage basins with the highest negative correlation of MBS area totals with STM snow-free dates it could be inferred that years with earlier than average snowmelt (i.e. in 2004, 2005, 2010, and 2015 at about 14 days earlier) was strongly associated with the high annual area burned at MBS. One drainage basin, Kangilipak Lake in the eastern Noatak River basin, showed a positive correlation (p < 0.05) of snow-free dates with MBS area. We found no drainage basins in the NOAT that had a significant (at p < 0.1) linear regression of HBS area versus average STM snow-free date.
Correlations over the years 2001 to 2015 of the average STM snow-free date with the total number of NDVI breakpoints detected by BFAST analysis each year (Figure 10) showed 96 drainage basins (comprising roughly one-third of the NOAT area) had significant (p < 0.05) linear regression coefficients of between +0.38 and +0.87 (r values). Only two drainage basins (out of 286 total), both located in the extreme southern outlet of the Noatak River at Shiliak Creek, showed a significant negative correlation of average STM date with the number of NDVI breakpoints detected each year. Those drainage basins with the highest positive correlation (r > 0.75) of NDVI breakpoints with the STM snow-free date were in the Kaluich Creek, Makpik Creek, Imelyak River, and Amakomanak Creek drainages of the eastern and southern NOAT. It could be inferred that in these drainage areas, later-than- average snowmelt dates were strongly associated with a high number of abrupt shifts in greenness cover. In the drainage areas highlighted in Figure 10 with high positive correlation values, late snow cover areas in years (such as 2003 and 2006) that persisted into mid-June was commonly detected as an anomalously low NDVI value for that time of year and counted as a breakpoint.

Landsat NDVI comparisons
Among all the relatively cloud-free Landsat images acquired since the year 2000 over the NOAT, with which one could make comparisons of NDVI just before and during years of extreme snowmelt dates, only those from path/row 80/12 on August 5 of 2014 and path/row 81/12 on July 30 of 2015 were significantly overlapping and within one week of date acquisition, one year apart.
Comparison of MODIS (500-m resolution) snowmelt timing (as the snow-free day of the year) in 2015 with the difference in mid-summer NDVI from Landsat images acquired in 2014 and 2015 ( Figure 11) showed that declines in NDVI of around −1000 units were commonly detected in areas of relatively early 2015 snow-free dates (early May), specifically in low elevation (below 300 m) parts of the Kagvik Creek, Nunaviksak Creek, Trail Creek, and Kaluktavik River drainages. Areas such as the Wrench Creek drainage showed increases in NDVI of around +1000 units and were commonly detected in areas with relatively late 2015 snow-free dates (late May). Highelevations areas at between 600 m and 900 m showed the latest snow-free dates (i.e. mid-June) and generally little difference in mid-summer NDVI between 2014 and 2015.

Discussion
Changes in tundra vegetation can have important effects on fire regimes (Racine, Johnson, & Viereck, 1987), permafrost stability, surface water and energy fluxes (e.g. Frost, Epstein, Walker, Matyshak, & Ermokhina, 2017;Kępski et al., 2017), as well as the abundance and distribution of wildlife populations (e.g. Fauchald, Park, Tømmervik, Myneni, & Helene Hauser, 2017;Horstkotte et al., 2017). 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 a large national preserve in western arctic Alaska. The NOAT covers a wide range of climate regimes, from lowland coastal to upland mountain, making it a uniquely valuable study area for vegetation dynamics in the Alaskan tundra. The western tundra region of Alaska is climatically distinct from the North Slope, in that it receives 2-3 times more annual precipitation, is an average of 5-10°C warmer than the North Slope, and is subject to synoptic patterns more conducive to lightning during fire season (Shulski & Wendler, 2007).  Results from NDVI time series analysis using BFAST showed that the years 2001 to 2003 had the highest numbers of extreme positive NDVI anomalies and also had the consistently latest (into early June) average snow-free dates over the NOAT. These findings together support the hypothesis that late snowmelt years in this region of the Alaskan arctic are associated with increases in that same year's mid-summer green cover, despite a typical 4-6 week delays in the early season (May-June) green-up period.
Conversely, previous studies (Bokhorst et al., 2011;Epstein et al., 2016) have reported that wintertime warming in the arctic tundra can rapidly melt snowpacks and begin to expose new plant growth to freezedamaging temperatures due to the loss of insulating snow. Moreover, early snow-melt results in soils drying out sooner and losing the insulating properties  offered by snowpacks (Bjerke et al., 2014(Bjerke et al., , 2017. Results from this analysis of NDVI breakpoints and anomalous values indicated that the transition in snow cover timing from the extreme late snowmelt year of 2003 to the extreme early snowmelt years of 2004 and 2005 had a significant impact on vegetation green cover within the Nimiuktuk River, Noatuk River, Middle Anisak River, and Makpik Creek drainages of the NOAT. The pervasive negative NDVI slope trends over the entire 18-year MODIS time series within these drainage basins, together with several years of extreme early snow-free dates, is evidence that these areas were most likely to have been subjected to vegetation cover browning by a frost-drought effect. Between 2000 and 2016, 22 large wildfires have occurred within the NOAT burning over 914 km 2 , with more than half of the total area burned occurring in 2010 (MTBS; Eidenshink, 2007). BFAST results detected all of these large wildfire events as significant negative NDVI breakpoint clusters within the study area. The majority of NDVI negative breakpoint dates in the NOAT were detected in the mid-summer (July-August) period in the high-fire frequency years of 2005 and 2015.
The results from correlation analysis of annual snowmelt timing Z-scores with the total area burned within major drainage basins of the NOAT was identical to that of O' Leary et al. (2016) from their analysis of the burned areas in the Northern Rocky Mountain region of Wyoming, North Dakota, Montana, and Idaho. This earlier study reported a significant correlation (p < 0.05) between early spring snowmelt and total annual area burned. Our correlation analysis for the NOAT is the first of its kind for tundra areas of Alaska, and went a step further and used burn severity categories (MBS and HBS) and their area totals each year from the MTBS Landsat record of wildfires. The results showed that six large drainage basins, all within the Eli River and Kugururok River systems of the western NOAT, had significant correlations between early spring snowmelt and annual area burned. As pointed out by O'Leary et al. (2016), early snowmelt is one way to achieve low fuel moisture and extended fire seasons leading to increased annual area burned.
Considering longer timeframes for fire returns, over the 60-year period of 1950 to 2010, Hu et al. (2015) reported that 14% of the tundra has burned in both the Noatak and Seward Peninsula ecoregions, corresponding to fire rotation periods of around 420 years. Six thousand years of fire history in the NOAT was examined by Higuera, Barnes, Chipman, Urban, and Hu (2011), who suggested that recent changes in vegetation may be linked to changes in fire occurrence. Tundra vegetation has responded to local-scale climate, which in turn has influenced the flammability of surrounding areas. New directions of these impacts may depend upon the specific makeup of future tundra vegetation, especially increasing densities of shrubs and deciduous trees.
The magnitude and timing of negative residual anomalies 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 tundra cover in the NOAT were detected during the growing seasons of 2001, 2005, 2007, and 2009. Example BFAST plots have shown that the magnitude of these LNR departures was typically between −2000 and −3000 NDVI units, which matched or exceeded the magnitude of the NDVI decline resulting from large wildfire burns in the NOAT.
BFAST analysis further revealed that a majority of LNR values for tundra cover of the NOAT were detected during the month of May, and the second highest frequency of LNR values detected during the month of June. One can hypothesize from this evidence that a substantial number of abrupt vegetation browning events have been associated with early spring thawing (such as in 2005) in low elevation drainages of the NOAT and with late spring snow cover (such as in 2001 and 2009) at higher elevation mountain areas. Further research is necessary to judge whether these types of anomalies are linked to changing growing season lengths in the Arctic tundra, and to more frequent frost-drought damage events.

Conclusions
Using the highest resolution MODIS data at 250-m resolution and an 18-yr time series to generate detailed map results over this vast protected area, new insights and metrics of change can be derived from BFAST statistical outputs. Results from NDVI time-series analysis showed that late snow cover area, in years such as 2001, 2003, and 2006, that persisted into mid-June was commonly detected with an anomalously low NDVI value for that time of year and counted as a breakpoint. The transition from an extreme late snowmelt year to an extreme early snowmelt year had significant impact on vegetation green cover within several northern river drainages of the NOAT. Negative NDVI slope trends over the entire 18 year MODIS time series within such drainage basins was an indicator that vegetation browning occurred by a frost-drought effect, but requires further investigation to verify this type of widespread tundra vegetation change.

Disclosure statement
No potential conflict of interest was reported by the authors.