Changes in growing season duration and productivity of northern vegetation inferred from long-term remote sensing data

Monitoring and understanding climate-induced changes in the boreal and arctic vegetation is critical to aid in prognosticating their future. We used a 33 year (1982–2014) long record of satellite observations to robustly assess changes in metrics of growing season (onset: SOS, end: EOS and length: LOS) and seasonal total gross primary productivity. Particular attention was paid to evaluating the accuracy of these metrics by comparing them to multiple independent direct and indirect growing season and productivity measures. These comparisons reveal that the derived metrics capture the spatio-temporal variations and trends with acceptable significance level (generally p < 0.05). We find that LOS has lengthened by 2.60 d dec−1 (p < 0.05) due to an earlier onset of SOS (−1.61 d dec−1, p < 0.05) and a delayed EOS (0.67 d dec−1, p < 0.1) at the circumpolar scale over the past three decades. Relatively greater rates of changes in growing season were observed in Eurasia (EA) and in boreal regions than in North America (NA) and the arctic regions. However, this tendency of earlier SOS and delayed EOS was prominent only during the earlier part of the data record (1982–1999). During the later part (2000–2014), this tendency was reversed, i.e. delayed SOS and earlier EOS. As for seasonal total productivity, we find that 42.0% of northern vegetation shows a statistically significant (p < 0.1) greening trend over the last three decades. This greening translates to a 20.9% gain in productivity since 1982. In contrast, only 2.5% of northern vegetation shows browning, or a 1.2% loss of productivity. These trends in productivity were continuous through the period of record, unlike changes in growing season metrics. Similarly, we find relatively greater increasing rates of productivity in EA and in arctic regions than in NA and the boreal regions. These results highlight spatially and temporally varying vegetation dynamics and are reflective of biome-specific responses of northern vegetation during last three decades.


Introduction
Boreal and arctic ecosystems cover 22% of the terrestrial surface and stretch over North America (NA) and Eurasia (EA). They play a crucial role in the Earth system by regulating energy-water-carbon exchanges between the land surface and the planetary boundary layer (Chapin et al 2000). During the last half-century, these regions have experienced temperature increases of 0.3°C-1.0°C per decade higher than anywhere else on the Earth, particularly during the winter and spring seasons (Solomon 2007). A changing thermal regime and its consequences on physical, hydrological and biogeochemical conditions such as snow depth, soil moisture, disturbance etc. have already affected northern vegetation structure and function (Walther et al 2002). For example, increasing shrub cover across a broad range of hemispheric tundra area, termed as 'shrubification', has been documented (Tape et al 2006) and changing tree growth has been observed in NA boreal forests . As these changes may feedback on regional and global climate, an accurate characterization of changes during the recent past and some idea of future changes is a critical topic of research.
As a way to diagnose vegetation response to climate change, monitoring growing season duration and productivity has drawn particular attention because these are sensitive and easily measurable indicators (Richardson et al 2013). Field studies have indicated that the growing season duration for northern vegetation has significantly increased over the past decades due to both an earlier start and delayed ending (Parmesan andYohe 2003, Menzel et al 2006). This is generally thought to result in a longer carbon assimilation period due to relaxation of low temperature limits on metabolism (Nemani et al 2003), and in turn an increase in primary productivity (Xu et al 2013, Forkel et al 2016. Indeed, ground observations confirm enhanced productivity from a lengthened photosynthetically active period (Richardson et al 2010, Keenan et al 2014. Satellite observations have been employed to monitor and understand changes in growing season duration and productivity at large spatial scales. Remote sensing data reveal widespread lengthening of the growing season and an increase in gross primary productivity, also called 'greening', both of which are associated with warmer air temperatures in the high latitudes during 1980s and 1990s (Myneni et al 1997). After this period, divergent responses in productivity between boreal (decrease in productivity called 'browning') and arctic (contiguous greening) vegetation (Goetz et al 2005, Piao et al 2011, Bjerke et al 2014, and a reduced or reversed rate of regional growing season changes were also reported (Høgda et al 2013, Wang et al 2015. Furthermore, asymmetric seasonal warming (Serreze et al 2000) and a multitude of drivers greatly complicate the characterization of variations in growing season duration and productivity. This complexity justifies the need for a comprehensive examination of the magnitude and direction of changes across the northern hemispheric landscape using the longest satellite data set currently available.
The primary objectives of this study are to: (1) evaluate the reliability of long-term growing season duration and productivity metrics inferred from satellite data, (2) investigate the spatiotemporal pattern and trend of changes in growing season duration and productivity, and (3) quantify changes across continents (EA and NA), biomes (arctic and boreal) and vegetation types. To achieve these objectives, we used a satellite dataset covering the northern high latitude region (>45°N) for the period 1982-2014 (33 years long). We first define pixel-wise growing season duration and productivity metrics, then introduce independent datasets to assess the reliability of metrics inferred from satellite data. Robust statistical tests and trend analyses are used to evaluate long-term vegetation dynamics.

Materials and methods
This study is focused on vegetation in the boreal and arctic regions depicted in figure S1. We define 12 subvegetation classes and 4 vegetation groups using the Moderate Resolution Imaging Spectroradiometer (MODIS) International Geosphere-Biosphere Programme land cover (Friedl et al 2010) and Circumpolar Arctic Vegetation Map (CAVM, Walker et al 2005). Details for vegetation map can be found in SI section 1. All data sets used in this study are briefly described in SI data section and their spatial resolutions are identically harmonized into 1/12°for comparison purpose.
2.1. Determination of long-term (33 year) growing season and productivity Normalized difference vegetation index (NDVI) is a radiometric measure of the amount of photosynthetically active radiation (∼400-700 nm) absorbed by vegetation. It is calculated from contrasting reflectances at near-infrared (ρ nir ) and red (ρ red ) bands: NDVI=(ρ nir −ρ red )/(ρ nir +ρ red ) (Tucker 1979). NDVI has been widely used in studies of phenology, productivity, biomass and disturbance monitoring as it has been proven to be a good surrogate of vegetation photosynthetic activity (Pettorelli et al 2005). Here, we used the latest version of Global Inventory Modeling and Mapping Studies (GIMMS) NDVI dataset (NDVI3g) which is spanning from July 1981 to December 2014 with a native resolution of 1/12°at bimonthly time steps (Pinzon and Tucker 2014).
The growing season summed NDVI (or, GSSNDVI) has been found to be a good proxy for vegetation gross primary productivity (Goward et al 1985, Wang et al 2004. We derived long-term GSSNDVI from 1982 to 2014 through the inferred corresponding growing season metrics: onset, end and length of growing season (SOS, EOS and LOS, respectively). Two preprocessing steps were first performed to maintain distinct seasonal vegetation trajectory and minimize spurious signals (e.g., cloud and snow): (1) implementing the Savitzky-Golay filter to smooth the NDVI3g time series (Chen et al 2004, Jönsson andEklundh 2004); (2) identifying background NDVI and replacing NDVI that varied irregularly during the winter period (e.g., Beck et al 2006). After that, we linearly interpolated the dataset to a daily time step. We also use the daily freeze/thaw (FT) state of the ground to define the photosynthetically active period because vegetation may remain green during the dormant season (SI section 2.9).
Based on the daily NDVI and FT time series, we define pixel-wise photosynthetically active growing season metrics as follows (figure S2, Zhu et al 2016): (a) SOS is the day when the NDVI value is greater than 0.1 and has increased by 25% of the growing season amplitude; (b) EOS is the day when the NDVI value is greater than 0.1 and has decreased by 25% of the growing season amplitude; (c) the ground should be in thawed state; (d) LOS is the duration between SOS and EOS. Based on extracted growing season, the pixelwise GSSNDVI for each grid (p) and year (y) can be calculated by cumulating daily NDVI ( f NDVI (t)) over LOS as below . Temperature based potential SOS (PSOS), EOS (PEOS), LOS (PLOS) and growing season summed warmth index (GSSWI) were also used (see SI section 2.10) as the temporal coverage, i.e. number of years, of other evaluation datasets was limited.
The cross-comparisons were performed at both site and continental scales. For site scale evaluation, we selected 109 sample sites based on the latest Benchmark Land Multisite Analysis and Intercomparison of Products (BELMANIP-2) scheme as it provided a good sampling across biomes and land surface types (Baret et al 2006). In flux tower versus satellite data comparisons, valid flux sites and data were ascertained as follows: (i) more than 95% of the days had daily GPP data, and (ii) the mean daily quality flag was more than 0.75 (Richardson et al 2010). For continental scale comparisons, all the metrics were converted to anomalies with respect to their common period and then spatially averaged over North America (NA), Eurasia (EA) and the entire circumpolar (CP) region.

Results and discussion
3.1. Evaluation of NDVI3g based growing season and productivity metrics The NDVI3g based metrics of growing season (SOS, EOS and LOS) and seasonal total productivity (GSSNDVI) agree well with corresponding metrics derived from other evaluation datasets (table 1 and figures 1 and 2). Table 1 provides a summary of comparison for the 109 BELMANIP-2 sites distributed over the northern vegetated lands. The correspondence between NDVI3g based metrics and those from the improved MODIS NDVI (MCD43C4) is good-R 2 and RMSE of 0.96 and 5.23 days for SOS, 0.77 and 9.23 days for EOS and 0.89 and 12.37 days for LOS, respectively. Similarly, reasonable agreement is seen for SOS (R 2 =0.74) and LOS (R 2 =0.59) between NDVI3g and the MODIS phenology product (MCD12Q2). However, EOS from NDVI3g tended to be much later (bias=22.42 days). This could be due to varying data-fitting techniques (Savitzky-Golay versus piecewise logistic) and/or detection methods (amplitude threshold versus curvature; White et al 2009, Ganguly et al 2010. GSSNDVI captures the spatiotemporal patterns of productivity metrics derived from the other datasets (table 1). The GSSNDVI explains more than 80% of the spatial variation in the MODIS GPP product (MOD17A3; R 2 =0.81). Additional comparisons with MTE-GPP indicate that GSSNDVI captures both spatial (R 2 =0.85) and temporal (R=0.52) variations in gross primary productivity.
At the continental scale, the SOS metrics from NDVI3g and temperature (PSOS) agree quite well (R=0.86 in EA, R=0.80 in NA, R=0.82 in CP) ( figure 1(A)). They show a gradual transition from positive to negative anomalies, thus demonstrating advancing onset of thermal and photosynthetic growing seasons during last three decades. Unlike SOS, the NDVI3g based EOS metric does not exhibit a close correspondence with temperature-based PEOS over last three decades ( figure 1(B)). Previous studies have noted that while SOS of northern vegetation is mostly controlled by preseason temperature, EOS has multiple driving factors such as photoperiod, temperature, nutrients, etc (White et al 1997, Gill et al 2015. Nevertheless, a close association between LOS and PLOS is seen (R=0.74 in CP; figure 1(C)). We also note wellsynchronized temporal variations with MODIS LOS metrics from three different MODIS datasets (MCD12Q2, MOD13C1 and MCD43C4). Overall, NDVI3g based growing season metrics reveal good temporal agreements with those of MODIS although we observed some deviations in later common period (2012∼), particularly for SOS and LOS. This divergence has been caused by the differences in NDVI/EVI response to vegetation growth between sensors, rather than by the processing methods (see detail explanations in figure S3).
GSSNDVI at the circumpolar scale provides a reasonable representation of the long-term MTE-GPP (R=0.67) and GSSWI (R=0.79). Statistically significant strong correlations indicate cumulative Table 1. Evaluation of NDVI3g based onset of growing season (SOS), end of growing season (EOS), length of growing season (LOS), and growing season summed normalized difference vegetation index (GSSNDVI) at site scale. Respective results of spatial (abbreviated as S) and temporal (abbreviated as T) evaluations are indicated by the coefficient of determination (R 2 ) and the correlation coefficient (R).  growing season temperature as the driver of interannual and long-term variations in growing season photosynthetic activity. GSSNDVI variations also agree with those seen in four different productivity proxies from MODIS data ( figure 1(D)). In particular, productivity proxies based on integral NDVI or EVI over the growing season (MOD13C1, MCD43C4 and MCD12Q2) show relatively stronger correlations than model-based GPP estimates (MOD17A3). This could be due to potentially inaccurate parameterizations in the models and/or uncertainties due to additional meteorological forcing data required by such models (Verma et al 2014). The long-term SOS, EOS, LOS and GSSNDVI anomalies reflect the impact of global climate events such as the eruption of Mt. Pinatubo in 1991 (shorter growing season and decreased productivity) (Lucht et al 2002) and the strong El Niño event in 1997-98 (longer growing season and increased productivity) (Buermann et al 2003). A particularly prominent feature in these metrics is the intense photosynthetic activity in NA during 2010, which is about three standard deviations above the mean GSSNDVI. This exceptional anomaly in NA is a consequence of the greatest warmth in 2010 (Blunden et al 2011) and it is also seen in metrics of MODIS data or other previous studies (Friedl et al 2014, Xia et al 2015. These matching characteristics features in metrics inferred from data from different sensors are particularly noteworthy. As for metrics from 36 FLUXNET sites (140 site years), NDVI3g explains 73%, 77% and 82% of variations in SOS, EOS and LOS, respectively (figure 2(A)). NDVI3g SOS and EOS estimates are, on average, 4.2 and 14.6 days later than those inferred from tower GPP data. This translates to a growing season that is 10.5 days longer. Still, NDVI3g data capture the large variation (60-260 days) in growing season across a range of vegetation types (mixed forests, evergreen needleleaf forests, grasses and tundra) (table S1). Similarly, NDVI3g data capture about 80% of the tower based variations in GPP. However, GSSNDVI tends to saturate and shows large variation when GPP is above 1.5 kgC m −2 yr −1 (figure 2(B)). This saturation is a well-known behavior of vegetation index data in dense and productive vegetation types (Sellers 1985, Myneni andWilliams 1994). The saturation has less impact in our study area because only 3.7% of the vegetation exhibits GSSNDVI greater than 150.

Long-term changes in growing season over northern lands
The growing season in the north has lengthened, on an average, by 8.58 days over the past 33 years (2.60 d dec −1 , p<0.05, table 2). The lengthening is greater in EA than in NA (3.04 versus 1.83 d dec −1 , p<0.05). Changes during spring contributed less than changes in autumn to this lengthening in the case of NA. The opposite is the case in EA. Interestingly, changes in growing season duration differed between the first two decades of the data record (1982-1999; 5.06 d dec −1 , p<0.05), which was an exceptionally warm period (Trenberth and Fasullo 2013) and the latter part of the data record (2000-2014; −1.08 d dec −1 , p>0.1) during which a warming hiatus was noted (table 2). This switch from a lengthened (i.e., advancing SOS and delaying EOS) to a shortened (i.e., delaying SOS and advancing EOS) duration was also reported in other studies  inter-annual variations (table S2). Interestingly, at least of the same signs are reported when observed abrupt divergence is ignored.
About 30.6% of northern vegetated land shows statistically significant (Vogelsang's t-PS_T test at 10% significance level) changes in SOS over the past 33 years (figure 3(A) and table S3). A majority of these (27.9%) shows an advancing SOS trend, that is, a trend towards earlier springtime greening. Only 2.7% show the opposite trend. The former is especially pronounced in EA while the latter is seen mostly in boreal NA. However, the degree of advancing trend in SOS over the boreal region (Max. PDFs in EA: −2.8 d dec −1 , NA: −3.0 d dec −1 ) is relatively higher than in the arctic (Max. PDFs EA: −2.5 d yr −1 , NA: −2.3 d dec −1 ). This pattern was reported by previous studies (Shen et al 2014(Shen et al , 2015 and it implies that an earlier SOS in a warmer region may have higher temperature sensitivity than those in a colder region. Reported less sensitive green-up response in arctic vegetation also possibly associates with increasing snowfall in winter/spring time which may hinder much earlier green-up in warmer arctic (e.g., Bieniek et al 2015).
About 21.9% of the study region displays a significant delay in autumn senescence (EOS) over the 33 year period of record (figure 3(B) and table S3). The opposite is seen in about 7.8% of the study area. Boreal regions in both NA and EA show predominant delaying EOS, however, the patterns vary between arctic regions in the two continents. Large proportions (>75% of significant changes) of arctic NA show the delayed EOS trend. In EA, this is observed in only about 25% of the vegetated arctic region showing significant changes. These trends in spring greening and autumn senescence resulted in nearly 33% of the northern vegetation experiencing a lengthened growing season (figure 3(C) and table S3). In most such regions, the longer growing season was due to earlier spring time greening. As shown in figure 3(C), trends in LOS over boreal regions in both continents (Max. PDFs in EA: 3.50 d dec −1 , NA: 3.75 d dec −1 ) have relatively greater lengthening rate than those in arctic regions (Max. PDFs in EA: 2.25 d dec −1 , NA: 3.25 d dec −1 ).

Long-term changes in productivity over northern lands
The analysis indicates that GSSNDVI, a measure of seasonal gross primary productivity, has increased by 2.97 dec −1 (p<0.01) over the circumpolar region. The rate of increase in NA (2.32 dec −1 , p<0.01) is less than in EA (3.34 dec −1 , p<0.01) since the early 1980s (table 2). GSSNDVI exhibits a continuously increasing trend throughout this period, unlike the growing season metrics which show opposite trends between the early (1982-1999) and later (2000-2014) periods of the data record. However, the GSSNDVI trend during the later period (1.87 dec −1 , p>0.1) is lower than in the earlier period (4.23 dec −1 , p>0.05). These results are concordant between AVHRR based NDVI3g data and MODIS NDVI data (table S2).
About 44.4% of the northern vegetated lands exhibit significant changes (p<0.1). 42.0% of the area experience increasing (greening) GSSNDVI trends ( figure 4(A) and table 3). Only a small proportion displays decreasing trend (browning, 2.5%). The greening is more prominently observed in North American mixed forests to the east and arctic coastal tundra and in Eurasian needle-leaf and mixed forests, shrub lands and tundra. A fragmented pattern of greening and browning, mostly over evergreen needle-leaf forest and the forest-shrub ecotone, is seen in the NA boreal region, unlike its counterpart in EA, which shows widespread contiguous greening. This fragmented browning in the interior NA has been reported as consequences of increasing drought stress and fire disturbance (Goetz et al 2005. For a large browning area located in the east Bering coast of Alaska ( figure 4(A)), according to Bieniek et al (2015), this may be linked to delayed snow melt due to increased snow depth in the late winter/early spring as well as increased cloud cover during midsummer. Another large patch of decreasing productivity is prominently seen in the central Siberian plateau (figure 4(A)) which is mostly composed by open larch forest, shrub and erect shrub tundra. This declined productivity is mostly due to the anthropogenic influence (i.e., Cu-Ni smelters) (Toutoubalina and Rees 1999). And smaller areas with such a decline are also found around similar smelters in Kola Peninsula in western part of Russia (Tømmervik et al 2003).  -term (1982-2014) growing season and productivity trends over continental scale. Trends over separated 1982-1999 and 2000-2014 periods are also calculated. The trends were evaluated by Vogelsang's t-PS_T test. CP, NA and EA are for circumpolar, North America and Eurasia regions, respectively. Note: *** : p<0.01, ** : p<0.05, * : p<0.1. Figure 3. Long-term (1982Long-term ( -2014 trends in vegetation growing season onset (SOS, (A)), end (EOS, (B)) and length (LOS, (C)). The trend was calculated using Vogelsang's t-PS_T test at 10% significance level. Non-vegetated pixels and pixels without significant trend were shown in white and gray, respectively. Probability density function (PDF) of change rate per decade for only significant positive and negative changes is also provided (all cases including insignificant or no changes can be found in figure S5). PDFs are normalized to total area showing significant changes in each continent and biome (table S3). NA and EA are for North America and Eurasia. In PDFs, green and red lines represent significant positive and negative changes. Solid and dash lines stand for arctic and boreal regions, respectively.
As shown in figure 4(B), arctic vegetation in both NA and EA (Max. PDFs EA: 5.0% dec −1 , NA: 6.5% dec −1 ) displays relatively greater greening rates (with respect to 1982) than boreal vegetation (Max. PDFs EA: 3.5% dec −1 , NA: 4.0% dec −1 ). The areal proportion of boreal browning is dominant (67.9% of browning area in CP) in the northern lands. In particular, North American boreal vegetation accounts for 55.6% of the browning area in the circumpolar region.
The seasonal maximum value of NDVI (MAX) determines the seasonal trajectory of photosynthetic activity. Thus, examining changes in MAX helps to better understand spatiotemporal changes in GSSNDVI. The spatial distribution of statistically significant increasing trends in MAX, shown in figure S4, closely resembles that of GSSNDVI ( figure 4(A)), especially in the coastal arctic regions. Whereas resembled trend pattern between GSSNDVI and growing season  Long-term (1982Long-term ( -2014 trend in vegetation productivity (GSSNDVI) over Northern vegetated area (A). The trend was calculated using Vogelsang's t-PS_T test at 10% significance level. Non-vegetated pixels and pixels without significant trend were shown in white and gray, respectively. (B) Probability density function (PDF) of GSSNDVI change rate per decade for only showing significant positive and negative changes (All cases including insignificant or no changes can be found in figure S5). PDFs are normalized to total area showing significant changes in each continent and biome (table S3). Green and red lines represent significant positive and negative PDFs. Solid and dash lines stand for arctic and boreal regions, respectively. (C) Trend in spatially aggregated GSSNDVI by grouped vegetation types from 1982 to 2014. Only significant greening and browning pixels were aggregated. For comparison purpose, the GSSNDVIs of all vegetation types were scaled to the GSSNDVI of tundra. NA and EA are for North America and Eurasia, respectively. Details of vegetation groups used in (C) can be found in figure S1. duration (figure 3(C)) can be found in relatively warmer vegetated area. This implies that the seasonal maximum productivity and growing season duration jointly control inter-annual variation and trend of GSSNDVI with differently characterized relative contributions (Xia et al 2015).
Figure 4(C) displays the GSSNDVI time series of four different vegetation groups. The GSSNDVI of forests, other woody vegetation and herbaceous vegetation are scaled to the GSSNDVI of tundra for comparison purposes. All four vegetation groups show increasing GSSNDVI trends with tundra exhibiting the largest trend (8.5% dec −1 ) and forests displaying the lowest (5.5% dec −1 ). This reflects the higher sensitivity of tundra vegetation productivity as compared to boreal forests (Verbyla 2008, Beck andGoetz 2011). There is considerable variation in the trajectory of these time series and the declining greening rate can be clearly seen after the late 1990s. These flattened or slowed change rates are coincident with recently observed warming deceleration (Trenberth and Fasullo 2013) and divergent vegetation growth responses imply differently characterized sensitivities to changing climate. For instance, continued warming may appear to no longer promote boreal forest growth, while the warming may benefit tundra growth Goetz 2011, Bi et al 2013).
Our analysis indicates that 42.0% of the total northern vegetated area shows a greening trend over past three decades (table 3). This translates to a 20.9% gain in productivity since 1982. In contrast, the 2.5% of browning regions resulted in a decrease of about 1.2% of gross primary productivity since 1982. Note that the quantities of productivity changes represent only regions showing significant directional changes. All forests, in particular the mixed and evergreen needleaf forests, contributed significantly to the observed gains in productivity. Equally noteworthy is the contribution of shrublands and the forest-shrub ecotone to productivity gains in view of the large greening extent observed in these vegetation types (table 3).

Concluding remarks
We investigated changes in metrics of growing season (SOS, EOS and LOS) and gross primary productivity (GSSNDVI) over the boreal and arctic lands using long-term satellite observations (GIMMS NDVI3g). An accurate derivation of growing season duration from satellite data is a challenging task. Also, the vegetation index data accumulated over the derived growing season must reflect gross primary productivity. In this sense, the main discriminating point of this work is threefold: (1) this study introduced the photosynthetically active growing season definition by combining optically measured vegetation greenness and ground freeze/thaw data to properly demonstrate photosynthetic activities in northern vegetation.
(2) Moreover, we incorporated yearly varying growing season to productivity characterization, thus we enable  N). Also, total area of each vegetation classes is given (abbreviated as T). For productivity change, productivity increase is abbreviated as I; productivity decrease is abbreviated as D. where NVC i is the total pixel number of the ith vegetation classes showing significant positive or negative changes, T p is the yearly common productivity trend (yr −1 ) of pixel p, A p is the area weight (unitless) of pixel p, G 1982 (=9.04×10 8 , unitless) and B 1982 (=5.89×10 7 , unitless) are the total GSSNDVI of greening pixels and total GSSNDVI of browning pixels in 1982 (table S4). Total vegetated area is about 26.02 million km 2 .
to understand the relative contribution of growing season and peak greenness on annual gross productivity variability.
(3) We evaluated retrieved growing season and productivity metrics using independent multiple direct and indirect measures, particularly eddy-covariance measurements encompassing a wide range of biomes and regions. Special emphasis was placed on assuring that the derived metrics were accurate by comparing them to several independent direct and indirect reference data sets. Overall, these inter-comparison and evaluation analyses reflected that the metrics derived from NDVI3g were reasonably accurate at a range of spatio-temporal scales. Statistical analyses presented in this article provided comprehensive information about patterns in interannual variations and long-term trends over the past three decades. At the hemispheric scale, we observed a significant advance in SOS, delay in EOS and lengthened LOS, all of which are concordant with thermal growing season variations. The longer growing season and increasing photosynthetic activity resulted in a predominant greening trend over 42.0% of the northern vegetated area. This translated to a 20.9% gain in gross primary productivity during last three decades. The GSSNDVI exhibited a continuously increasing trend throughout the 1982-2014 period, unlike the growing season metrics which showed opposite trends between the early (1982)(1983)(1984)(1985)(1986)(1987)(1988)(1989)(1990)(1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999) and later (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014) periods of the data record. The arctic and boreal regions showed surprisingly different variationsgreater rate of productivity change and smaller rate of growing season duration change in the arctic versus the opposite in the boreal vegetation-perhaps reflective of the biome-specific temperature sensitivity of the vegetation. Together these results document largescale spatio-temporal changes happening in the northern vegetation.