Changes in greening in the high Arctic: insights from a 30 year AVHRR max NDVI dataset for Svalbard

Satellite-aided studies of vegetation cover, biomass and productivity are becoming increasingly important for monitoring the effects of a changing climate on the biosphere. With their large spatial coverage and good temporal resolution, space-borne instruments are ideal to observe remote areas over extended time periods. However, long time series datasets with global coverage have in many cases too low spatial resolution for sparsely vegetated high latitude areas. This study has made use of a newly developed 30 year 1 km spatial resolution dataset from 1986 to 2015, provided by the NOAA AVHRR series of satellites, in order to calculate the annual maximum NDVI over parts of Svalbard (78°N). This parameter is indicative of vegetation productivity and has therefore enabled us to study long-term changes in greening within the Inner Fjord Zone on Svalbard. In addition, local meteorological data are available to link maximum NDVI values to the temporal behavior of the mean growing season (summer) temperature for the study area. Over the 30 year period, we find positive trends in both maximum NDVI (average increase of 29%) and mean summer temperature (59%), which were significantly positively correlated with each other. This suggests a temporal greening trend mediated by summer warming. However, as also recently reported for lower latitudes, the strength of the year-to-year correlation between maximum NDVI and mean summer temperature decreased, suggesting that the response of vegetation to summer warming has not remained the same over the entire study period.


Introduction
Time series of satellite based normalized difference vegetation index (NDVI) data have significantly improved our understanding of intra-and inter-annual variations in vegetation properties from a regional to global scale [1,2], with wide implications for our predictive understanding of ecological processes in time and space [3,4]. The NDVI index is defined as where R nir and R red are the measured reflectance in the near infrared and red channels respectively. Hence, the NDVI is represented on the scale 0 to 1, and lies typically in the interval 0.05-0.7 for increasing amounts of green vegetation [5]. This index is known to be highly correlated with the absorbed fraction of photosynthetically active radiation and therefore gross photosynthesis processes [6][7][8]. NDVI is also found to be correlated with above-ground tundra biomass [9,10] and has been adopted and applied as a proxy for leaf area index [11,12] and net primary production [13]. NDVI is thus a very useful parameter by which to measure vegetation growth for large and remote land areas.
Several different global NDVI data series exist. Since 2000 the MODIS sensor on board the Terra satellite has provided daily NDVI data at 250 m spatial resolution, with aggregated products such as the 16 day maximum NDVI data with 1 km spatial resolution MOD13A2. However, even longer time-series exists. The National Oceanic and Atmospheric Administration (NOAA) Advanced Very High Resolution Radiometer (AVHRR) series of instruments offers the longest possible datasets dating back to 1981. The most widely used AVHRR dataset is the Global Inventory Modelling and Mapping Studies (GIMMS) [14]. This dataset comprises 15 day maximum value composites at 8 km spatial resolution now covering the July 1981-December 2015 time period. However, for very sparsely vegetated areas such as those at high latitudes, this is in many cases insufficient spatial resolution.
There are very few field observation sites with long time series of growing season and plant productivity parameters in the high Arctic. On Svalbard, only a few field phenological and biomass studies have been carried out [15][16][17], and from these point observations and maps [18,19] it is difficult to detect long-term trends and the year-to-year spatial variability on the archipelago. However [20], used MODIS data for the 2000-2013 period and mapped the onset of the growing season on Svalbard, and revealed large annual variability, but no clear trend during the 14 year period as a whole. This indirectly indicates that the patterns recently observed at lower latitudes, with weakening greening trends in response to increased summer temperatures [21], also occur in the northernmost regions of the Arctic.
The aim of this study is twofold; (1) to develop the longest possible NOAA-AVHRR based NDVI time series for Svalbard; and (2) to study temporal patterns of greening trends and their link to local measurements of summer temperature. Remote sensing based mapping of the growing season in the Arctic zone is, however, a challenging task. This is due to the large local differences created by local topography, snow cover and hydrology. This variability is difficult to detect due to the limited spatial resolution of most of the current satellite sensors that have high temporal resolution. Frequent cloud cover, short growing season and insufficient response from scattered vegetation cover further increase the challenges. In particular, frequent cloud cover during the short time period from snowmelt to onset of the growing season makes this phenology parameter difficult to measure accurately from the NDVI data. In this study we have therefore focused on studying the growing season maximum NDVI value (hereafter MaxNDVI). The MaxNDVI is the highest summer NDVI value, representing peak vegetation photosynthetic activity, and serves as an indicator of tundra biomass [22,23]. This is the parameter that is least influenced by noise and cloud cover due to the relatively longer time period available during the entire growing season.

Study area
The archipelago of Svalbard is located in the Arctic part of Norway extending roughly from 76°30′ to 81°N (figure 1) and covers an area of approximately 62 000 km 2 . A large part of the archipelago is dominated by glaciers, amounting to more than 60% of the total land area. The Nordenskiöld Land peninsula, located in the western-central part of the main island Spitsbergen, is the largest ice-free area in the archipelago and is also the area with the warmest summer temperatures. The climate station in the administration center Longyearbyen (figure 1) exhibits a mean July temperature of 6.5°C (during the standard normal period . Approximately half of the ice-free area is without vegetation, and the largest areas with dense vascular plant dominated plant communities are also found in the central part [24]. Vegetation is diverse [25] and includes closed tundra, fens, mires and even high-Arctic steppes in the warmest parts. 2.2. Data processing 2.2.1. Satellite data AVHRR produces imagery at 1.1 and 4 km spatial resolution. The 4 km or global area coverage (GAC) product is derived from the 1 km product by onboard sampling taking only one line out of every three and averaging every four of five adjacent samples along the scan line. Time series such as the GIMMS dataset is obtained from the GAC data, whose actual resolution is further reduced to 8 km×8 km for the final result of the NDVI time-series.
The local area coverage (LAC) data set is recorded onboard at original resolution (1.1 km) for part of an orbit and later transmitted to earth. The high resolution picture transmission data is real-time 1.1 km downlinked data. Since the launch of MetOp-A in 2006 a fourth data type, Full resolution area coverage (FRAC 1.1 km), is available daily for the entire globe.
However, the 1 km LAC record is not continuous. Its availability depends upon prior arrangements made by NOAA, or on the proximity of a local receiving station that can capture the data directly from the satellite. Full resolution data were recorded on an onboard tape for subsequent transmission during a station overpass. The average instantaneous field-ofview of 1.4 milliradians yields a LAC ground resolution of approximately 1.1 km at the satellite nadir from the nominal orbit altitude of 833 km. A total of 13060 AVHRR LAC scenes were downloaded from the archive to cover the Svalbard snow-free season (June through September) for the period 1986-2015 (see table A1 in appendix). Prior to 1986 there is not enough data available to cover the full growing season.
To cover the entire 30 year period, data from 9 satellites were used (NOAA-9, 11,12,15,16,17,18,20,21). Data from band 1 (0.55-0.68 μm) and band 2 (0.73-1.1 μm) were used to calculate the NDVI. The individual channel digital counts were transformed to top-of-atmosphere reflectances using the models provided by [26] for NOAA-9 and NOAA-11, [27] for NOAA-12 and [28] for NOAA-15. These models provide time-varying calibration coefficients, compensated for sensor degradation. For satellites from NOAA-17 and later, the updated calibration coefficients available in the header of the data files were utilized. The channels on the various AVHRR sensors cover slightly different spectral ranges. Hence, the NDVI values calculated from different sensors will vary slightly. Therefore, the NDVI values from the various AVHRR sensors were cross-calibrated to NOAA-9 by utilizing the models specified by [29,30].
The Mt. Pinatubo eruption in 1991 injected large quantities of aerosols into the Earth's stratosphere. These aerosols, along with smoke from biomass burning and dust from soil erosion and other factors, can introduce variability in the AVHRR NDVI record. However, the simulations published by [31] indicates that the aerosols only affects observed radiance for latitudes between approximately 20°S and 40°N. Accordingly, we did not apply a specific atmospheric correction for the period affected by Mt. Pinatubo stratospheric aerosol (1991)(1992)(1993). Note, however, the low MaxNDVI recorded in 1992 (figure 2), which may be due to natural variation rather than aerosols affecting radiance per se. NDVI is also sensitive to the periodic variations in solar illumination angle and sensor view angles induced by the satellite orbits. We have not performed any correction to the NDVI time series for varying solar zenith angle (SZA), since it has been reported that at high northern latitudes the SZA component represents only a very small part of the NDVI signal [32].
The finished NDVI products are first geo-rectified onto a map projection using the coordinate grid that accompanies the data sets, using the method described in NOAA KLM User's Guide and NOAA Polar Orbiter Data User's Guide. However, for the earliest satellites (prior to year 2000), the coordinate grid provided is very imprecise and deviations of more than 20 km are not uncommon [33]. Although precise orbital parameters are available for all NOAA satellites (http:// celestrak.com/NORAD/archives/), the AVHRR data set does not contain information on the satellite's attitude. Hence, we cannot utilize the precise orbital information to post-process the AVHRR data to obtain a more precise geo-location. Instead, we have performed a manual co-registration of every scene used in the NDVI time series. Prior to this work, a subset of the best scenes (i.e. minimum clouds and free of visual artefacts and scan line errors) was manually selected to span the complete 30 year time period reducing the dataset to 760 scenes (see table A1 in appendix).
The spatial resolution of the final product was set to 1 km, which is similar to the ground resolution of AVHRR at nadir. This is higher spatial resolution compared with many previous long time studies of AVHRR data, often based on GAC data that has typically been resampled to 8 km spatial resolution [34]. The increased spatial resolution is necessary due to the extremely sparse vegetation on Svalbard, but comes at the cost of increased noise and pixel correlation of off-nadir pixels. Finally, two composite NDVI images were generated for each month, taking the maximum pixel-wise NDVI value in the period from day 1 to day 15 and from day 16 to the end of the month, respectively. Although this method does not eliminate atmospheric and surface bidirectional effects, it does act to minimize such effects without an explicit atmospheric correction being required [35]. In addition, scan-angle effects, cloud contamination, and shadow effects due to varying SZAs and SZAs are minimized [36].

NOAA-AVHRR masks on Svalbard
We have applied a mask that covers most of the central part of Svalbard to study the temporal variations in NDVI where there is most dense vegetation during the summer months. This Inner Fjord Zone (IFZ) is the area with the supposed highest summer temperatures, and includes most of the Middle Arctic Tundra Zone on Svalbard [37]. In order to construct the final mask, the mean NDVI was determined from averaging over the period 4 July-3 August using a MODIS-based dataset from 2000 to 2013 [20]. We have thresholded the mask such that only pixels with mean NDVI values above 0.2 remained in the final mask. For the trend in max NDVI and correlation computations we averaged the NDVI values within this area. We have selected the IFZ since this area has the densest vegetation cover of vascular plants and highest biomass [19,24]. This area also excludes some of the large flat areas in Svalbard with moss-tundra, since moss-tundra (with less influence of vascular plants) may exhibit high NDVI values just after snowmelt. In such cases the NDVI could be attributed to moisture rather than temperatures [20,38].

Temperature data
Daily mean temperature data were obtained for the Svalbard airport observation site (located in the middle of the IFZ) using the meteorological data portal eklima.met.no, a website which provides access to a climate database maintained by the Norwegian Meteorological Institute. We downloaded the data for the years from 1986 through to the end of August 2015. Earlier studies of Arctic tundra have found that NDVI values increase from late May and peak in early August [39], and that early season temperatures contribute most to NDVI annual variability. Although the greening of vegetation does not start before mid June [20], the snow may melt well before the greening [16]. Thus, we have chosen the period from May 16th through to August 1st over which to average the daily temperature data in order to obtain a temperature that is representative of the growing season on Svalbard.

Extracting trends in time series
Because of the potential sources of data noise described above, some outlier years with a biased estimate of max NDVI (e.g. due to long periods of cloud cover during peak season) can be expected in such a long time-series. To partly account for the presence of such outliers, we performed our analyses also on subsets of the time-series. For all possible pairs of years, we removed the corresponding temperature and maximum NDVI values from the time series such that 28 years remained. The linear Pearson correlation coefficient was calculated for each set of 28 years (as well as the full time-series) in order to identify which pair of years, when removed from the dataset, produced the best correlation between mean summer temperature and maximum NDVI for the IFZ. In order to determine the temporal trend in temperature and maximum NDVI we have subsequently performed linear regression on the subset of the time series.
Furthermore, we used the obtained linear fits to de-trend the datasets. By removing the underlying trends we can extract the inter-annual variations in mean summer temperature and maximum NDVI. To investigate in greater depth the relationship between maximum NDVI and mean summer temperature throughout the 30 year period, we have carried out a similar method as performed by [21]. In this work the authors used the de-trended time series and calculated the partial correlation coefficient using a 15 year moving window in order to remove statistically the effects of other parameters (for example cloud cover and precipitation) that may influence the apparent correlation between the two variables. By calculating the partial correlation coefficient in this way, we were also able to examine whether the observed sensitivity of maximum NDVI to mean summer temperature remained constant throughout the period of study.

Results
Maximum NDVI and mean summer temperature exhibited positive linear trends of 0.0017 yr −1 (29% increase; t=4.02, P<0.05) and 0.065°C yr −1 (59% increase; t=5. 15 Qualitatively, the temporal variations of mean summer temperature and maximum NDVI follow each other rather closely. Thus, with all years included, the correlation coefficient between temperature and maximum NDVI was 0.62 (p=10 −4 ). The outlier removal procedure suggested that elimination of 1989 and 1995 gave the highest correlation coefficient (r=0.716, p=10 −5 ). In the case of the outliers, for 1989 the temperature decreases (with respect to 1988) to the lowest on record, while NDVI increases slightly. For 1995 we observe an increase in mean summer temperature (with respect to that of 1994) to the highest on record at that time, but there is almost no change in maximum NDVI.
The correlation coefficient for the de-trended time series of maximum NDVI and mean summer temperature (figure 3) can be seen to have declined over the time interval, indicating that the sensitivity of maximum NDVI to mean summer temperature has weakened since the start of the period of study. However, the weakening in correlation over time can be most likely attributed to greater spread in values rather than a nonlinear relationship between NDVI and temperature (figure 4).

Discussion
In this study we have taken advantage of a 30 year, 1 km-resolution AVHRR dataset to study temporal patterns in the yearly values of maximum NDVI for the IFZ on Svalbard, and how these patterns relate to summer warming. Our results indicate that while both mean summer temperature and maximum NDVI have increased over the study period, the rate of increase for maximum NDVI is lower over the second half of the time-series when compared with the rate of increase for the 30 year period as a whole ( figure 2). Accordingly, the overall positive correlation between maximum NDVI and temperature decreases from year to year during the course of the study ( figure 3).
Similar increases in NDVI and, hence, greening trends have been reported in early studies of trends in seasonal NDVI at latitudes greater than 65°N [2], where the Pathfinder and GIMMS AVHRR time series products were used to show that the seasonal NDVI increased by 8%/9 yr (27%/30 yr) and 14%/8 yr (52%/30 yr) respectively. In that study, however, the time period covered by the dataset ranged only from 1981 to 1991. Likewise, our finding that the relationship between mean summer temperature and maximum NDVI is weakening is also similar to earlier published work on the GIMMS dataset, showing that the strength of the correlation between the inter-annual variability of NDVI and temperature became weaker over the period 1982-2011 [21]. We note however that those results were far more generic and applied to a much greater and diverse area of study, spanning all latitudes greater than 30°N. One explanation for the weakening in correlation, as suggested by [21], is that the response of photosynthesis to temperature may be nonlinear, but we did not find evidence for such a nonlinear relationship. Other explanations may be related to shrub expansion over grass-dominated tundra [21], which is less likely at high latitudes such as Svalbard, where only dwarf shrubs such as Salix polaris and Cassiope tetragona are found. Although the role of mosses in our maximum NDVI estimates also remains largely unknown, their rate of photosynthesis is independent of temperature and solar radiation [38,40] since they have no roots, thus meaning that they grow significantly better than vascular plants during poor summers with lowered temperatures and cloudy weather. Hence the  variation between species and functional groups in the strength of delayed effects and sensitivity to temperature per se [17] may indeed have contributed to the observed changes in trend patterns over time.
However, some authors have identified certain situations that may result in a negative feedback in vegetation growth. For example, a short and intense warming event during the winter/spring period would lead to earlier snow melt, but this may also result in exposure of vegetation to colder temperatures earlier in the year, thereby leading to damage to buds and possible reduction in greening as well as damage on dwarf shrubs as observed in Fennoscandia [41][42][43].
Other studies have also documented experimental evidence that illustrate the profound effects of anomalously warm winter spells on Svalbard. These spells can result in intense rain-on-snow (ROS) events, which may contribute to the lowered value of maximum NDVI as reported on by [44] and formation of ice layers. These can in turn affect the subsequent summer vegetation growth and maximum NDVI [45] and lead to damage of lichens [46]. Such damage has been observed at different locations in the study area and elsewhere on Svalbard during the years 2012-2015, especially in dwarf shrubs such as Cassiope tetragona. Ongoing work using high resolution satellite imagery and Unmanned Aircraft Systems (UAS) combined with ground investigations is in progress in order to examine this.
The frequency of such extreme winter warming events in parts of the Arctic has increased [47,48] and can perhaps partly explain the tendencies for 'Arctic browning' [49] and our observed weakening in correlation between maximum NDVI and mean summer temperature as also reported elsewhere [21,50,51]. Thus, the possible long-lasting effects of extreme winter warming events as outlined by previous studies cannot be ruled out as a possible contributor to the slower rate of greening exhibited in our dataset.
Alternatively, a decline in the rate of increase in maximum NDVI while mean summer temperature appears to increase may also be explained by changes in the onset and length of the growing season. AVHRR and MODIS data used to study trends in growing season parameters report that growing season length has increased over the last three decades at high latitudes due to earlier onset and delayed end of season, however with mostly less increase or even regionally a decreasing trend since approximately the millennium shift [50][51][52][53][54][55][56][57]. Moreover, the increase of the growing season productivity for herbaceous, shrub and tundra vegetation was reported to have flattened out during the last decade of the period 1982-2014, with the most pronounced decline for tundra vegetation after 2010 [51].
Another possible explanation for our observed temperature-max NDVI relationship could be related to the results of [50] who used GIMMS data to study max. NDVI, reporting positive trends for the period 1982-2011 but with weaker trends in the most recent period, as observed in our study. They also examined trends in the summer warmth index (SWI) and time integrated NDVI and noted declines in both parameters in line with changes in sea level pressure in the Arctic. These sea level pressure changes consequently resulted in large-scale circulation changes, leading to increased cloud cover and therefore lower SWI and lower max NDVI over Eurasia. Nonetheless, in the case of Svalbard, our results indicate a steady increase in summer temperatures in this region, which makes this mechanism a less likely candidate explanation for our observed trends.

Conclusion
Over the 30 year study period covered by the NOAA AVHRR NDVI data, a steady summer warming has resulted in a positive trend in maximum NDVI in high Arctic Svalbard and, hence, an overall greening trend similar to those found at lower latitudes in the Arctic. However, the rate of greening during the second half of the study period was observed to have slowed relative to the first half of the study period while the rate of increase in summer temperature remained relatively constant. This invariably resulted in a reduced year-to-year correlation between maximum NDVI and mean summer temperature over the duration of the study period. Although the mechanisms for this reduced sensitivity to summer temperatures remain unknown, our findings have great implications for our understanding of how current and future climate change affects ecosystem structure and dynamics in the high Arctic [58,59].

Acknowledgments
We gratefully acknowledge the contribution of the late Senior Researcher Inge Lauknes during the development of the processing software, Cecilie Sneberg Grøtteland and Lise Revheim Sandal for continued software development and data processing and to Stine Skrunes for the tedious work of georeferencing the oldest datatsets. This work was partly financially supported by the Research Council of Norway's projects 'Predicting effects of climate change on Svalbard reindeer population dynamics: a mechanistic approach' (POLARPROG grant 216051) and 'The effects on population dynamical key-processes of a changing climate in arctic ecosystems' (NORKLIMA grant 178561/S30), and the ArcticBiomass project (grant 227064/E10), as well as the Polish-Norwegian programme of the EEA Norway Grants (WICLAP project, ID 198571).