Increased wetness confounds Landsat-derived NDVI trends in the central Alaska North Slope region, 1985–2011

Satellite data from the circumpolar Arctic have shown increases in vegetation indices correlated to warming air temperatures (e.g. Bhatt et al 2013 Remote Sensing 5 4229–54). However, more information is needed at finer scales to relate the satellite trends to vegetation changes on the ground. We examined changes using Landsat TM and ETM+ data between 1985 and 2011 in the central Alaska North Slope region, where the vegetation and landscapes are relatively well-known and mapped. We calculated trends in the normalized difference vegetation index (NDVI) and tasseled-cap transformation indices, and related them to high-resolution aerial photographs, ground studies, and vegetation maps. Significant, mostly negative, changes in NDVI occurred in 7.3% of the area, with greater change in aquatic and barren types. Large reflectance changes due to erosion, deposition and lake drainage were evident. Oil industry-related changes such as construction of artificial islands, roads, and gravel pads were also easily identified. Regional trends showed decreases in NDVI for most vegetation types, but increases in tasseled-cap greenness (56% of study area, greatest for vegetation types with high shrub cover) and tasseled-cap wetness (11% of area), consistent with documented degradation of polygon ice wedges, indicating that increasing cover of water may be masking increases in vegetation when summarized using the water-sensitive NDVI.


Introduction
The Arctic is the focus of much attention due to the fact that it is warming at twice the rate of lower latitudes (Overland et al 2014). Remote sensing has proved valuable for monitoring the Arctic and recording changes to land cover and vegetation of these remote landscapes. The normalized difference vegetation index (NDVI) calculated from satellite data has been especially effective in the treeless arctic, showing a close relationship to biomass for large areas (Walker et al 2003. Although the relationship between arctic biomass and NDVI is robust at large spatial scales, plot-level studies have shown variation depending on the vegetation type and the time of year sampled (Hope et al 1993, Riedel et al 2005. Variation also occurs when comparing the relationship of biomass to NDVI from sensors with different spatial resolutions, such as ground-sampled NDVI and satellite data (Laidler et al 2008). Thus, it is not surprising that trends have been found to vary between satellite sensors with different resolutions (Guay et al 2014), raising questions about what is actually happening to the vegetation and landscapes at finer scales.
For the North Slope of Alaska, the longest continuous NDVI time series is the 12.4 km pixel resolution AVHRR GIMMS NDVI3g data (Global Inventory Modeling and Mapping Studies, 3rd generation series from the Advanced Very High Resolution Radiometers onboard a series of NOAA satellites). Trends from these data showed continued increases from 1982 to 2011, though the rate of increase was lower towards the coast and slower in the last three years of that record compared to earlier years (Bhatt et al 2013, NOAA Arctic Report Card 2015. The slowing of the increase in NDVI occurred despite continued warming trends (Overland et al 2014). The proposed mechanism for the greening-the response of cold-limited vegetation and especially shrubs to warmer temperatures, has been documented by warming experiments (Elmendorf et al 2012a) and long-term plot data (Elmendorf et al 2012b). However, spatial heterogeneity in the response of tundra plants to warming has been noted (e.g. Myers-Smith et al 2015), with responses differing depending on summer temperatures, soil moisture conditions, and plant growth forms (Elmendorf et al 2012a). Finer-spatial resolution analyses using Landsat data for a 823 km 2 area of the foothills of the Brooks Range found that significant increases in NDVI were occurring over a relatively small proportion of the landscape (5%), when examined at Landsat 30 m pixel resolution . This landscape-scale heterogeneity in NDVI trends was also found in a Landsat time-series analysis from northern Canada (Fraser et al 2011).
Landsat TM data not only offer greater spatial resolution, but also greater spectral coverage than AVHRR data, with several more bands available for analysis. The tasseled-cap transformation is a way of incorporating the variation in all the Landsat bands into a few indices that can be interpreted based on surface attributes (Kauth and Thomas 1976). The transformation produces three major indices that account for successively less of the variation, tasseled-cap index 1, 2, and 3, related to brightness, greenness and wetness (Crist 1985, Huang et al 2002. Looking at Landsat reflectance trends with tasseled-cap indices, researchers were able to identify different types of vegetation change in northern Canada, including wildfire, fluvial and thermokarst processes (Fraser et al 2012, Olthof and Fraser 2014, Fraser et al 2014b. The goal of this research was to better understand trends in vegetation and plant biomass on the North Slope of Alaska, using Landsat data to investigate landscape-scale patterns. This coastal region has continuous permafrost, with an average ice volume of 77% (Kanevskiy et al 2013). Ice wedges between polygons are very common, with a volumetric content of about 11% for the region. Warming permafrost and warm summer temperatures have resulted in deeper annual thaw and melting of the tops of ice wedges (Jorgenson et al 2006. Figure 1 shows an example of ice-wedge degradation from the study area. A recent paper showed that this type of change is common throughout the Arctic in ice-wedge polygon landscapes (Liljedahl et al 2016). The study area has also been affected by the exploration and development of oil and gas resources. Changes in vegetation in this area are important to wildlife habitat, carbon and energy balances, industrial infrastructure, and arctic communities.

Methods
The study area is a 6075 km 2 portion of the Coastal Plain of northern Alaska, bordered by the Colville River floodplain on the west, the Shaviovik River on the east, and by the Beaufort Sea to the north (approx. 70 to 70.5°N and 147 to 151.5°W) (figure 2). The study area encompasses most of the oil development in arctic Alaska, including the original Prudhoe Bay Oilfield and adjacent oilfields (Orians et al 2003). This coastal region has a steep summer temperature gradient. Mean July temperatures at the coast are 4°C-5°C, exceed 7°C a few kilometers inland, and reach 10°C-12°C 100 km inland (Haugen 1982). Most of the region has nonacidic tundra vegetation that grows on calcareous loess deposited downwind of the major rivers (Walker and Everett 1991). The western and eastern parts of the study area have patches of acidic tussock tundra that offer a contrast to the nonacidic areas. It is mostly wetlands, with little topographic relief, ranging from 0 m a.s.l. in the north to 130 m elevation in the southeast, where the foothills of the Brooks Range begin. Mean annual temperatures at the Deadhorse airport average −7.7°C, with 108 mm of annual precipitation, mostly in the form of snow (Prudhoe Bay, WRCC 1986-1999).
Landsat scenes were downloaded from the USGS EROS GLOVIS website. We selected scenes from Path/Row 74/10 and 74/11 with 30 m pixel resolution from Landsat 4 or later (TM, ETM+), during the growing season (23 June-27 August), that were relatively cloud free over the study area (table 1). Landsat 2 and 3 MSS data were excluded because of the poor quality of the spectral data and the coarser pixel  Scenes were not available for the 1990s due to data relay problems and limited tasking during that decade (Goward et al 2006). Possible bias due to seasonal phenology was tested using AVHRR NDVI3g data to compare the bi-weekly value including the Landsat scene dates for the study area to the study area mean maximum NDVI for the year. AVHRR NDVI3g corresponded to >95% of maximum annual NDVI values (one outlier, 23 June 2008, 81%). The 2008 data were included because this was one of the few cloud-free scenes, and scenes were available from both the year before and the year after, minimizing detrimental effects on the trend analysis.
The two path/rows were mosaicked and clipped to the study extent. Gaps due to the scan line corrector malfunction were masked for the 2005 scene. Clouds were manually masked for all years except those that were cloud-free: 2001, 2002, 2005 and 2008. Top-ofatmosphere reflectance was calculated using specification from the scenes' metadata and (Chander et al 2009).
The NDVI was calculated as:
/ Tasseled-cap indices were calculated using ENVI version 5.3 software (Exelis Visual Information Solutions, Boulder, Colorado). The tasseled-cap transformation creates three approximately orthogonal indices, based on the structure of the data. It is an invariant transformation which can be applied to any scene from the same sensor (Crist 1985). Weightings used to create each index for TM and ETM band data are listed in table S1 (Crist 1985, Huang et al 2002. Tasseled-cap index 1 (TC1), 'brightness', is related to the total reflectance in all bands. Index 2 (TC2), 'greenness', is related to the amount of vegetation on the ground. Index 3 (TC3), 'wetness', is related to both surface water and soil moisture (Crist 1985).
Trends in NDVI and tasseled-cap indices were calculated using TheilSen robust regression and tested for significance using the 'mblm' package in 'R' statistical program (R 3.2.2, R Development Core Team). All unmasked pixels in each year were evaluated. Pixel locations for which there were fewer than 5 scenes with data were not included in the trend analysis. This only occurred for 21.7 km 2 in the northwest corner of the study area, which is ocean ( figure 3).
A Landsat-scale map of the vegetation of the North Slope produced for the North Slope Science Initiative (NSSI) (Ducks Unlimited Inc. 2013) was used to categorize the results. The analysis excluded 204.4 km 2 of river water including the hydrologically connected lakes of the Colville Delta and 730.3 km 2 of lake water, as mapped by Jorgenson and Heiner (2003).
The 17 NSSI vegetation types were grouped into six types that matched previous mapping of the Prudhoe Bay region (Walker et al 1980), eliminating categories with few pixels in the study area (figure 2, table S2). Barren or partially vegetated (BA) included coastal mud flats that were either barren or partially colonized by saline graminoids, as well as barren or partially vegetated riverine gravel deposits and sand dunes. Aquatic graminoid tundra (AG) included lake margins and very wet areas with standing water. Wet sedge, moss tundra (WS) covered most of the inter-lake areas of the thaw lake plain. Moist sedge, prostrate dwarfshrub, moss tundra (MS) was the most common vegetation type in the study area, occurring on hills and better-drained areas between thaw lakes. Moist tussock sedge, erect dwarf-shrub, moss tundra (TS) occurred on hillsides along the southern edge of the study area, and west of the Colville River on sandy substrates. Low to tall riparian shrubland (RS) occurred along rivers, becoming more common with distance from the coast, with the largest patches occurring along the Colville River.
NDVI and tasseled-cap index trends were examined at a finer scale by using detailed mapping in the Prudhoe Bay Oilfield area (figure 4, S1, S2). The geobotanical mapping of this area includes pre- The reason for looking at these small areas in more detail was to make use of the wealth of mapping and ground data available to help interpret the index trend analyses in terms of known changes in land cover and infrastructure.
In order to relate the satellite index results to a quantifiable vegetation characteristic, we compared Landsat NDVI and tasseled-cap greenness to aboveground plant biomass. Biomass was sampled in 2006 along the Dalton Highway in common community types (Raynolds et al 2008). These community types occur in repeating combinations based on patternedground features such as polygon centers, rims, trough and frost boil centers and rims. Biomass for specific Landsat pixels was calculated based on the proportional mapping of community types within 10×10 m grids (Raynolds et al 2008). NDVI values were extracted from the 24 August 2007 Landsat scene, the scene closest in time to the 2006 biomass sampling date.

Results
The limited amount of cloud-free summer Landsat data available due to the short growing season, frequent cloud cover, and lack of archived Landsat imagery over northern Alaska in the 1990s resulted in patchy trend values (figure 3). We compensated for this by presenting only significant trend data (p<0.05) and by summarizing over large areas to deemphasize local variability (Chen et al 2013). We focused on comparisons of different vegetation types, landform types and disturbance types, looking at both NDVI and the tasseled cap indices.
There was no significant trend in the time series of the average Landsat NDVI value for the whole study area. However, the per-pixel trend analysis using robust regression showed a slightly negative average trend of NDVI for the land area (−0.011±0.041 NDVI units/decade) ( figure 3). The per-pixel average significant NDVI trend (p<0.05) was negative overall (−0.028±0.074 NDVI units/decade) (figure 5), with values normally distributed around the mean. Only 7.3% of the study area had significant trends in NDVI, mostly along rivers and streams, lakes, and areas of industrial development (figure 5 zoom, table S3).
The two most extensive vegetation types in the study area were moist sedge (MS) and wet sedge (WS) (together over 75% of the land area, table 2). NDVI trend analysis by land-cover type showed two vegetation types with positive trends, MS and tussock tundra (TS) ( figure 6). However, only a small percentage of these units had significant increases (4% and 3%, respectively). Most vegetation types had less than 10% of their area with significant changes. AG showed the strongest decrease in NDVI, with 23% of its area having significant trends. Barren areas (BA) also had a relatively large proportion with significant trends (28%).
Trends in TC1 Brightness were generally negative, averaging −0.66±1.81 units yr -1 for the land portion of the study area. Significant trends occurred on only 2.7% of the land, with an average significant trend of −4.50±4.05 units yr -1 from 1985 to 2011 (table 2). Trends in TC2 Greenness were positive, averaging 0.59±0.77 units yr -1 . Significant trends occurred on 56.4% of the study area, with an average significant trend of 0.60±0.53 units yr -1 . Trends in TC3 Wetness were also positive, averaging 0.54±1.06 units yr -1 , with significant trends occurring on 11.4% of the area, with an average trend of 0.93±1.15 units yr -1 .
Significant decreasing TC1 brightness trends occurred more commonly on BA than other cover types, but still occurred on less than 10% of this cover type (table 2). The largest decreases were in TS. Increases in TC2 greenness occurred on most of the study area, but on only about one third of the BA and AG areas. Increases were greatest for RS and TS. Significant increases in TC3 wetness were most common in areas with AG and RS vegetation. Increases were least for WS and greatest for TS.
A graph displaying the three tasseled-cap indices along separate axes showed how distinctive the distribution of trends was for each vegetation type within the three-dimensional space (figure 7). For example, RS stood out as having both high TC2-greenness and high TC3-wetness trends. In contrast, MS showed strong decreasing TC1-brightness, but low positive trends in both TC2-greenness and TC3wetness.
NDVI trends within the three detailed map areas of the central Prudhoe Bay Oilfield were negative, similar to the whole study area (table S4). The only positive significant NDVI trend for the mapped surface form types was for pingos (table S5). TC1-brightness decreases were greatest for pingos and floodplain areas and for the forb-barren vegetation that occurs on floodplains. TC2-greenness increases were greatest for floodplains and thermokarst areas, and least for actively eroding banks. TC3-wetness increases were greatest for landforms transitioning to water (floodplains and actively eroding banks) and for vegetation types of these areas, and also aquatic marsh on the edges of eroding lakes. Most landforms (residual surfaces [19% of study area], drained thaw-lake basins [19%], and lakes [16%]) showed less change than the more active riparian landforms (rivers, streams, and active floodplains), with the exception of residual surfaces, which showed a large decrease in TC1-brightness (table S6).
Tasseled-cap indices were especially helpful in identifying changes in land cover in sparsely-vegetated or non-vegetated areas, where NDVI values were very  . In zoomed-in portion, '1' shows Alpine Oilfield development with decreased NDVI, '2' shows revegetation of drained lake with increased NDVI, '3' shows riparian erosion (red-orange decreases NDVI) and vegetation colonization and succession (green increased NDVI). (b) Graph of mean NDVI value for each year for the entire study area (excluding water). Table 2. Land-cover types within Alaska North Slope oilfield study area, area (and percent of study land surface), significant tasseled-cap trends (and percent of each type with significant trends), (units yr -1 , p<0.05). low. A composite image of the three indices (TC1=red, TC2=green, TC3=blue) for a coastal area in the northeast Colville River delta showed dark blue areas along the coastal mud flats due to either ocean transgression or changing delta deposition/erosion patterns, a new artificial gravel island as yelloworange, and shifting barrier island deposits as a mix of dark blue erosion and yellow-orange deposition (figure 8). On land, drained lakes stood out sharply, making it easy to identify nine large lakes that had drained over the study period, and nine others with orange shores indicating that they had partially drained. In contrast, most of the lakes in the study area expanded slightly, as indicated by shores showing decreased NDVI and dark blue tones in the tasseledcap composite image. A close look at the heavily developed eastern area of the Prudhoe Bay Oilfield (Map C) showed that most of the significant trends in tasseled-cap indices were related to industrial development, with changes on and adjacent to gravel pads, roads and airstrips ( figure 9, table S7). Similar trends occurred on Map A and Map B, which can be found in the supplemental material (figures S4, S5). Only TC2 Greenness showed significant trends in areas away from infrastructure.
Biomass and NDVI showed a strong linear relationship, with R 2 =0.51 and p=0.0088 (figure 10). The NDVI values were very closely related to TC2 greenness (R 2 =0.91, p<0.0001). This is not surprising as the TC greenness formula contrasts the same Landsat bands (Bands 3 and 4), though it applies no normalization (see table S1 for tasseled-cap transformation weightings). The relationship between TC greenness and biomass showed a somewhat stronger linear relationship, with R 2 of 0.65 and p=0.0016 (figure 10). Using the NDVI/biomass relationship, a change of −0.001 NDVI units/year for the plots sampled would be approximately equivalent to a decrease of 3.7 g m −2 yr −1 .

Regional trends
This study found no significant overall trend in Landsat NDVI for the study area from 1985 to 2011. This is likely due to the scarcity of the data (10 scenes over 26 years for an index that is known to fluctuate from year to year). The trend calculated using the GIMMS NDVI3g bi-weekly data that most closely matched the Landsat dates was also non-significant, whereas when all years were included there was a trend. GIMMS NDVI3g maximum annual NDVI trend for the study area for 1982-2010 was a very slight increase: 0.00276 NDVI units/decade, or a 0.51% increase/decade. In general, Landsat trends in this area have shown less increase in NDVI than the GIMMS NDVI3g (northern Brooks Range Raynolds et al 2013, northern Canada andAlaska Ju andMasek 2016). The contrast may be due to the scarcity of summer data from these cloudy areas, or due to the particular years for which Landsat data were acquired, and may also be partly due to difference between the AVHRR and Landsat sensors and the intercalibrations required to create a time series (Guay et al 2014).
The greatest changes in the Landsat trends occurred as expected, in riparian areas due to vegetation succession and erosion, with 28% of the barren or partially vegetated areas (BA) showing significant change (both positive and negative). Vegetation types occurring in the southern portion of the study area had positive trends in NDVI. This confirms other research on Alaska's North Slope, showing higher increases in NDVI in the foothills of the Brooks Range, and lesser increases on the Coastal Plain, with the least change near the coast (Stow et al 2007, Bhatt et al 2013. This result is consistent with the hypothesis that shrubs are important for greening, as shrub cover decreases with latitude and proximity to the coast (Walker 1987, Kade et al 2005. It is also consistent with the hypothesis that the most rapid changes in vegetation will take place in the southern tundra, where boreal species are already established . The other trend evident across the region was decreasing NDVI around many lake margins and decreases in the NDVI of AG vegetation commonly found around lakes. This is likely due to erosion of lake shores from the slow, continual process that has been documented by remote sensing studies in the Prudhoe Bay Oilfield Trends in tasseled-cap indices emphasized the different types of changes occurring in different vegetation types. Brightness decreases were greatest and NDVI trends were positive in vegetation types with the greatest amount of non-vascular plant cover. Mosses and lichens are often less green than vascular plants, so a decrease in their cover relative to sedges and shrubs could account for these changes. Moist tussock sedge, erect dwarf-shrub, moss tundra (TS) commonly includes an understory of mosses and reindeer lichen (Cladina spp.) (Walker et al 1994). Erect dwarf-shrub, moss tundra (MS) commonly has non-sorted circles (frost boils) covered with crustose and other types of lichens (Ducks Unlimited Inc. 2013). Lichen and moss decreases have been noted as a consistent result of warming (Elmendorf et al 2012a), and have been associated with increases in vascular plant cover (Cornelissen et al 2001), particularly shrubs (Fraser et al 2014a). These changes in plant community composition would be expected to decrease TC1brightness and increase NDVI and TC2-greenness, as was found for TS and MS types.
TC2 greenness showed general increases, but the largest increases were for vegetation types with high shrub cover-riparian shrublands (RS) and TS which has both birch and willow shrub cover. The lowest increases in TC2 greenness were seen in vegetation types with little to no shrubs: wet sedge tundra (WS) and aquatic graminoid tundra (AG). These results support the hypothesis that vegetation types with shrubs are the ones best able to respond to warming temperatures with increased vegetation growth (Myers-Smith et al 2011, Jorgenson et al 2015b. TC3 wetness increased the most in TS. This vegetation type occurs on ice-rich soils (Walker et al 1994), as evidenced by subsidence and erosion seen after fire (Jones et al 2015). Increases in TC3 Wetness in BA and RS were likely due to erosion in braided river floodplains. The smallest increases in wetness were seen in vegetation types that are already wet (WS and WG).

Landscape-scale trends
Landsat trend results were patchy across the study area due to missing data in specific years (cloud cover and scan-line corrector malfunction). However, historical and ecological ground-based explanations match the regional results presented above, and also are confirmed by the finer-resolution analysis discussed here.
Tasseled-cap indices were especially effective for identifying trends in non-vegetated areas, such as river deltas and barrier islands. The multi-spectral combination of the three tasseled-cap index trends revealed ocean transgression or loss of delta deposits along  (units/year, 1985 to 2011) for vegetation types within the Alaska North Slope oilfield study area. Origin in lower right corner has least change, far upper corner has most change. AG=aquatic graminoid tundra; BA=barren or partially vegetated; MS=moist sedge, prostrate dwarf-shrub, moss tundra; RS=low to tall riparian shrub tundra; TS=moist tussock-sedge, erect dwarf-shrub, moss tundra; WS=wet sedge, moss tundra. Vegetation types were from Ducks Unlimited Inc. (2013), reclassified into modified physiognomic categories (Walker et al 1980). much of the coast of the Colville River Delta. Barrier islands showed the expected movement of sediments westward and towards the coast (Hopkins and Hartz 1978). Both these processes were also seen in the Lena Delta in Russia (Nitze and Grosse 2016). While regional trends were all due to natural and climaterelated changes, industrial activities were clearly evident at the landscape scale. Both NDVI and tasseledcap trends were significant where construction of gravel pads and roads occurred on tundra (decreases in all indices except increased TC1 brightness).
Within the detailed maps in the central Prudhoe Bay oilfield (Maps A, B, C), results based on vegetation types were consistent with the regional trends. NDVI decreases were greatest for AG. This type occurred along lake shores, so its loss was consistent with the erosion of lake shores reported in this area . NDVI decreases were also evident in areas of industrial disturbance including excavations, flooded areas, and vehicle tracks. TC2 Greenness increases were greatest for revegetating areas, both natural succession along rivers and streams and revegetation of areas of industrial disturbance. TC3 Wetness showed the greatest increases due to eroding river gravel banks and in areas of thermokarst.
TC1 brightness decreases were greatest for residual surfaces experiencing thermokarst resulting from permafrost degradation. These areas have the oldest undisturbed soil, with the most time for ice-segregation permafrost features to develop. The thawing of polygonal ice wedges both adjacent to infrastructure and in undisturbed tundra in this area has resulted in standing water and lusher vegetation in troughs between the polygons (figure 1) . Though these changes occur on a scale smaller than a Landsat pixel, both the increased standing water and the increased vegetation within a pixel would reduce the TC1 brightness. This ice-wedge degradation process has been documented from ground studies within the study area (Jorgenson et al 2015a), and documented throughout the Arctic (Liljedahl et al 2016).

Conclusions
Landsat time-series analysis and accompanying tasseled-cap transformation analysis were useful in unraveling the various changes occurring to vegetation on the North Slope, despite limitations due to sparse data. Regional trends were due mostly to natural and climate-related changes in vegetation and ice-wedge thermokarst, though in the most heavily developed parts of the study area most changes were industryrelated.
TC greenness trends were consistently positive, while NDVI trends were mostly negative. Both indices were closely related to aboveground plant biomass. We suggest that the increase in the TC Wetness is the likely reason that regional NDVI trends were negative while regional TC greenness trends were positive. NDVI is known to be sensitive to water cover, which has a NDVI value of zero, and to be affected (directly and also indirectly through the vegetation) by soil moisture (Laidler et al 2008, Lin et al 2011). As a result, a pixel with thermokarst degradation, with increasing wetness due to standing water and increasing greenness due to lush sedge growth in the troughs, could well have a negative trend in NDVI. In this study area, biomass was more closely related to TC Greenness than to NDVI, so vegetation biomass may be increasing slightly, as indicated by the AVHRR NDVI trends for the area, even though Landsat NDVI showed a small negative overall trend.
Individual ice-wedge thermokarst pits are smaller than the pixel size of Landsat images so it is difficult to positively ascribe the cause of negative NDVI changes  solely to thermokarst. In this study, we could identify a likely cause because of the abundance of high-resolution air photos, VHR satellite imagery, and a historical GIS database that traced the landscape history of the region. A more firm conclusion regarding the cause of NDVI changes will require consistent, long-term monitoring using very-high resolution satellite imagery combined with regular ground sampling of vegetation composition and biomass. There also remains a need to examine the regional consequences of industrial development in more detail.