Temporal Changes in Coupled Vegetation Phenology and Productivity are Biome-Specific in the Northern Hemisphere

: Global warming has greatly stimulated vegetation growth through both extending the growing season and promoting photosynthesis in the Northern Hemisphere (NH). Analyzing the combined dynamics of such trends can potentially improve our current understanding on changes in vegetation functioning and the complex relationship between anthropogenic and climatic drivers. This study aims to analyze the relationships (long-term trends and correlations) of length of vegetation growing season (LOS) and vegetation productivity assessed by the growing season NDVI integral (GSI) in the NH (>30 ◦ N) to study any dependency of major biomes that are characterized by different imprint from anthropogenic inﬂuence. Spatial patterns of converging/diverging trends in LOS and GSI and temporal changes in the coupling between LOS and GSI are analyzed for major biomes at hemispheric and continental scales from the third generation Global Inventory Monitoring and Modeling Studies (GIMMS) Normalized Difference Vegetation Index (NDVI) dataset for a 32-year period (1982–2013). A quarter area of the NH is covered by converging trends (consistent signiﬁcant trends in LOS and GSI), whereas diverging trends (opposing signiﬁcant trends in LOS and GSI) cover about 6% of the region. Diverging trends are observed mainly in high latitudes and arid/semi-arid areas of non-forest biomes (shrublands, savannas, and grasslands), whereas forest biomes and croplands are primarily characterized by converging trends. The study shows spatially-distinct and biome-speciﬁc patterns between the continental land masses of Eurasia (EA) and North America (NA). Finally, areas of high positive correlation between LOS and GSI showed to increase during the period of analysis, with areas of signiﬁcant positive trends in correlation being more widespread in NA as compared to EA. The temporal changes in the coupled vegetation phenology and productivity suggest complex relationships and interactions that are induced by both ongoing climate change and increasingly intensive human disturbances.


Introduction
Vegetation is of great importance in the interaction of the biosphere and the atmosphere in the Earth system through biophysical (e.g., albedo and water cycle) and biogeochemical (e.g., carbon cycle) feedbacks to modulate regional and global climate [1][2][3][4]. Vegetation metrics like phenology (i.e., the timing of start/end of vegetation growing season and its duration) and photosynthetic activity (commonly measured by the seasonally integrated "greenness" of vegetation) are critical biological characteristics of vegetation, and are sensitive to climate change [5][6][7][8][9][10]. Moreover, even small changes in vegetation photosynthesis could have large effects on the role of terrestrial ecosystems as a carbon sink [11][12][13][14]. Thus, in the context of ongoing global warming, monitoring long-term vegetation activity changes at regional and global scales can help to quantify and understand the effects of climate change on various terrestrial ecosystems and their feedbacks [3,15].
Long-term time-series satellite data have been widely used to study vegetation activity at various spatial and temporal scales [6,[16][17][18][19], such as the NDVI (Normalized Difference Vegetation Index) dataset that is produced by the Global Inventory Monitoring and Modeling Studies (GIMMS) group from the Advanced Very High Resolution Radiometer (AVHRR) of National Oceanic and Atmospheric Administration (NOAA). The GIMMS3g NDVI (starting in July 1981) has been widely used as a proxy for vegetation productivity by the global change research community. GIMMS NDVI-based studies have shown that the length of vegetation growing season (LOS) has increased in most areas of the Northern Hemisphere (NH) from 1982 to the late 1990s through an advancing start of growing season (SOS) and/or delaying end of growing season (EOS) due to climate change [6,[16][17][18]20,21]. However, since the end of the 1990s, this advancing trend of SOS has been shown to be weakened, and even reversed, in some regions [6,7,14,[22][23][24][25][26], although there is no consensus about when the reversal happened. On the other hand, long-term time-series NDVI dataset from AVHRR have also revealed significant greening trends of vegetation (based on seasonal NDVI integrals) in the NH from the early 1980s to the late 1990s [7,17,18,20], which have been interpreted as an increase in photosynthesis and vegetation productivity through terrestrial ecosystem models [20,[27][28][29]. Thus, the observed 'northern greening' has been central for an increasing carbon sink in the northern ecosystems over the past 30 years [16,[30][31][32][33][34].
A linkage between changes of LOS and changes in vegetation productivity would be expected, as changing LOS implies changing the duration available for vegetation photosynthesis. In the context of global warming, this changing ecosystem functioning is observed and is expected to impact biomass accumulation of terrestrial ecosystem, and thereby regional and global carbon balances [35][36][37][38][39][40][41][42]. It has been suggested that longer LOS may extend the period of net carbon sequestration and eventually enhance vegetation productivity in northern ecosystems [38]. Moreover, Richardson et al. [41] found that forest phenology derived from carbon flux measurements was closely correlated to both gross and net ecosystem productivity. However, remotely sensed based studies of changes in vegetation phenology and biomass have primarily been conducted independently; thereby disregarding that synthesizing the two research themes might reveal novel information on changes in plant functional traits, ultimately to be associated with imprints from anthropogenic land management. Whether a given area is characterized by either converging or diverging trends in LOS and productivity may be induced by the interplay between human disturbances (land use and land cover changes) and climatic drivers. Analyzing the combined spatio-temporal dynamics of such trends can therefore improve our current understanding of the expected complex relationship between anthropogenic and climatic drivers on changes in vegetation functioning during recent decades. We focused on northern hemisphere ecosystems (>30 • N) (Figure 1), which are characterized by vegetation with clear seasonality. The response of vegetation to global warming is highly sensitive in this region, and northern ecosystems are expected to continue growing at a faster rate and thus accumulating more carbon in a warmer future [16,20]. This study investigates the relationships (long-term trends and correlations) between LOS and vegetation productivity in the NH (>30 • N) using third generation GIMMS NDVI dataset for a 32-year period (1982-2013) by analyzing converging/diverging patterns in trends in LOS and the growing season NDVI integral as a function of major land biomes at hemispheric and continental scales.

Study Area and Biomes
Eight major land biomes were aggregated from the Land Cover Type 2 (University of Maryland (UMD) scheme) of the MODIS MCD12C1 Land Cover Type Climate Modeling Grid product in 2001 at a spatial resolution of 0.05°. This product is derived from observations spanning a year's input of MODIS Terra and Aqua data using the same algorithm that produces the MCD12Q1 Global 500 m Land Cover Type product. To examine the difference of changes between natural and human managed vegetation systems, this study included croplands, although its seasonal variability is greatly influenced by human disturbances. This study analyzed only pixels with a mean annual NDVI of greater than 0.1 for the GIMMS3g NDVI dataset, since lower NDVI values usually indicate sparse vegetation, which are more likely to be contaminated by the signal from bare soil or atmospheric/sensor related noise [6,18].

NDVI Data
The GIMMS3g NDVI dataset [43] is analyzed for the period of January 1982 to December 2013. Data is provided in 1/12 degree resolution and is composited into 15-day periods using the Maximum Value Compositing (MVC) technique [44]. The data were obtained from the AVHRR instrument onboard the NOAA satellite series 7, 9, 11, 14, 16, 17, 18, and 19. These data have been calibrated for sensor shifts, and corrected to remove the effects of sensor degradation, satellite orbit drift, solar zenith angles, and other factors, such as the effects of stratospheric aerosol loadings from the El Chichon and Mount Pinatubo eruptions in April 1984 and June 1991 [18,45]. This dataset has been widely used for detecting vegetation phenological and photosynthetic changes, and the current version has improved data quality especially in high latitudes with shorter growing season when compared to previous versions of the GIMMS NDVI datasets [46].

Estimation of Vegetation Phenology and Productivity Metrics
The phenology metrics of LOS, SOS, and EOS of GIMMS3g are calculated per pixel from time series parameterization by a Savitsky-Golay filter in the TIMESAT software [47,48]. The growing season NDVI integral (GSI) is used here as a proxy for vegetation productivity (greenness) [49], and is also computed from TIMESAT. The advantage of this approach is to minimize any influence from year-to-year variations in the length of seasonal snow cover on the calculated NDVI trends [49]. The Savitsky-Golay filter is a moving filter to fit values from least square to polynomial [48] and the degree of adaptation to the upper envelope can be tuned to the desired level [50]. The start and end of the growing season (SOS and EOS) are determined, respectively, for individual years from the per-pixel polynomial fit based on a parameterization of the fitted seasonal NDVI curve. The

Study Area and Biomes
Eight major land biomes were aggregated from the Land Cover Type 2 (University of Maryland (UMD) scheme) of the MODIS MCD12C1 Land Cover Type Climate Modeling Grid product in 2001 at a spatial resolution of 0.05 • . This product is derived from observations spanning a year's input of MODIS Terra and Aqua data using the same algorithm that produces the MCD12Q1 Global 500 m Land Cover Type product. To examine the difference of changes between natural and human managed vegetation systems, this study included croplands, although its seasonal variability is greatly influenced by human disturbances. This study analyzed only pixels with a mean annual NDVI of greater than 0.1 for the GIMMS3g NDVI dataset, since lower NDVI values usually indicate sparse vegetation, which are more likely to be contaminated by the signal from bare soil or atmospheric/sensor related noise [6,18].

NDVI Data
The GIMMS3g NDVI dataset [43] is analyzed for the period of January 1982 to December 2013. Data is provided in 1/12 degree resolution and is composited into 15-day periods using the Maximum Value Compositing (MVC) technique [44]. The data were obtained from the AVHRR instrument onboard the NOAA satellite series 7, 9, 11, 14, 16, 17, 18, and 19. These data have been calibrated for sensor shifts, and corrected to remove the effects of sensor degradation, satellite orbit drift, solar zenith angles, and other factors, such as the effects of stratospheric aerosol loadings from the El Chichon and Mount Pinatubo eruptions in April 1984 and June 1991 [18,45]. This dataset has been widely used for detecting vegetation phenological and photosynthetic changes, and the current version has improved data quality especially in high latitudes with shorter growing season when compared to previous versions of the GIMMS NDVI datasets [46].

Estimation of Vegetation Phenology and Productivity Metrics
The phenology metrics of LOS, SOS, and EOS of GIMMS3g are calculated per pixel from time series parameterization by a Savitsky-Golay filter in the TIMESAT software [47,48]. The growing season NDVI integral (GSI) is used here as a proxy for vegetation productivity (greenness) [49], and is also computed from TIMESAT. The advantage of this approach is to minimize any influence from year-to-year variations in the length of seasonal snow cover on the calculated NDVI trends [49]. The Savitsky-Golay filter is a moving filter to fit values from least square to polynomial [48] and the degree of adaptation to the upper envelope can be tuned to the desired level [50]. The start and end of the growing season (SOS and EOS) are determined, respectively, for individual years from the per-pixel polynomial fit based on a parameterization of the fitted seasonal NDVI curve. The parameterization is calculated from a percentage NDVI threshold of total seasonal NDVI amplitude. The LOS is determined by the duration between the start and end of the growing season in each year. The GSI representing growing season vegetation productivity is calculated by the area between the fitted function and the base level from the start to end of the growing season. The parameters that are applied for the TIMESAT analysis are as follows: seasonal parameter = 0.5, number of envelope iterations = 2, adaptation strength = 1, Savitzky-Golay window size = 2, amplitude season start and end = 20%. For the case of croplands with multiple seasonal peaks, we have chosen to make the fitting procedure always output a primary season with maximum annual amplitude.

Linear Trend Estimation
Pearson Product-moment linear correlation test on trends is employed to estimate the direction and strength of linear development in the long-term vegetation metrics (i.e., LOS, SOS, EOS, and GSI). We used a simple least square linear regression model with time as the independent variable and vegetation metrics derived from NDVI as the dependent variables. Linear correlation maps out the Pearson Product-moment linear correlation coefficients (r-values ranging from −1 to 1) between the values of each pixel over time and a perfectly linear series [51].

Trend Comparison
Significance maps (r-values) of LOS and GSI time series trends were reclassified into seven classes: positive/negative trends at confidence levels of 90% (p < 0.1), 95% (p < 0.05), and 99% (p < 0.01), and not significant (p > 0.1) ( Figure 2C,D). These maps were combined into a single trend difference map by aggregating all significant positive and negative trends, respectively, to a single class covering pixels that are significant at the 90% level (p < 0.1) (one for positive and one for negative trends) before being combined, yielding nine possible combinations ( Figure 3 and Table 1. Specifically, combinations of consistent significant positive/negative trends in LOS and GSI are defined as converging trends, whereas combinations of opposing significant trends in LOS and GSI are termed diverging trends (Table 1). Additionally, a distinction is made between converging trends type 1 (both trends significantly positive) and converging trends type 2 (both trends significantly negative). The same distinction also applies for diverging trends. A significance level of 90% was chosen to maintain the clearest possible spatial patterns of regional-scale clusters of pixels characterized by different trend combinations. The total numbers of pixels of areas with significant trends are summarized for major land biomes at hemispheric and continental scales.
parameterization is calculated from a percentage NDVI threshold of total seasonal NDVI amplitude. The LOS is determined by the duration between the start and end of the growing season in each year. The GSI representing growing season vegetation productivity is calculated by the area between the fitted function and the base level from the start to end of the growing season. The parameters that are applied for the TIMESAT analysis are as follows: seasonal parameter = 0.5, number of envelope iterations = 2, adaptation strength = 1, Savitzky-Golay window size = 2, amplitude season start and end = 20%. For the case of croplands with multiple seasonal peaks, we have chosen to make the fitting procedure always output a primary season with maximum annual amplitude.

Linear Trend Estimation
Pearson Product-moment linear correlation test on trends is employed to estimate the direction and strength of linear development in the long-term vegetation metrics (i.e., LOS, SOS, EOS, and GSI). We used a simple least square linear regression model with time as the independent variable and vegetation metrics derived from NDVI as the dependent variables. Linear correlation maps out the Pearson Product-moment linear correlation coefficients (r-values ranging from −1 to 1) between the values of each pixel over time and a perfectly linear series [51].

Trend Comparison
Significance maps (r-values) of LOS and GSI time series trends were reclassified into seven classes: positive/negative trends at confidence levels of 90% (p < 0.1), 95% (p < 0.05), and 99% (p < 0.01), and not significant (p > 0.1) ( Figure 2C,D). These maps were combined into a single trend difference map by aggregating all significant positive and negative trends, respectively, to a single class covering pixels that are significant at the 90% level (p < 0.1) (one for positive and one for negative trends) before being combined, yielding nine possible combinations ( Figure 3 and Table 1. Specifically, combinations of consistent significant positive/negative trends in LOS and GSI are defined as converging trends, whereas combinations of opposing significant trends in LOS and GSI are termed diverging trends (Table 1). Additionally, a distinction is made between converging trends type 1 (both trends significantly positive) and converging trends type 2 (both trends significantly negative). The same distinction also applies for diverging trends. A significance level of 90% was chosen to maintain the clearest possible spatial patterns of regional-scale clusters of pixels characterized by different trend combinations. The total numbers of pixels of areas with significant trends are summarized for major land biomes at hemispheric and continental scales.

Estimation of Correlations between LOS and GSI and Their Trends
All of the vegetation metrics were detrended to avoid possible spurious regressions of variables with a co-varying trend to allow for focused analysis of correlations in interannual variability. The

Estimation of Correlations between LOS and GSI and Their Trends
All of the vegetation metrics were detrended to avoid possible spurious regressions of variables with a co-varying trend to allow for focused analysis of correlations in interannual variability. The

Estimation of Correlations between LOS and GSI and Their Trends
All of the vegetation metrics were detrended to avoid possible spurious regressions of variables with a co-varying trend to allow for focused analysis of correlations in interannual variability. The linear trend derived from the least squares method was removed for both LOS and GSI time series within the entire research period (1982-2013) and each 16-year window (sub-period analysis described below) [9]. Correlations between LOS and GSI were assessed by means of parametric correlations using the Pearson coefficient for (1) the entire period and (2) for 17 sub-periods calculated from a 16-year window to be able to detect changes on correlation over time. A 16-year period was selected as a compromise between sufficient numbers of years to calculate the correlation on the one hand, and to have enough observations (17) to study changes over time on the other hand. Finally, the trends in the correlations between LOS and GSI were calculated from the 17 time-series correlation maps (for each 16-year window) by linear temporal trend analysis.

Trends in Vegetation Phenology and Growing Season Integral
More pixels of the NH with significant linear trends in LOS ( Figure 2C) (42.45%) were observed as compared to SOS (Figure 2A) (34.83%), and EOS ( Figure 2B) (41.37%) ( Table 2). Moreover, most areas show opposite spatial patterns in trends in SOS and LOS (Figure 2A,C), but consistent pattern in EOS and LOS trends ( Figure 2B,C). The quantitative assessment and spatial patterns indicate that extending LOS in the NH is caused by both advancing SOS and delaying EOS in most regions, whereas shortening LOS is due to both delaying SOS and advancing EOS.  Areas of pixels of the NH with significant positive trends in GSI ( Figure 2D) (51.22%) are considerably larger than those in LOS (23.24%) ( Table 2), whereas the areas of pixels with negative trends in GSI (9.95%) are much smaller than those in LOS (19.21%) ( Table 2). This suggests that areas of significant greening trends are dominating the NH, with LOS appearing comparably mixed in lengthening and shortening trends over the past three decades, indicating that some areas are experiencing increasing vegetation productivity without significant changes in LOS.
At the continental scale, North America (NA) is dominated by significant decreasing trends in LOS, but Eurasia (EA) is dominated by increasing trends ( Figure 2C and Table 2). This pattern can be explained by the predominance of significant increasing trends in SOS in NA (28.05%) combined with relatively large areas of decreasing trends in EOS (21.56%) (comparable to increasing trends (23.8%) in EOS) ( Table 2). In contrast, EA is dominated by significant decreasing trends in SOS (25.68%), and has slightly larger areas with positive trends in EOS (20.43%) when compared to negative trends (19.11%). For GSI ( Figure 2D), the areas of pixels with significant positive trends are a little smaller in NA than those of EA, whereas pixels of NA with negative trends are slightly larger than those in EA ( Table 2). The difference between pixels of GSI with significant positive and negative trends in EA (53.52%/8.30%) is considerably larger than that in NA (46.22%/13.51%), suggesting that the spatial pattern of trends in vegetation productivity is more divergent in NA than in EA.

Converging/Diverging Trends in LOS and GSI
A direct comparison of trends in LOS and GSI of GIMMS3g (1982-2013) ( Figure 3) is made for the pixels of significant trends in Figure 2C,D. Pixels with converging trends in LOS and GSI (consistent significant positive/negative trends in LOS and GSI) (24.53%, (a + d) Table 1) cover considerably larger areas than pixels with diverging trends (opposite significant trends in LOS and GSI) (5.78%, (b + c) Table 1), indicating that the NH is dominated by converging trends in LOS and GSI. However, areas of pixels of NA with converging trends (18.70%) are considerably smaller than those of EA (27.22%), whereas areas of pixels of NA with diverging trends (11.08%) are considerably larger than those of EA (3.33%) ( Table 1). This suggests that the vegetation metrics in NA are characterized by more diverging trends than EA. More specifically, areas of pixels of the NH with consistent significant positive trends in LOS and GSI (converging trends type 1) (19.46%, (d) Table 1) are considerably larger, when compared to pixels with consistent significant negative trends (converging trends type 2) (5.07%, (a) Table 1), indicating that the NH is dominated by converging trends type 1. However, areas of pixels of NA with converging trends type 1 (12.23%) are considerably smaller than those of EA (22.80%), whereas areas of pixels in NA with converging trends type 2 are somewhat larger than those of EA ( Table 1), suggesting that the dominance of converging trends type 1 is more pervasive in EA. In contrast, only 0.41% pixels of the NH are characterized by significant positive trends in LOS/significant negative trends in GSI (diverging trends type 1) ((c) Table 1), whereas the opposite pattern (diverging trends type 2) ((b) Table 1) covers considerably larger areas (5.37%), suggesting that diverging trends in the NH are dominated by diverging trends type 2. At the continental scale, areas of pixels of NA with diverging trends type 1 are somewhat larger than those of EA (Table 1), while pixels of NA with diverging trends type 2 (10.28%) are considerably larger than those of EA (3.11%), thus indicating that the dominance of diverging trends type 2 is more pervasive in NA.
Spatially, large areas with converging trends type 1 are mainly found in central and eastern NA, northern and eastern Europe, north and south Russia, north Mongolia, and north and northeast China ((+1, +1) in Figure 3). On the contrary, areas with converging trends type 2 are primarily observed in northeast Canada, central America, east and central Russia and northeast Tibetan Plateau ((−1, −1) in Figure 3). In contrast, areas with diverging trends in LOS and GSI can be seen in high latitudes of NA and EA, central NA, northern Europe, Kazakhstan, and eastern Russia ((−1, +1) and (+1, −1) in Figure 3), suggesting that diverging trends are mainly in high latitudes and arid/semi-arid areas.

Converging/Diverging Trends for Major Biomes
High percentages of pixels with converging trends (C in Figure 4A) are mainly observed for all forest biomes (highest for mixed forest) and croplands, whereas low percentages of pixels with converging trends are primarily seen for non-forest biomes (shrublands, savannas, and grasslands). In contrast, high percentages of pixels with diverging trends (D in Figure 4A) are mostly observed in non-forest biomes, whereas low percentages of pixels with diverging trends are mainly observed for evergreen needleleaf forest, deciduous broadleaf forest, mixed forest and croplands ( Figure 4A), suggesting that diverging trends are predominant in non-forest biomes and are less frequently observed for biomes of forest and croplands. However, at the continental scale, converging trends are more frequent for biomes of deciduous broadleaf forest, mixed forest, and croplands, and less for biomes of non-forest and needleleaf forest (both evergreen and deciduous) in NA, whereas EA has the same pattern as for the NH ( Figure 4A). Diverging trends are more frequent for biomes of shrublands and grasslands in NA and less common for forest biomes, while in EA more diverging trends are observed for biomes of non-forest and needleleaf forest as compared to biomes of deciduous broadleaf forest, mixed forest, and croplands ( Figure 4A). Moreover, the high ratio (Ratio C/D, Figure 4B) indicates that areas with converging trends are considerably larger than areas with diverging trends for all biomes. Non-forest biomes show the lowest ratios in the NH and both continents, whereas the highest ratios are observed for evergreen needleleaf forest, deciduous broadleaf forest, mixed forest, and croplands in the NH, EA, and all forest biomes in NA ( Figure 4B). This further suggests that non-forest biomes are characterized by more diverging trends, when compared to biomes of forest and croplands primarily showing converging trends.
Remote Sens. 2017, 9,1277 8 of 18 are more frequent for biomes of deciduous broadleaf forest, mixed forest, and croplands, and less for biomes of non-forest and needleleaf forest (both evergreen and deciduous) in NA, whereas EA has the same pattern as for the NH ( Figure 4A). Diverging trends are more frequent for biomes of shrublands and grasslands in NA and less common for forest biomes, while in EA more diverging trends are observed for biomes of non-forest and needleleaf forest as compared to biomes of deciduous broadleaf forest, mixed forest, and croplands ( Figure 4A). Moreover, the high ratio (Ratio C/D, Figure 4B) indicates that areas with converging trends are considerably larger than areas with diverging trends for all biomes. Non-forest biomes show the lowest ratios in the NH and both continents, whereas the highest ratios are observed for evergreen needleleaf forest, deciduous broadleaf forest, mixed forest, and croplands in the NH, EA, and all forest biomes in NA ( Figure 4B). This further suggests that non-forest biomes are characterized by more diverging trends, when compared to biomes of forest and croplands primarily showing converging trends.  Converging trends type 1 (C1 in Figure 4C) is primarily observed in mixed forest, croplands, and shrublands for the NH and EA (mixed forest, croplands for NA), while converging trends type 2 (C2 in Figure 4C) is mainly found in non-forest biomes for the NH, NA, and EA. In contrast, the two types of diverging trends are both dominating in non-forest biomes for the NH (Figure 4E). At the continental scale, areas with diverging trends type 2 are primarily observed in shrublands and grasslands for NA and in non-forest biomes for EA, whereas areas with diverging trends type 1 are mainly found in non-forest biomes for NA and in shrublands and grasslands for EA.
High percentages of pixels with converging trends type 1 per biome (C1 in Figure 4D) are mainly observed for evergreen needleleaf forest, deciduous broadleaf forest, mixed forest, and croplands in the NH, NA, and EA, whereas high percentages of pixels with converging trends type 2 (C2 in Figure 4D) are primarily found for deciduous needleleaf forest and non-forest biomes in the NH and both continents. This suggests that biomes of forest and croplands are dominated by converging trends type 1, whereas converging trends type 2 are pervasive in deciduous needleleaf forest and non-forest biomes. In the case of diverging trends, high percentages of pixels with diverging trends type 2 (D2 in Figure 4F) are mainly seen for deciduous needleleaf forest and non-forest biomes in the NH, while high percentages of pixels with diverging trends type 1 (D1 in Figure 4F) are primarily found for non-forest biomes in the NH and NA.

Correlations between LOS and GSI (1982-2013)
The interannual variation of GSI is extensively positively correlated with LOS in the NH ( Figure 5) ( Table 3, 83.09%). Nevertheless, areas with low (non-significant) and significant negative correlation can be observed, primarily in high latitudes/arctic region of NA and EA, central NA, central Asia, Mongolia, and north China, indicating that LOS and GSI are less or negatively correlated in high latitudes/arctic and arid/semi-arid regions. However, differences in the spatial coverage of these correlations between LOS and GSI are distinct between NA and EA (Table 3). EA shows larger areas with significant positive correlation (85.31%) than those of NA (78.29%), whereas pixels with negative correlation in NA (1.73%) are larger than those in EA (0.44%).
Remote Sens. 2017, 9, 1277 9 of 18 Converging trends type 1 (C1 in Figure 4C) is primarily observed in mixed forest, croplands, and shrublands for the NH and EA (mixed forest, croplands for NA), while converging trends type 2 (C2 in Figure 4C) is mainly found in non-forest biomes for the NH, NA, and EA. In contrast, the two types of diverging trends are both dominating in non-forest biomes for the NH (Figure 4E). At the continental scale, areas with diverging trends type 2 are primarily observed in shrublands and grasslands for NA and in non-forest biomes for EA, whereas areas with diverging trends type 1 are mainly found in non-forest biomes for NA and in shrublands and grasslands for EA.
High percentages of pixels with converging trends type 1 per biome (C1 in Figure 4D) are mainly observed for evergreen needleleaf forest, deciduous broadleaf forest, mixed forest, and croplands in the NH, NA, and EA, whereas high percentages of pixels with converging trends type 2 (C2 in Figure  4D) are primarily found for deciduous needleleaf forest and non-forest biomes in the NH and both continents. This suggests that biomes of forest and croplands are dominated by converging trends type 1, whereas converging trends type 2 are pervasive in deciduous needleleaf forest and non-forest biomes. In the case of diverging trends, high percentages of pixels with diverging trends type 2 (D2 in Figure 4F) are mainly seen for deciduous needleleaf forest and non-forest biomes in the NH, while high percentages of pixels with diverging trends type 1 (D1 in Figure 4F) are primarily found for non-forest biomes in the NH and NA.

Correlations between LOS and GSI (1982-2013)
The interannual variation of GSI is extensively positively correlated with LOS in the NH ( Figure  5) (Table 3, 83.09%). Nevertheless, areas with low (non-significant) and significant negative correlation can be observed, primarily in high latitudes/arctic region of NA and EA, central NA, central Asia, Mongolia, and north China, indicating that LOS and GSI are less or negatively correlated in high latitudes/arctic and arid/semi-arid regions. However, differences in the spatial coverage of these correlations between LOS and GSI are distinct between NA and EA (Table 3). EA shows larger areas with significant positive correlation (85.31%) than those of NA (78.29%), whereas pixels with negative correlation in NA (1.73%) are larger than those in EA (0.44%).  Pixels with significant positive correlation between LOS and GSI ( Figure 6A) are mostly distributed in biomes of shrublands, grasslands, and croplands, but are rarely observed in any of the  Pixels with significant positive correlation between LOS and GSI ( Figure 6A) are mostly distributed in biomes of shrublands, grasslands, and croplands, but are rarely observed in any of the forest biomes, except for mixed forest. In contrast, pixels with significant negative correlation are primarily distributed in non-forest biomes for the NH and NA (shrublands and grasslands for EA), but are rarely observed in forest biomes and croplands.
forest biomes, except for mixed forest. In contrast, pixels with significant negative correlation are primarily distributed in non-forest biomes for the NH and NA (shrublands and grasslands for EA), but are rarely observed in forest biomes and croplands. Biomes having the lowest percentages of pixels with significant positive correlation between LOS and GSI ( Figure 6B) are mainly observed for non-forest biomes for the NH and NA (shrublands and grasslands for EA), whereas forest biomes and croplands show the highest percentages for the NH, NA, and EA. On the contrary, biomes having the highest percentages of pixels with significant negative correlation between LOS and GSI are primarily found in non-forest biomes for the NH and NA, whereas forest biomes and croplands show the lowest percentages.

Trends in Correlations between LOS and GSI
Linear trends in correlations between LOS and GSI calculated from correlations between LOS and GSI for each 16-year window from 1982-1997 to 1998-2013 (Figure 7) show larger areas (39.76%) of significant increasing strength of the correlations between LOS and GSI than areas with decreasing correlations (31.14%) for the NH. It is evident that most areas in high latitudes of NA show positive trends and areas of pixels with significant positive trends in NA (46.92%) are considerably larger than those with negative trends (25.92%). Contrastingly, areas of pixels with a significant positive trend in EA (36.45%) are only a little larger than those with negative trends (33.55%) ( Table 4). This suggests that EA shows more divergent pattern of trends in correlations between LOS and GSI, as compared to NA with dominating increasing trends.  Biomes having the lowest percentages of pixels with significant positive correlation between LOS and GSI ( Figure 6B) are mainly observed for non-forest biomes for the NH and NA (shrublands and grasslands for EA), whereas forest biomes and croplands show the highest percentages for the NH, NA, and EA. On the contrary, biomes having the highest percentages of pixels with significant negative correlation between LOS and GSI are primarily found in non-forest biomes for the NH and NA, whereas forest biomes and croplands show the lowest percentages.

Trends in Correlations between LOS and GSI
Linear trends in correlations between LOS and GSI calculated from correlations between LOS and GSI for each 16-year window from 1982-1997 to 1998-2013 (Figure 7) show larger areas (39.76%) of significant increasing strength of the correlations between LOS and GSI than areas with decreasing correlations (31.14%) for the NH. It is evident that most areas in high latitudes of NA show positive trends and areas of pixels with significant positive trends in NA (46.92%) are considerably larger than those with negative trends (25.92%). Contrastingly, areas of pixels with a significant positive trend in EA (36.45%) are only a little larger than those with negative trends (33.55%) ( Table 4). This suggests that EA shows more divergent pattern of trends in correlations between LOS and GSI, as compared to NA with dominating increasing trends.
Remote Sens. 2017, 9,1277 10 of 18 forest biomes, except for mixed forest. In contrast, pixels with significant negative correlation are primarily distributed in non-forest biomes for the NH and NA (shrublands and grasslands for EA), but are rarely observed in forest biomes and croplands. Biomes having the lowest percentages of pixels with significant positive correlation between LOS and GSI ( Figure 6B) are mainly observed for non-forest biomes for the NH and NA (shrublands and grasslands for EA), whereas forest biomes and croplands show the highest percentages for the NH, NA, and EA. On the contrary, biomes having the highest percentages of pixels with significant negative correlation between LOS and GSI are primarily found in non-forest biomes for the NH and NA, whereas forest biomes and croplands show the lowest percentages.

Trends in Correlations between LOS and GSI
Linear trends in correlations between LOS and GSI calculated from correlations between LOS and GSI for each 16-year window from 1982-1997 to 1998-2013 (Figure 7) show larger areas (39.76%) of significant increasing strength of the correlations between LOS and GSI than areas with decreasing correlations (31.14%) for the NH. It is evident that most areas in high latitudes of NA show positive trends and areas of pixels with significant positive trends in NA (46.92%) are considerably larger than those with negative trends (25.92%). Contrastingly, areas of pixels with a significant positive trend in EA (36.45%) are only a little larger than those with negative trends (33.55%) ( Table 4). This suggests that EA shows more divergent pattern of trends in correlations between LOS and GSI, as compared to NA with dominating increasing trends.    Pixels with significant positive trends in correlations between LOS and GSI ( Figure 8A) are mainly distributed in croplands and non-forest biomes for the NH, but are rarely found in forest biomes, except for mixed forest. Interestingly, a similar distribution pattern is observed for pixels with negative trends. Biomes having high percentages of pixels with significant positive trends in correlations between LOS and GSI ( Figure 8B) are mainly observed for savannas, shrublands, evergreen needleaf forest, and mixed forest for the NH, whereas low percentages of these trends are primarily found for deciduous broadleaf and needleleaf forest, croplands, and grasslands. Complementarily, biomes having high percentages of pixels with significant negative trends in correlations between LOS and GSI ( Figure 8B) are mostly observed for deciduous broadleaf and needleleaf forest, croplands, and grasslands, while low percentages of these trends are mainly found for savannas, shrublands, evergreen needleleaf forest, and mixed forest.  Pixels with significant positive trends in correlations between LOS and GSI ( Figure 8A) are mainly distributed in croplands and non-forest biomes for the NH, but are rarely found in forest biomes, except for mixed forest. Interestingly, a similar distribution pattern is observed for pixels with negative trends. Biomes having high percentages of pixels with significant positive trends in correlations between LOS and GSI ( Figure 8B) are mainly observed for savannas, shrublands, evergreen needleaf forest, and mixed forest for the NH, whereas low percentages of these trends are primarily found for deciduous broadleaf and needleleaf forest, croplands, and grasslands. Complementarily, biomes having high percentages of pixels with significant negative trends in correlations between LOS and GSI ( Figure 8B) are mostly observed for deciduous broadleaf and needleleaf forest, croplands, and grasslands, while low percentages of these trends are mainly found for savannas, shrublands, evergreen needleleaf forest, and mixed forest.

Trends in Vegetation Phenology and Growing Season Integral
The calculated trends in LOS, including SOS and EOS (Figure 2A-C), are consistent with the spatial patterns of phenological trends in Zhao et al. [52], who extracted vegetation phenological changes above 40°N of the NH during 1982-2013 from GIMMS3g NDVI data. They also employed the TIMESAT approach, but with different parameter settings as compared to this study, indicating the robustness of the method for phenology retrieval and consequently the results obtained in this study. Moreover, the calculated trends in GSI ( Figure 2D) are in agreement with the result from Fensholt et al. [53], although their study period  was two years shorter than this study.
This study suggests that EA is dominated by significant increasing trends in LOS, whereas decreasing trends are more common in NA ( Figure 2C and Table 2). These results are supported by Barichivich et al. [54] who showed that the overall increasing trend in LOS was stronger and more significant in EA than that in NA. Garonna et al. [55] found that areas with greatest growing season lengthening are concentrated in pockets of continental Europe, western Russia, and southern Fennoscandia, whereas regions with growing season shortening are mainly found in the Po valley,

Trends in Vegetation Phenology and Growing Season Integral
The calculated trends in LOS, including SOS and EOS (Figure 2A-C), are consistent with the spatial patterns of phenological trends in Zhao et al. [52], who extracted vegetation phenological changes above 40 • N of the NH during 1982-2013 from GIMMS3g NDVI data. They also employed the TIMESAT approach, but with different parameter settings as compared to this study, indicating the robustness of the method for phenology retrieval and consequently the results obtained in this study. Moreover, the calculated trends in GSI ( Figure 2D) are in agreement with the result from Fensholt et al. [53], although their study period  was two years shorter than this study.
This study suggests that EA is dominated by significant increasing trends in LOS, whereas decreasing trends are more common in NA ( Figure 2C and Table 2). These results are supported by Barichivich et al. [54] who showed that the overall increasing trend in LOS was stronger and more significant in EA than that in NA. Garonna et al. [55] found that areas with greatest growing season lengthening are concentrated in pockets of continental Europe, western Russia, and southern Fennoscandia, whereas regions with growing season shortening are mainly found in the Po valley, western France, and around the Caspian Sea for Europe from 1982 to 2011, which is also consistent with the result of this study ( Figure 2C). However, the findings that an increasing LOS is due to both advancing SOS and delaying EOS for most areas, but decreasing LOS is caused by both delaying SOS and advancing EOS in most regions are contradictory to Wu et al. [56] who showed that EOS was positively correlated with SOS in China's temperate regions during 1999-2013 (using SPOT-VGT data). This difference may be explained by the different periods of analysis and data sources between the two studies. Moreover, the observed dominance of greening trends in EA ( Figure 2D) is also supported by Piao et al. [26], who reported a significant positive trend in average growing season NDVI for temperate and boreal EA during 1982-2006.
It has been shown that changes in ecosystem properties can be better characterized by separating the whole study period into individual time-periods representing specific vegetation conditions through time [57,58]. From such approach, both increasing LOS and vegetation greenness have been shown to be slowing down or have experienced abrupt changes around 2000 in the NH using break point analysis techniques [6,7,26,52,[58][59][60]. Nevertheless, a recent study by Wang et al. [61] suggested that the advancing SOS is unlikely to have slowed down or have abruptly changed at the hemispheric scale of the NH over the last three decades, although it could have happened at a local/regional scale because of changes in fire regimes or winter chilling.

Converging/Diverging Trends and Correlations between LOS and GSI
This study shows that a quarter area of the NH is covered by dominant converging trends in LOS and GSI, although EA shows more areas of convergence than NA ( Figure 3). Meanwhile, GSI is predominantly positively correlated with LOS in the NH, with significant positive correlations being more widespread in EA than in NA ( Figure 5). Regions with converging trends in LOS and GSI are generally corresponding to areas with high positive correlations between LOS and GSI, which could be expected. However, the fact that trends in LOS and GSI are of similar direction does not a priori cause a high correlation as it depends on the magnitude and statistical strengths of the individual trends. This pattern is in accordance with the general assumption that increasing LOS that is caused by pervasive warming in middle and high latitudes of the NH has led to increasing vegetation productivity, simply because more days are available for carbon assimilation and biomass accumulation [38,[40][41][42][62][63][64]. For example, Dragoni et al. [65] showed that LOS in a temperate deciduous forest in Midwestern US extended at the rate of about three days/year during 1998-2008 (because of delayed autumn senescence) was accompanied by an increasing trend of annual net C uptake. A comparable result for a temperate deciduous forest in Denmark was reported by Pilegaard et al. [66], who showed that an increasing trend in C uptake during the period of 1996-2009 was partially attributed to a concurrent increasing trend in the duration of photosynthetic activity. Moreover, highly consistent with the converging trends that were shown in this study (Figure 3), Keenan et al. [67] found that carbon uptake through photosynthesis has increased in temperate forests of eastern US owing to warming-induced increasing LOS.
Diverging trends in LOS and GSI only cover about 6% areas of the NH, mainly in high latitudes and arid/semi-arid regions (Figure 3), and larger areas of divergence are observed in NA than EA. Additionally, GSI is characterized by low or negative correlations, with LOS primarily in high latitudes and arid/semi-arid areas and regions of significant negative correlations are more frequent in NA than in EA ( Figure 5). This indicates that areas with diverging trends are corresponding to regions with low or negative correlations between LOS and GSI in the NH. This pattern may be due to distinct responses of vegetation phenology and productivity to the ongoing warming in high northern latitudes [54,[68][69][70]. In agreement with the result from this study, Park et al. [70] found that vegetation productivity has continuously increased during 2000-2014 compared to decreasing LOS in the boreal and arctic regions of NH. Moreover, Ivits et al. [71] showed that vegetation productivity and phenology reacted differently to drought events depending on ecosystem and land cover at the continental scale of Europe. This was further confirmed by Xia et al. [69], who showed that vegetation phenology and productivity responded differently to heatwave and had contrasting recovery trajectories after a wildfire disturbance. In addition to the direct impacts from a changing climate or climate extremes on vegetation phenology and productivity, an increasingly abundant deciduous shrub cover in tundra landscapes has been reported in high latitudes and arctic regions [72]. Such changes in vegetation cover may also be a driver of divergent trends and changing correlations between LOS and GSI in some high latitude regions that were found in this study.
It has been observed that forest biomes and croplands are more frequently characterized by converging trends than non-forest biomes in the NH ( Figure 4A,B), whereas non-forest biomes are more diverging than croplands and forest biomes ( Figure 4F), suggesting biome-specific temporal changes in the coupled vegetation phenology and productivity. This was also confirmed by several studies at the regional scale [56,69,73]. Wu et al. [56] found that both SOS and EOS had location-dependent impacts on annual GPP. This is further supported by Xia [69] et al. who revealed that the inter-biome variations in annual GPP can be better explained by the variations in seasonal maximal photosynthetic capacity than LOS. In addition to the impact from ongoing global warming and natural climate variability, the relationship between vegetation phenology and productivity could also be affected by changes in both/either vegetation metrics induced by anthropogenic land use/cover changes (e.g., deforestation/afforestation, land clearing, irrigation and fertilization, and changes in land management practices) [53,60,[74][75][76]. Whereas, it could be assumed that changes in land management may also cause divergence between LOS and GSI (e.g., a change from one crop-type to another with different LOS/GSI characteristics), it is noteworthy that both converging trends and pixels with significant positive correlations between LOS and GSI are more pronounced for mixed forest than for the three natural forest biomes, which is also the case for croplands as compared to the three natural non-forest biomes ( Figures 4A and 6B). These differences may reflect anthropogenic management in the form of fertilization and irrigation that could alleviate climatic constraints on vegetation growth influencing concurrently on both phenology and productivity [74,76]. Moreover, it should be noted that areas of diverging trends type 2 in croplands are much smaller than it is the case for the three natural non-forest biomes ( Figure 4F). This further indicates that anthropogenic land management is likely to be able to alleviate adverse climatic effects on vegetation phenology and productivity for both farming, and forestry and thereby weaken the observed diverging trends in LOS and GSI. Contrastingly, it could be expected that changes in land use (e.g., changes in crop types or land abandonment) could induce a differential response on LOS and GSI, as observed in Horion et al. [60], who revealed widespread changes in ecosystem functioning in the northern EA agricultural regions after the collapse of the Soviet Union in 1991 due to widespread agricultural land abandonment. Such extreme case of large scale land abandonment is likely to influence the results of divergent (converging/diverging) trends and changing correlations between LOS and GSI for major land biomes of the NH found here (Figures 3-8). Overall, this suggests the existence of complex relations and interactions between vegetation phenological changes and interannual variability of vegetation productivity induced by the interplay between both natural climate change and land management.

Changing Correlations between LOS and GSI
This study shows widespread changes in the strength of correlations between LOS and GSI (70.90% of the study area) during the past three decades in the NH (Figure 7). Regions that are covered by significant positive correlation between LOS and GSI are becoming more dominant over time, as compared to areas of decreasing correlations. However, by comparing Figures 5 and 7, it can be observed that grasslands in the mid-latitude arid/semi-arid areas with low/negative correlations (e.g., central NA, central Asia, Mongolia, and north China) show significant decreasing strength in the correlation between LOS and GSI, whereas high-latitude tundra regions with low/negative correlations (e.g., arctic area of NA and EA) exhibit a significant increasing correlation. Since both regions are less affected by direct human influence, this difference may reflect latitude-and/or biome-dependent changes in ecosystem functioning in response to long-term climate change or climate extremes/natural disturbance. In contrast to significant changes in these non-forest biomes, both positive/negative trends in correlations between LOS and GSI are rarely found in forest biomes ( Figure 8A). This is also confirmed by Chen et al. [77] who revealed non-significant correlation between LOS and plant productivity in most of southern Canadian Arctic tundra and further suggests biome-specific changes in long-term correlations between LOS and GSI.

Conclusions
Global warming has stimulated vegetation growth through both extending growing season and promoting photosynthesis in the Northern Hemisphere (NH), but the impact is not necessarily uniform across different biomes that are characterized by differences in ecosystem functional types and different degrees of influence from human management. The complementary information about changes in vegetation growing season length and productivity might thereby provide new insights to be used for an improved understanding of the coupled impact from climate change and land management (land use and land cover changes) on changes in vegetation functioning and traits. Here, we analyzed spatio-temporal patterns of observed converging/diverging trends in vegetation phenology and biomass over recent decades (1982-2013) as a function of major biomes covering the NH by using GIMMS3g NDVI data and MODIS land cover information.
A quarter area of the NH was found to be covered by converging trends in length of vegetation growing season (LOS) and growing season NDVI integral (GSI), while diverging trends covered about 6% regions. The results showed spatially-distinct and biome-specific patterns between the continental land masses of Eurasia and North America. Areas of diverging trends were mainly observed in high latitudes and arid/semi-arid areas, whereas forest biomes and croplands showed more widespread signs of converging trends than non-forest biomes. Whereas, the areas of diverging trends in LOS and GSI in high latitude biomes is likely to be a sign of warming induced changes in vegetation species the converging trends for cropland areas may reflect anthropogenic management in the form of fertilization and irrigation that could alleviate climatic constraints on vegetation growth influencing concurrently on both phenology and productivity. The temporal changes in the coupled vegetation phenology and productivity suggest complex relationships and interactions between vegetation phenological changes and the interannual variability of vegetation productivity and future studies should emphasize detailed analyses of hot-spot areas of diverging trends to improve our understanding of the forcing mechanisms of the observed changes in terrestrial ecosystem functioning.