Changes in vegetation phenology and productivity in Alaska over the past two decades

Understanding trends in vegetation phenology and growing season productivity at a regional scale is important for global change studies, particularly as linkages can be made between climate shifts and the vegetation’s potential to sequester or release carbon into the atmosphere. Trends and geographic patterns of change in vegetation growth and phenology from the MODerate resolution Imaging Spectroradiometer (MODIS) satellite data sets were analyzed for the state of Alaska over the period 2000 to 2018. Phenology metrics derived from the MODIS Normalized Difference Vegetation Index (NDVI) time-series at 250 m resolution tracked changes in the total integrated greenness cover (TIN), maximum annual NDVI (MAXN), and start of the season timing (SOST) date over the past two decades. SOST trends showed significantly earlier seasonal vegetation greening (at more than one day per year) across the northeastern Brooks Range Mountains, on the Yukon-Kuskokwim coastal plain, and in the southern coastal areas of Alaska. TIN and MAXN have increased significantly across the western Arctic Coastal Plain and within the perimeters of most large wildfires of the Interior boreal region that burned since the year 2000, whereas TIN and MAXN have decreased notably in watersheds of Bristol Bay and in the Cook Inlet lowlands of southwestern Alaska, in the same regions where earlier-trending SOST was also detected. Mapping results from this MODIS time-series analysis have identified a new database of localized study locations across Alaska where vegetation phenology has recently shifted notably, and where land cover types and ecosystem processes could be changing rapidly.


Introduction
Recent studies have shown that many regions of Alaska are experiencing a high level of inter-annual variability in the seasonality of vegetation growth [1][2][3]. For example, the northern and interior Brooks Range Mountains have experienced the greatest year-to-year variation in snowmelt timing dates of any sub-region of Alaska over the past 20 years [4]. Farther north, the date of spring ground thaw occurred eight days earlier and the length of the vegetation growing season was shown to be 11 days longer on the Arctic Coastal Plain in 2013 compared to the 1970s [5].
High magnitudes of inter-annual variability in springtime snow-free conditions may be having strong impacts on many aspects of ecosystem functioning [6], including the dynamics of Alaskan wildlife populations [7] and on the area of tundra severely burned each year by wildfires [8]. Over the past 30 years, the expansion of woody shrub cover in northern Alaska has been reported from both remote sensing and field studies [9][10][11][12][13]. In some boreal forest sub-regions of the state, evergreen tree cover has decreased, while deciduous tree cover has increased, along with the expansion of shrub and herbaceous vegetation above tree lines [13], patterns which have been confirmed by field measurements of shrub encroachment and greening trends [12].
Satellite remote sensing of land-surface phenology can be used to characterize large-scale variations in the vegetation growing season from one year to the next [14,15]. Across Alaska, integrated monthly "greening", observed as positive trends in the satellite-derived normalized difference vegetation index (NDVI), can indicate increased ecosystem productivity and/or higher seasonal density of vegetation cover, while negative integrated NDVI trends or "browning" have been construed as reduced productivity in vegetation cover due to wildfire, insect damage, drought stress, and/or thawing permafrost [8,[16][17][18][19][20][21].
The time-integrated NDVI (TIN) phenology metric [22] can be treated as a measure of seasonal plant growth because it represents the daily interpolation of green cover detected by the MODIS satellite sensors above the baseline for the entire duration of the growing season. As such, shifts in TIN can serve as proxies for change in annual primary productivity of an ecosystem from year to year [23][24][25][26]. Moreover, tundra vegetation biomass (shrubs, forbs, graminoids, lichens, and mosses) in northern Alaska has been found to correlate closely (R 2 = 0.8) and increase exponentially with maximum summer NDVI (MAXN) [27]. Comparison of TIN to the MAXN level can reveal important attributes of a seasonal growth curve.
The start of the season timing (SOST) date can serve as a proxy for the early-season initiation of vegetation growth in relation to snow-free conditions and thawing of rooting zone active layers [28]. Satellite data analyses have shown significant and widespread advances of SOST over the northern hemisphere during the 1980s and 1990s, trends which were associated with strong springtime air temperature warming over these decades [29,30]. The ongoing advance in SOST in Alaska has been among the highest in the northern latitude zones, averaging greater than three days per decade [31].
In this study, we combined time-series analysis results for the SOST, TIN, and MAXN metrics from the full 19 year Moderate Resolution Imaging Spectroradiometer (MODIS) NDVI phenology record [22] to test the hypothesis that recent vegetation responses to inter-annual variability in springtime snowmelt dynamics and growing season temperatures in Alaska are highly dependent on regional land cover types. Specifically, such responses may be related to different physiological adjustments to local climate and to elevation in boreal forests versus tundra shrub/tussock communities. Digital maps of large wildfire perimeters for Alaska [32] were overlaid on the MODIS time-series results to separate the effects of vegetation burning from all other factors of change in the NDVI metrics.

MODIS Vegetation Index and Phenology Metrics
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 [33,34]. The MODIS Collection 6 [35] NDVI data set provides consistent spatial and temporal profiles of vegetation canopy greenness according to the equation: For this study, we used the data set from Zhu et al. [22] who produced numerous phenology metric products for the state of Alaska from MOD13 NDVI data sets, spanning the years 2000 to 2018. Zhu et al. [22] used the seven-day interval eMODIS data set developed by Reed et al. [36] and the MODIS Quality Layer for calculation of NDVI phenology metrics. In the time series of the pixels, only "good" (cloud free) = 0 and "snow" = 4 data are retained; all other values are interpolated or replaced. The phenology metric algorithms then stack these quality-checked layers into a multiple-band data set, then interpolate and smooth the data using a weighted-least-square algorithm on the interpolated time series.
In areas where relatively low levels of annual TIN have been consistently associated with relatively high levels of MAXN, there is a high probability that vegetation cover has experienced growth limitations or stresses of some kind during a considerable portion of the growing season [37]. Under this assumption, a relatively low ratio value (less than 50:1) of annual TIN to annual MAXN results would indicate greater overall growth stress for the vegetation at a given location.

Sub-Regions, Land Cover, and Elevation Map Layers
The MLRA (Major Land Resource Areas; [38]) product was used to delineate major sub-regions of the Brooks Range Mountains in northern Alaska. Digital elevation (in vertical meters) was derived from United States Geologic Survey (USGS) mapping at 300 m ground resolution. Vegetation cover types were mapped for the state at 30 m ground resolution from the 2011 National Land Cover Dataset (NLCD) of Alaska [39,40] (available at www.mrlc.gov/nlcd11_leg.php). Large river basin boundaries and sub-basin names for the study area were obtained from the USGS hydrologic unit database [41].

Statistical Analysis Methods
Time-series analysis of MODIS annual NDVI phenology metrics over the period 2000 to 2018 was carried out using TrendRaster algorithms developed by Forkel et al. [42]. The significance level (p value) and slope of the trend for each MODIS NDVI time series was estimated by the Mann-Kendall (M-K) nonparametric test [43] using the TrendRaster function of the "green-brown" R software package [44]. The M-K test analyzes the difference in signs between earlier and later data values in a time series. If a trend is present, the sign values will tend to consistently increase or decrease. Every data value is compared to every value preceding it in the time series, which gives a total of n × (n-1) / 2 pairs of data, where "n" is the number of observations in the set. This nonparametric test is applicable for data sets that do not meet the assumption of normality.
For localized correlation analysis between annual SOST and TIN values, we applied a 7 × 7 MODIS pixel (3 km 2 ) moving window linear regression model. This was justified as the smallest local area over which an adequate number of MODIS pixels (N = 49) would be included in the correlation coefficient determination. For regional-scale linear regressions of SOST or TIN values with elevation or forest cover, we expanded the moving window to a 51 × 51 MODIS pixel (162 km 2 ) area in order to realistically capture the scale of mountainous topography and natural variations in vegetation classes in Alaska.

Start of the Growing Season Date
We found that trends in the SOST over the period of 2000 to 2018 showed significantly (p < 0.05) earlier seasonal vegetation greening (at more than one day per year) across sections of the northeastern region of the Brooks Range Mountains into the Richardson Mountains and Mackenzie River Delta of Canada, on the Yukon-Kuskokwim coastal plain, and most notably in the Ahklun Mountains across to Bristol Bay and the Cook Inlet lowlands in the southern coastal areas (Figure 1a,b). Conversely, a trend of significantly later vegetation greening (at more than one day per year) was detected mainly within areas burned by wildfires since 1999 in the Yukon-Kuskokwim and Interior Alaska Highlands. Notable areas outside of fire boundaries that showed significantly later vegetation greening since the year 2000 were detected in the Porcupine River, Tanana River, and Lower Kuskokwim River valleys.

Growing Season Green Vegetation Cover
Trends in both TIN ( Figure 1c) and MAXN (map not shown) over the period of 2000 to 2018 were similar, and both showed the same significantly (p < 0.05) increasing summer vegetation greening in nearly all boreal forest areas burned by wildfires in the Interior Alaska and Yukon-Kuskokwim Highlands (Figure 1d), commonly at rates of greater than TIN = +0.5 units per year. Outside of wildfire boundaries, equivalent vegetation greening rates and distribution were detected in unburned lowland areas of the Seward Peninsula, the Noatak, Colville, Kobuk, Chandalar, and Sheenjek River valleys, and across much of the Arctic Coastal Plain.
In contrast, there was a notable collection of unburned locations (as indicated in Figure 1d) with significantly (p < 0.05) decreasing TIN and MAXN trends over the period of 2000 to 2018, mainly detected in the Ahklun Mountains, from the watersheds of Bristol Bay southward along the Alaska Peninsula, and in the Cook Inlet lowlands of southwestern Alaska, in the same regions where earlier-trending SOST was also detected. There were also isolated areas of decreasing TIN and MAXN detected in the Yukon Flats and the Koyukuk Flats of the Interior Alaska Lowlands that showed no recent (since 1999) records of large wildfires.
We found that a high negative correlation between SOST and TIN in any given year was most frequently detected on the Arctic Coastal Plain, in the southern Brooks Range Mountains, the Interior Alaska Mountain Range, and the Ahklun Mountains, commonly below 1500 m elevation (Figure 2a), indicating that an earlier start of the growing season resulted in higher seasonal green cover and plant production. Conversely, locations where a relatively late start of the growing season resulted in higher integrated green cover were most frequently detected in the northern Brooks Range Mountains above 1500 m elevation, in the Tanana River valley, the Yukon-Kuskokwim highlands, the Cook Inlet lowlands, on the Alaska Peninsula, and along the Gulf of Alaska Coast. The variability from year to year in the correlation between SOST and TIN was lowest across the Brooks Range Mountains and other mountainous areas above elevations of 1500 m ( Figure 2b). Elevation explained more than half of the geographic variation (based on a 162 km 2 moving window) in the mean correlation between annual SOST and TIN over just 5.1% of the Alaskan vegetated land area (Figure 2c), whereas forest cover explained more than half of the geographic variation in the mean correlation between annual SOST and TIN over only 2.4% of the Alaskan vegetated land area ( Figure 2d). Any positive associations shown in Figure 2c,d indicated that the higher the elevation or forest cover, the stronger the correlation between annual SOST and TIN in this area, whereas any negative associations shown in Figure 2c,d indicated that the higher the elevation or forest cover, the weaker the correlation between annual SOST and TIN in this area.

Growing Season Green Vegetation Cover
Trends in both TIN ( Figure 1c) and MAXN (map not shown) over the period of 2000 to 2018 were similar, and both showed the same significantly (p < 0.05) increasing summer vegetation greening in nearly all boreal forest areas burned by wildfires in the Interior Alaska and Yukon-Kuskokwim Highlands (Figure 1d), commonly at rates of greater than TIN = +0.5 units per year. Outside of wildfire boundaries, equivalent vegetation greening rates and distribution were detected in unburned lowland areas of the Seward Peninsula, the Noatak, Colville, Kobuk, Chandalar, and Sheenjek River valleys, and across much of the Arctic Coastal Plain.
In contrast, there was a notable collection of unburned locations (as indicated in Figure 1d) with significantly (p < 0.05) decreasing TIN and MAXN trends over the period of 2000 to 2018, mainly detected in the Ahklun Mountains, from the watersheds of Bristol Bay southward along the Alaska Peninsula, and in the Cook Inlet lowlands of southwestern Alaska, in the same regions where earliertrending SOST was also detected. There were also isolated areas of decreasing TIN and MAXN detected in the Yukon Flats and the Koyukuk Flats of the Interior Alaska Lowlands that showed no recent (since 1999) records of large wildfires.
We found that a high negative correlation between SOST and TIN in any given year was most frequently detected on the Arctic Coastal Plain, in the southern Brooks Range Mountains, the Interior Alaska Mountain Range, and the Ahklun Mountains, commonly below 1500 m elevation (Figure 2a), indicating that an earlier start of the growing season resulted in higher seasonal green cover and plant production. Conversely, locations where a relatively late start of the growing season resulted in higher integrated green cover were most frequently detected in the northern Brooks Range Mountains above 1500 m elevation, in the Tanana River valley, the Yukon-Kuskokwim highlands, the Cook Inlet lowlands, on the Alaska Peninsula, and along the Gulf of Alaska Coast. The variability from year to year in the correlation between SOST and TIN was lowest across the Brooks Range Mountains and other mountainous areas above elevations of 1500 m (Figure 2b). Elevation explained more than half of the geographic variation (based on a 162 km 2 moving window) in the mean correlation between annual SOST and TIN over just 5.1% of the Alaskan vegetated land area ( Figure  2c), whereas forest cover explained more than half of the geographic variation in the mean correlation between annual SOST and TIN over only 2.4% of the Alaskan vegetated land area (Figure 2d). Any positive associations shown in Figure 2c and 2d indicated that the higher the elevation or forest cover, the stronger the correlation between annual SOST and TIN in this area, whereas any negative associations shown in Figure 2c and 2d indicated that the higher the elevation or forest cover, the weaker the correlation between annual SOST and TIN in this area.  where forest cover (respectively) accounted for greater than 50% of the geographic variation in the mean correlation between SOST and TIN over the years 2000 to 2018, based on a regional-scale (51 × 51 MODIS pixels; 162 km 2 ) moving window linear regression. Blue pixels indicated areas where either elevation or forest cover positively accounted for geographic variation (at spatial R 2 > 0.5), whereas red pixels indicated areas where either elevation or forest cover negatively accounted for a high fraction (> 50%) of the geographic variation in the association between annual SOST and TIN.

Vegetation Growth Stress Index
Regions where the MODIS phenology stress index was pronounced (annual TIN:MAXN ratio was consistently lowest at less than 35:1) were identified in the northern Brooks Range Mountains, on the Arctic Coastal Plain, and from the central Yukon-Kuskokwim Highlands southward to the Bristol Bay drainages and along the Alaska Peninsula (Figure 3a). The locations where this stress index was less pronounced (annual TIN:MAXN ratio was relatively high at greater than or equal to where forest cover (respectively) accounted for greater than 50% of the geographic variation in the mean correlation between SOST and TIN over the years 2000 to 2018, based on a regional-scale (51 × 51 MODIS pixels; 162 km 2 ) moving window linear regression. Blue pixels indicated areas where either elevation or forest cover positively accounted for geographic variation (at spatial R 2 > 0.5), whereas red pixels indicated areas where either elevation or forest cover negatively accounted for a high fraction (> 50%) of the geographic variation in the association between annual SOST and TIN.

Vegetation Growth Stress Index
Regions where the MODIS phenology stress index was pronounced (annual TIN:MAXN ratio was consistently lowest at less than 35:1) were identified in the northern Brooks Range Mountains, on the Arctic Coastal Plain, and from the central Yukon-Kuskokwim Highlands southward to the Bristol Bay drainages and along the Alaska Peninsula (Figure 3a). The locations where this stress index was less pronounced (annual TIN:MAXN ratio was relatively high at greater than or equal to 70:1) were most

Regional Case Study Locations
Surface air temperatures in northern Alaska have increased rapidly during the past several decades, according to recent field studies in the Colville and Chandler River basins of the Arctic Coastal Plain [45,46], and have resulted in site measurements of increased of plant productivity and expansion of arctic shrub cover. Consequently, areas of the Colville River around Kikak Creek (69.76°N, −151.52°W) extending south into the Lower Chandler River basin were examined more closely here as a case study of tundra greening trends (increasing TIN) in this region (Figure 4a

Regional Case Study Locations
Surface air temperatures in northern Alaska have increased rapidly during the past several decades, according to recent field studies in the Colville and Chandler River basins of the Arctic Coastal Plain [45,46], and have resulted in site measurements of increased of plant productivity and expansion of arctic shrub cover. Consequently, areas of the Colville River around Kikak Creek   Vegetation greening trends in the boreal region of interior Alaska were also detected in unburned (since the year 2000) mixed woodland and herbaceous wetland covered landscapes of the Lower Tanana River basin (Figure 4b). Both TIN and MAXN have shown steady increases over the past two decades at the representative location shown in Figure 5b, and monthly air temperatures recorded at the closest weather station have showed steady increases in the winter minima, from around −23 • C prior to 2012 to winter minima commonly above −20 • C since 2013. As in the Colville River example, precipitation totals have been trending at or above the long-term average since 2013.  Figure 5b, and monthly air temperatures recorded at the closest weather station have showed steady increases in the winter minima, from around −23 °C prior to 2012 to winter minima commonly above −20 °C since 2013. As in the Colville River example, precipitation totals have been trending at or above the long-term average since 2013.
Conversely, significant vegetation browning trends in the Yukon Flats sub-region outside of any recently burned areas were detected in boreal forest stands of the Hodzana River drainage ( Figure  4c). The TIN at this representative browning location declined rapidly starting in 2009 and has not recovered to early 2000s TIN levels, despite periodic increases in MAXN (Figure 5c). Several recent years (2015 and 2017) of extremely early SOST dates have occurred in this area, which were accompanied by increases in winter minima air temperatures, from lower than −25 °C prior to 2012 to winter minima commonly warmer than −20 °C since 2013. There have been repeated extremely dry seasons over the years 2010 to 2018 in the Yukon Flats sub-region as well, which may have contributed to declining surface moisture conditions over the past decade.
Because of the extensive NDVI phenology trends toward significantly lower TIN levels and an earlier SOST across the largest watersheds of Bristol Bay (Figure 1), the Nushagak River basin was examined further as a vegetation browning case study. Along with the Kvichak River, the Nushagak River is one of the most productive sockeye salmon (Oncorhynchus nerka) systems in the world and supports the third largest king salmon (Oncorhynchus tshawytscha) run in the world [47]. These drainages are covered by mixed forests, tundra, lakes, and marshlands. White spruce (Picea glauca) and mixed spruce-birch (Betula papyrifera) forests, as well as muskeg and willow-alder (Salix and Alnus spp.) thickets, grow up to 300 m in elevation. At lower elevations, wet tundra or marshland is common (as shown in Figure 4d), and a large tidal marsh exists at the mouth of the Nushagak River. Conversely, significant vegetation browning trends in the Yukon Flats sub-region outside of any recently burned areas were detected in boreal forest stands of the Hodzana River drainage (Figure 4c). The TIN at this representative browning location declined rapidly starting in 2009 and has not recovered to early 2000s TIN levels, despite periodic increases in MAXN (Figure 5c). Several recent years (2015 and 2017) of extremely early SOST dates have occurred in this area, which were accompanied by increases in winter minima air temperatures, from lower than −25 • C prior to 2012 to winter minima commonly warmer than −20 • C since 2013. There have been repeated extremely dry seasons over the years 2010 to 2018 in the Yukon Flats sub-region as well, which may have contributed to declining surface moisture conditions over the past decade.
Because of the extensive NDVI phenology trends toward significantly lower TIN levels and an earlier SOST across the largest watersheds of Bristol Bay (Figure 1), the Nushagak River basin was examined further as a vegetation browning case study. Along with the Kvichak River, the Nushagak River is one of the most productive sockeye salmon (Oncorhynchus nerka) systems in the world and supports the third largest king salmon (Oncorhynchus tshawytscha) run in the world [47]. These drainages are covered by mixed forests, tundra, lakes, and marshlands. White spruce (Picea glauca) and mixed spruce-birch (Betula papyrifera) forests, as well as muskeg and willow-alder (Salix and Alnus spp.) thickets, grow up to 300 m in elevation. At lower elevations, wet tundra or marshland is common (as shown in Figure 4d), and a large tidal marsh exists at the mouth of the Nushagak River.
Starting in 2014, significant vegetation browning trends (in both TIN and MAXN) were detected at a representative mixed tundra-wetland location in the Kvichak sub-basin of the Nushagak River drainage, outside of any recently burned areas (Figure 5d). After 2011, monthly air temperatures recorded at the closest weather station showed an abrupt increase in winter minima temperatures of up to +10 • C and in summer maxima of around +4 • C. Monthly precipitation totals in the region have been trending slightly above the long-term average since 2014.

Discussion
Our results from new time-series analysis of nearly 20 years of MODIS phenology metrics can be best interpreted in light of published field studies of vegetation change across all of Alaska. For instance, patterns of localized SOST-TIN correlations revealed in this study were shown to be broadly consistent with the findings of Raynolds et al. [48], who reported that lowland locations of the Brooks Range Mountains showed significant positive TIN trends, despite a delayed start to the yearly green-up period in some years, while more upland locations were detected with a significant negative TIN trend. These relationships were supported by observations that increasing green vegetation cover in this region of Alaska's arctic tundra was due mainly to changes in hydrology and nutrient availability, namely as a result of thawing of ice-rich permafrost in the greening areas. The causes for decreases in seasonal NDVI in small patches have been linked to localized disturbance events, such as erosion and thermokarst (Raynolds et al., [48]).
Warming and snow drift manipulation experiments at Toolik Lake (Alaska) in the lowland (900 m elevation) Kuparuk River basin of the northern Brooks Range foothills have shown that growth of many arctic plants has increased with both temperature warming and snow cover [49]. This relationship may be the result of deeper snow packs that can provide additional soil moisture during the warm summer months and the effect of later snow melt, pushing rapid plant growth into warmer periods. At the same study sites, Wahren et al. [50] found that snow additions had an important effect on microclimate by insulating vegetation from winter wind and temperature extremes, modifying winter soil temperatures and increasing spring run-off. The MODIS NDVI phenology-climate associations taken together from our current study results likewise imply that tundra greening in northern Alaska has been enhanced by longer, warmer growing seasons and sustained surface moisture conditions over the past decade (2008 to 2018). Gonsamo et al. [51] examined several satellite datasets (MODIS NDVI included) to study the influence of spring plant growth on summer soil moisture content at a northern Canada tundra-taiga interface. Their results suggested that higher vegetation growth in spring leads to deceased summer soil moisture and greenness. It is not yet clear whether the observations of NDVI variations in this more easterly part of North America are indicative of those in the tundra ecoregions of Alaska as well.
Across Yukon River sub-basins of interior Alaska, Potter [52] has reported significant vegetation greening trends outside of recent burned areas from MODIS NDVI monthly time-series analysis, most notably in the Yukon Flats National Wildlife Refuge and in the Susulanta Hills sub-basin of the Nowitna River drainages. These NDVI slope lines had significant (p < 0.05) deseasonalized greening trends over the period 2000 to 2018. Greening land cover for both of these sub-regions was predominantly shrubland with scattered coverage of wooded wetlands. It was hypothesized in this study that temperature warming has generally enhanced the rates of shrubland and young woodland vegetation growth across most of interior Alaska over the past 20 years. In support of this hypothesis, evidence of shrub (alder, birch, and willow) expansion over the last half century has been documented through repeat photography of the past 50 years in Alaska (Tape et al. [10]).
The results from the remote sensing studies by Potter [52,53] for large areas burned by wildfires across all of Alaska clearly showed that the NDVI will commonly recover to pre-burn levels in less than 10 years, based on trends for deseasonalized greening from MODIS 250 m resolution data sets. Findings from these studies of greenness recovery in Alaska following high-severity burns support the supposition that vegetation cover density in these ecosystems, albeit possibly with a higher fraction of annual herbaceous cover than was present before the fire, will make a nearly complete recovery often within five years post-fire. Later (SOST) seasonal greening in severely burned areas could be due to a shift to sparser coverage of evergreen forest vegetation types in the boreal ecoregion, but that is not yet widely verified with recent field observations. Putting our notable MODIS vegetation browning results reported for the first time for the Southern Alaska Mountains region into a broader eco-geographic field context, the region around Bristol Bay has experienced warming winters and more rain than snow over the past years, leading to rising water tables [54]. In watersheds like the Nushagak River basin that drain into Bristol Bay, erosion is accelerating and occurs primarily through "slump blocking," in which unconsolidated sediments from windblown glacial deposits are carried by runoff cutting into riverbanks [52]. More extreme freeze-thaw cycles and accelerated surface runoff can promote landslides, which further erode the soils and sediments.
Spruce beetle outbreaks have become widespread across southwestern Alaska and Cook Inlet forests, and have been linked to the development of warmer, longer summers, compared to historic average conditions [55,56]. Warmer winters can increase the survival rate of spruce beetle populations, particularly in ecosystems located above the snow line [57]. These insect outbreaks may be an associated indicator of increasing stress on many diverse vegetation communities in southwestern Alaska and the Cook Inlet region, especially during the warmest periods of the year. In the lowlands of the western Kenai Peninsula on the Cook Inlet, Berg et al. [58] have reported accelerating invasion of woody vegetation into wetlands. Historical aerial photography for 11 wetland sites showed that herbaceous areas have contracted significantly. Peatlands in this ecoregion have become dry enough to host black spruce only since the 1850s, and abundant shrubs only since the 1970s.
In studies at the western end of the Chugach-St. Elias Mountains ecoregion near Anchorage, Dial et al. [59] reported that tall shrub cover has increased substantially above 600 m elevation to about 750 m, mostly due to the recruitment of tall willows that were absent in 1972. These authors further postulated that rapid upward advance of woody vegetation may be contributing to regional declines in Dall's sheep (Ovis dalli) populations. Taller shrubs may have shifted upward into formerly wet meadows, reducing abundance of food plants for sheep.
For future remote sensing studies of greening trends, it is worth noting that, in a review of high-latitude NDVI studies, Myers-Smith et al. [60] have cautioned that trends in NDVI data produced from different satellite sensors or using different composting methods do not always correspond at a given location. Therefore, it remains challenging in some regions of the Arctic to distinguish between true ecological change and differences due to satellite sensor and platform-related issues when interpreting localized greening or browning trends. More to this point for arctic Alaska, Pattison et al. [61] have postulated that Advanced Very-High-Resolution Radiometer (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. We have used MODIS phenology metrics instead for this study, for which similar biases to those found in AVHRR data sets have yet to be reported widely.
There has been a recognition of uncertainties in estimation of the SOST index from satellite image data in arctic and boreal ecoregions of the world as well. These sources of uncertainty are mainly due to the presence of snow that can significantly affect the NDVI signal (Sekhon et al. [62]) and to difficulties in capturing seasonal variation in canopy greenness of evergreen vegetation cover (Jin et al. [63]). Using an improved physically-based plant phenology index (PPI) with a nearly linear relationship to green leaf area index (LAI), it was reported by Karkauskaite et al. [64] that PPI SOST preceded NDVI SOST by several days in most high-latitude shrublands and herbaceous cover, where the strongest trends in earlier SOST for Alaska were detected in the MODIS data we have presented in this study (Figure 1b). For boreal (needleleaf) forest cover areas, Karkauskaite et al. [64] reported that PPI SOST was generally detected later than NDVI SOST by 7-16 days, suggesting that the trends toward later SOST shown in our Figure 1b across Interior Alaska boreal cover may be underestimating this change in the onset of the growing season over the past two decades.

Conclusions
The evidence presented in this study supported our hypothesis that recent vegetation responses to inter-annual variability in springtime snowmelt dynamics and growing season temperatures in Alaska are, in certain ways, dependent on regional land cover types. Strong trends in SOST and TIN were associated with localized climate variations, such as rising winter temperature minima, in both northern and southern tundra regions of Alaska, but not as strongly related to merely elevation gradients or to forest versus tundra cover types. The findings in this study will help guide future research plans aimed at explaining localized changes in land cover and vegetation community responses to changing climate conditions in Alaska. At the statewide scale, regions were identified with extensive trends in earlier seasonal vegetation greening, such as the northeastern Brooks Range Mountains and in the southern coastal areas of Alaska, and with increasing growing season integrated production, such as across the western Arctic Coastal Plain and within the perimeters of most large wildfires of the Interior boreal region that burned since the year 2000. In other regions, growing season integrated production has decreased notably, particularly in watersheds of Bristol Bay and in the Cook Inlet lowlands of southwestern Alaska, in most of the same areas where an earlier-trending start of the growing season was also detected. If changes in vegetation cover types (from forest to shrubland or from wetland to grassland) have been occurring in the areas where the strongest trends in MODIS phenology metrics coincide with air temperature warming patterns, then imagery from higher-resolution remote sensing can be used to document these cover changes. There may be numerous cryptic impacts on the arctic and boreal landscapes of Alaska indicated by abrupt decreases in NDVI phenology metrics, such as insect damage, drought stress, thawing permafrost, erosion, and thermokarst events. The trend maps provided from the time-series analysis of this study can narrow the vast expanses of the Alaskan wilderness to a manageable set of local study locations where phenology has already recently shifted most notably.
Author Contributions: C.P. designed the study. C.P. and O.A. processed the data, analyzed the results, prepared display items, and wrote the manuscript. All authors have read and agreed to the published version of the manuscript.
Funding: This research received no external funding