Relative Contribution of Growing Season Length and Amplitude to Long-Term Trend and Interannual Variability of Vegetation Productivity over Northeast China

Abstract

In the context of global warming, the terrestrial ecosystem productivity over the Northern Hemisphere presents a substantially enhanced trend. The magnitude of summer vegetation maximum growth, known as peak growth, remains only partially understood for its role in regulating changes in vegetation productivity. This study aimed to estimate the spatiotemporal dynamics of the length of growing season (LOS) and maximum growth magnitude (MAG) over Northeast China (NEC) using a long-term satellite record of normalized difference vegetation index (NDVI) for the period 1982–2015, and quantifying their relative contribution to the long-term trend and inter-annual variability (IAV) of vegetation productivity. Firstly, the key phenological metrics, including MAG and start and end of growing season (SOS, EOS), were derived. Secondly, growing season vegetation productivity, measured as the Summary of Vegetation Index (VIsum), was obtained by cumulating NDVI values. Thirdly, the relative impacts of LOS and MAG on the trend and IAV in VIsum were explored using the relative importance (RI) method at pixel and vegetation cover type level. For the entire NEC, LOS, and MAG exhibited a slightly decreasing trend and a weak increasing trend, respectively, thus resulting in an insignificant change in VIsum. The temporal phases of VIsum presented a consistent pace with LOS, but changed asynchronously with MAG. There was an underlying cycle of about 10 years in the changes of LOS, MAG, and VIsum. At a regional scale, VIsum tended to maintain a rising trend in the northern coniferous forest and grassland in western and southern NEC. The spatial distribution of the temporal trends of LOS and MAG generally show a contrasting pattern, in which LOS duration is expected to shorten (negative trend) in the central cropland and in some southwestern grasslands (81.5% of the vegetated area), while MAG would increase (positive trend) in croplands, southern grasslands, and northern coniferous forests (16.5%). The correlation index for the entire NEC suggested that LOS was negatively associated with MAG, indicating that the extended vegetation growth duration would result in a lower growth peak and vice versa. Across the various vegetation types, LOS was a substantial factor in controlling both the trend and IAV of VIsum (RI = 75%). There was an opposite spatial pattern in the relative contribution of LOS and MAG to VIsum, where LOS dominated in the northern coniferous forests and in the eastern broadleaf forests, with MAG mainly impacting croplands and the western grasslands (RI = 27%). Although LOS was still the key factor controlling the trend and IAV of VIsum during the study period, this situation may change in the case peak growth amplitude gradually increases in the future.

1. Introduction

Global change has led to significant changes in the trends and interannual variability of terrestrial ecosystem productivity [1,2]. Exploring the specific mechanisms driving these changes in vegetation gross primary productivity is needed for more accurately predicting how the carbon cycle will adjust in future climate scenarios [3,4]. In the Northern Hemisphere, climate warming affects vegetation productivity in two main ways. One is the direct increase in photosynthetic activity and plant growth potential caused by rising air temperature, and the other is the indirect effect of changes in plant phenology, soil moisture deficit, increase in soil nitrogen mineralization, etc. [5]. Global warming has changed vegetation phenology significantly with complex spatiotemporal patterns. Understanding the impacts of phenology on vegetation productivity remains limited both at global and regional scales. Meanwhile, the controlling effects of environmental factors (temperature, precipitation, and solar radiation) on IAV (inter-annual variability) of vegetation productivity have been widely reported for various geographical extents [6,7,8], but the biological mechanisms underlying the changes in carbon cycle are far from clear.

Vegetation phenology refers to the cyclical life rhythm of plants, such as flowering and leaf unfold, which directly reflect vegetation growth and the effects of climate change [6]. Prior studies have shown that phenological changes have a more significant effect on vegetation productivity than warming [5,9]. Currently, a prolonged the growing season, which is primarily caused by an earlier spring onset and a delayed autumn senescence, is widely reported as the main factor driving the enhanced carbon uptake across various biomes [1,5,10]. The timing shifts in phenological events substantially reshape ecosystem function and structure [11,12]. In addition, the peak growth of vegetation, defined as the maximum growth amplitude in vegetation greenness, productivity or leaf area, is also a pivotal physiological factor affecting the seasonal characteristics of terrestrial ecosystem productivity and atmospheric CO₂ concentration [13,14]. Earlier vegetation greenup in spring may induce water and nutrient stress to vegetation peak growth during summer, resulting in a decline of vegetation productivity [15,16]. A recent study has suggested that peak vegetation growth plays a more critical role than climatic factors in regulating the IAV of terrestrial ecosystem carbon sequestration [13]. Thus, growing season length and peak vegetation growth are the key ecological factors influencing vegetation productivity. In the case of the ongoing greening of the globe [17], how the two factors contribute to carbon uptake remains region-specific.

Northeast China (NEC), located in the northern mid-latitudes, is a typical transect zone in the International Geosphere Biosphere Program (IGBP) [18]. It is a sensitive region with a typical temperate continental climate and a hotspot in global change studies [17,19]. Plants in NEC have distinct seasonal characteristics and therefore they are suitable for monitoring phenophase changes using satellite remote sensing. Recent studies on phenological changes in NEC mainly focused on changing trends in the timing of phenological events, in which advanced spring onset is widely documented but with different spatial patterns depending on vegetation type [20,21]. The driving factors of the phenological shifts are primarily environmental, such as temperature, precipitation and soil moisture [22,23]. Net primary productivity (NPP) of croplands and grasslands in NEC has been found to be slightly and negatively correlated with growing season length, but positively associated with the spring greenup date [24]. The combined effects of growing season length and peak growth on controlling vegetation productivity in NEC have not been addressed in previous studies. Exploring changes in vegetation phenology and the impact on IAV productivity in NEC can contribute to understanding carbon sink processes in mid-high latitude terrestrial ecosystems.

This study aimed to investigate the multi-year dynamic and controlling effect of growing season length and peak growth on vegetation productivity over NEC. The most current long-term satellite vegetation greenness data (i.e., Global Inventory Monitoring and Modeling Studies (GIMMS) group Normalized Difference Vegetation Index dataset, referred as NDVI3g, spanning the period 1982–2015) was used to derive key phenological metrics, including the start date, end date, length and maximum magnitude of growing season (SOS, EOS, LOS, MAG) and their long-term trends were analyzed. Using the relative importance method, the contributions of LOS and MAG to growing season production (accumulated NDVI) were quantitatively investigated and the underlying mechanism and uncertainties were discussed.

2. Datasets and Methods

2.1. Study Area and Vegetation Greenness Data

The study area focused on NEC (Northeast China) that includes four provinces: Heilongjiang, Jilin, Liaoning and part of Inner Mongolia (Figure 1). The entire NEC covers an area of 144.8 × 10⁴ km², accounting for about 15% of the area of China. Spatial range of NEC spans 115°E~135°E and 38°~56°N. There is a varied topography over NEC that is primarily composed of mountains, hills, and plains (Figure 1a) and the corresponding altitude ranges from 0 to 2667 m. The two mountainous areas, Daxing’anling Mountain and Changbai Mountain are mainly located in the northern and eastern parts, respectively. The central and northeast sections are dominated by the Songnen Plain and the Sanjiang Plain, respectively. A large portion of NEC belongs to a temperate zone characterized by a continental monsoon climate with four distinct seasons, especially including the cold winters and warm summers. Annual rainfall over NEC ranges from 400 to 800 mm, with a spatial pattern decreasing from east to west [25]. The relative humidity presents a similar spatial pattern. NEC has a high proportion of vegetation coverage, composed of temperate forests, grasslands, and agriculture (Figure 1b). The spatial distribution of NEC vegetation cover types is obtained from the vegetation map of the People’s Republic of China (1:1,000,000) [26]. The coniferous forest is mainly distributed in Daxing’anling Mountain area, while broad-leaved forests are spread in Changbai Mountain area. The central and eastern plains are farmlands covered by cultivated plants. Grassland is primarily distributed in the northwestern part of NEC.

NDVI is a commonly used vegetation greenness that is tightly related to vegetation photosynthetic activity [27,28] and vegetation gross primary productivity [29,30]. Here vegetation greenness from NDVI3g V1.0 was applied as a proxy of vegetation productivity [28]. This dataset has a spatial resolution of 1/12 degree (approximately 8 km along the equator) with a 15-day temporal frequency from 1982 to 2015. The NDVI3g V1 product has been widely used since being published and is the most consistent long-term satellite vegetation dataset currently available [31,32]. In NDVI3g V1.0, effects of satellite orbital drifts, inter-sensor calibration and aerosols impacts from volcanic eruption have been corrected [16,33,34].

The multi-year satellite NDVI images were stacked in chronological order and the time series was smoothed with a Savitzky–Golay filter [35,36]. The phenological parameters were extracted from the smoothed time series and their statistics for diverse vegetation types were calculated. As a comparing verification, EVI (Enhanced Vegetation Index) dataset from MODIS (MOD13C1) was utilized to derive the same phenological metrics over NEC based on the same method of NDVI. Spatiotemporal patterns for the phenological metrics were visually compared.

2.2. Determination of Phenological Metrics

To quantitatively analyze the influence and contribution of vegetation phenological parameters to summer productivity, several key phenological metrics were derived: start date (SOS) and end date (EOS) of growing season, summer peak NDVI (PEAK), and length of growing season (LOS), determined as the difference between EOS and SOS. These phenological indices provide key information for understanding vegetation growth process and carbon cycles. SOS and EOS were identified as the point of curvature change of seasonal growth curve [37,38] (Figure 2). The phenological parameter extraction formula is:

$$NDV{I}_{\left(t\right)}=\left(m_{NDVI}-w_{NDVI}\right)\left(\frac{1}{1+{\mathrm{e}}^{-mS\left(t-S\right)}}+\frac{1}{1+{\mathrm{e}}^{-mA\left(t-A\right)}}-1\right)+w_{NDVI}$$

(1)

where NDVI(t) is the NDVI value at time t, and w_NDVI and m_NDVI denote winter NDVI and the maximum NDVI, respectively. S and A represent the inflection points on the growth curve for SOS and EOS, respectively. mS and mA are the increasing and decreasing rates at the inflection points. w_NDVI is derived with the following equation:

$$w_{NDVI}=\sqrt{\max \left(NDV{I}_{10}\right)\times \max \left(NDV{I}_{11}\right)}$$

(2)

where NDVI₁₀ and NDVI₁₁ are the NDVI values for October and November. With the exception of W_NDVI, other parameters were obtained by iteratively applying the nonlinear least squares method [37].

Based on the vegetation growth curve and the phenological events points, vegetation production during the growing season was approximated using accumulated NDVI. First, the multi-year averaged growing season NDVI was defined as the baseline, and the growing amplitude (MAG) was set as the difference between PEAK and the baseline value. Next, the growing season production (VIsum) was obtained by integrating NDVI under the growth curve and above the baseline (hatched section in Figure 2). According to previous studies, the vegetation growth process can be decomposed into growing season length and amplitude to reflect the carbon assimilation capacity in the vegetation carbon cycle [39]. The baseline values were assumed to exclude vegetation production during non-growth season and reduce the effects of plant respiration on production as well as other potential forces, such as noise signals or sparse vegetation. To better understand plant growth process over NEC, vegetation growing curves for various vegetation types are delineated in Figure A1.

2.3. Data Standardization and Definition of Variability

The three key phenological parameters, LOS, MAG, and VIsum, have different units (days, NDVI value and cumulative NDVI value, respectively) and magnitudes. To reduce the influence of large outliers and facilitate the correlation analysis between the three variables, we used the z-score standardization method (Equation (3)):

$$X^{*}=\frac{x-\mu }{\sigma }$$

(3)

where μ is the mean of all samples (x), σ is the standard deviation of all samples, and X* is the normalized value of the phenological parameter.

According to previous studies, the inter-annual variability of vegetation production is defined as the information in the time series of the data after detrending [13,39]. The specific process consists of two steps: 1) performing a least squares linear fitting on the time series data to determine the trend in the data, 2) removing the trend information from the original data, and obtaining the inter-annual variability.

2.4. Relative Importance Analysis

The relative importance (RI) approach is used to measure the relative contribution of vegetation growth amplitude and length to vegetation productivity. Specifically, a multiple linear regression analysis was first performed between the independent variables (LOS, MAG) and the dependent variable (VIsum). Then, the relative contribution of LOS and MAG to VIsum was determined in percentage [40]. The detailed regression model can be expressed as follows:

$$y_i=k_0+k_1{x}_{i1}+k_2{x}_{i2}+{\epsilon }_i$$

(4)

In the formula, y_i denotes annual VIsum, x_{ip, p = 1,2} represents LOS or MAG in each year, k_{ip, p = 1,2} is the linear regression coefficient, and ε_i is the residual term accounting for the contribution of other factors to vegetation production. The regression coefficient k_{ip, p = 1,2} was determined by minimizing the sum of the squares of ε_i. If ${\widehat{k}}_p$, ${\widehat{y}}_i$ represent the regression coefficient and the result obtained by the multiple linear regression model, respectively, then the determined coefficient R² can be expressed as:

$$R^2=\frac{{{\displaystyle \sum }}_{i=1}^n{\left(\widehat{y}-\overline{y}\right)}^2}{{{\displaystyle \sum }}_{i=1}^n{\left(y_i-\overline{y}\right)}^2}$$

(5)

The relative importance is then given as follows:

$$RI\left(x_p\right)=\frac{1}{2}{\displaystyle \sum }_{j=0}^1\left({\displaystyle \sum }_{i=1}^1\frac{R^2\left(M\right)-R^2\left(S\right)}{\left(\begin{array}{c}1\\ i\end{array}\right)}\right)$$

(6.1)

$$\begin{array}{c}S\subseteq \left\{x_1,x_2\right\}\backslash \left\{x_p\right\}\\ M\subseteq \left\{x_1,x_2\right\}\\ n\left(S\right)=j\end{array}$$

(6.2)

where x_{ip, p = 1,2} represents MAG or LOS, $RI\left(x_p\right)$ is the relative contribution of each variable, M denotes the samples of variables related to MAG and LOS in the regression model, S is the subset of samples, and R² (M) and R² (S) are the coefficients of determination for M and S [40].

2.5. Trend Analysis

Here, the Theil-Sen trend detection method was used to quantify the trends in vegetation production and in key phenological metrics in NEC. Unlike the trend estimation with a simple linear regression, the Theil-Sen trend analysis is a nonparametric regression method and is insensitive to significant outliers in the samples. The most important reason for using the Theil-Sen method in this study is that it does not require the assumption of normal distribution for time series data. The Theil-Sen method has been extensively used for detecting trends of long-term variables in meteorological, hydrological, and ecological datasets [41]. When dealing with temporal data with a skewed distribution or high heteroscedasticity, the accuracy of Theil-Sen trend estimation is much higher than that of the simple linear regression approach, which is not robust enough.

The steps for detecting the Theil–Sen trend of n sample points are as follows:

• The observations are arranged in chronological ascending order, and the Sen slope (Q_k) is calculated using the following formula:

  $$Q_k=\frac{x_j-x_i}{j-i}$$

  (7)

  where x_i and x_j represent the data values corresponding to the time points t_i and t_j (i > j), k = 1,…,N.
• The size of the Sen vector matrix is N = (n(n − 1))/2, where n is the number of time intervals. All Q_k are arranged in ascending order, and then the Sen slope is obtained by calculating the median (Equation (8)). The direction of the time series data is determined by the parameter Q_med. A positive value indicates an increasing trend, while a negative value implies a decreasing trend, and the values reflect the trend magnitude:

  $$Q_{med}=\{\begin{array}{c}{Q}_{\left[N+1\right]/2}, N \mathrm{is} \mathrm{odd}\\ \frac{Q_{\left[N\right]/2}+Q_{\left[N+1\right]/2}}{2}, N \mathrm{is} \mathrm{even}\end{array}$$

  (8)
• The statistical significance of the Theil-Sen trend is tested using the Mann-Kendall method. The Mann-Kendall method is a robust nonparametric statistical test method, which can reduce the influence of outliers in time series data [42]. In this study, the significance level was set to 0.05.

3. Results

3.1. Inter-Annual Varibilities in Vegetation Productivity and Phenological Metrics

During the past thirty-four years, the overall vegetation productivity (or vegetation greenness) in NEC did not show a significant trend. LOS slightly shortened during the past three decades (declining trend), while MAG presented a relatively strong increasing trend (Figure 3). The inter-annual variability in detrended VIsum, LOS and MAG (Figure 3b) was generally in agreement with the trend in the time series (Figure 3a). There was a potential periodicity in the variability of LOS, MAG and VIsum of approximately ten years (such as in the time range 1985–1995). To accurately quantify the period of the cycle, we determined the dominant frequency of VIsum, LOS, and MAG from a spectral analysis of the time series using the R package ‘forefast’ [43]. This procedure revealed that the dominant frequencies of LOS and VIsum were eleven and ten years, respectively. MAG had a relatively large inter-annual fluctuation and did not exhibit a remarkable period, but visually a shorter period of around 8–12 years could be discerned (but not identified by spectral analysis). The directions of the amplitudes in LOS (VIsum) and MAG showed an opposite pattern, in other words, MAG decreased with increasing LOS. Generally, LOS kept a consistent phase with the changes in VIsum.

Further, in different vegetation types, there were significant differences in spatially-averaged trends for the three variables. The increasing rate of MAG and VIsum in croplands was the highest, with a Sen slope of 0.17 and 0.11, while LOS significantly decreased (Sen slope = −0.15). MAG, LOS and VIsum in needleleaf forests and grasslands tended to increase over time (the forest Sen slope was 0.15, 0.28, 0.09, and grasslands Sen slope was 0.15, 0.06, 0.10). The multi-year averaged MAG, LOS and VIsum in the eastern broad-leaved forests were higher than in other vegetation types, but the rising rate of MAG and VIsum were small (Sen slope was 0.11 and 0.06), while LOS presented a decreasing tendency (Sen slope = −0.08).

At pixel level, the long-term trends for MAG, LOS, and VIsum were estimated using Theil-Sen slope. The spatial distribution of these trends showed a widespread diversity in different vegetation areas (Figure 4). MAG exhibited a significant positive trend in the central cropland area as well as a weak positive tendency in the northern and eastern forests (about 81.5% of the vegetated areas showed positive trends). However, in the southwestern grassland and the southeastern forest, MAG had a negative trend. LOS with positive trends was mainly distributed in the northern coniferous forest on Daxing’anling Mountains and in a small patch of broad-leaved forest in the southern portion (around 16.5% of the vegetated area showed positive trends). A vast cluster of croplands in the central part showed a negative trend in LOS. There were fewer significant areas of change in summer productivity VIsum than of MAG and LOS. VIsum with positive trends were mainly distributed in the central and southern croplands (about 49% of the vegetated area showed positive trends). VIsum in southwestern grasslands (around 45°N) and eastern forests had a decreasing tendency. In comparison, the spatial segmentation of long-term trends in MAG and VIsum presented a consistent pattern as both of them exhibited positive trends in the central croplands, and negative trends in southwestern grassland and eastern forest. This phenomenon may imply that vegetation growth amplitude plays an essential role in positively affecting vegetation production.

During the past three decades, the Theil-Sen trend of SOS and EOS over the entire NEC was 0.122 days/yr and −0.137 days/yr, respectively, indicating that spring onset was delayed and autumn senescence came earlier. At pixel scale, SOS trends were dominated by positive trends (72% pixels), prominently emerging in central and eastern cultivated vegetation and eastern forest areas (Figure 5a). On the other hand, negative SOS was discretely distributed in the Daxing’anling coniferous forest area, western grassland and a small patch in the southern forest area. In contrast, the majority of EOS had negative trends distributed in central croplands and western grasslands. Positive areas only accounted for 18.3% of the total EOS trends, which significantly occurred in the Daxing’anling coniferous forests, central meadows and in a small patch in the southern forests. The SOS and EOS delayed (positive trend) regions accounted for 72% and 18.3% of the vegetated area, respectively. Overall, the green-up dates for NEC vegetation has been delayed, and the EOS period has arrived earlier.

3.2. Relative Contribution of SOS and EOS to Vegetation Productivity

As mentioned above, SOS and EOS play a key role in controlling the seasonal growth cycle of vegetation, and therefore it is important to investigate their impacts on growing season productivity. First, the temporal correlation of SOS and EOS with VIsum over the thirty-four years study period is determined using the Spearman correlation coefficient at pixel scale. Then the spatial characteristics of these correlation coefficients are represented with two raster layers. To compare the magnitude of the relationships between SOS and EOS with VIsum, the absolute values of the two correlation coefficient layers are subtracted (|r_EOS| minus |r_SOS|) to derive the spatial distribution of the differences in the contribution of SOS and EOS to vegetation productivity (Figure 6). According to the sign of the subtraction, vegetation productivity of the northern coniferous forest is more correlated with EOS (positive values), while a vast portion of croplands in NEC have a significant link with SOS (negative values). There are no significant differences between the impacts of SOS and EOS on vegetation productivity in the eastern broad-leaved forest and in the coniferous and broad-leaved mixed forest. In parts of the steppe grassland in the East, SOS is pronouncedly correlated to vegetation productivity. The differences in the effects of the phenological indicators on productivity over different vegetation types can be explained by plant physiological processes and by the trends in the phenological indicators. Herbaceous crops are usually characterized by slow growth and rapid physiological maturity due to their genetic type and artificial management [44], thus EOS of crops will occur earlier than that of natural vegetation. In addition, EOS of crops in this study area shows a negative trend (i.e., it is advanced), so these factors may result in the weak relationship between EOS and vegetation productivity (negative delta values). In contrast, the EOS in the forests of the study area has a delayed trend, so its controlling effect on vegetation productivity is relatively significant.

3.3. Spatial Characteristics of Relative Importance of MAG and LOS to Vegetation Productivity

In this paper, the relative importance (RI) approach was used to detect the contribution of MAG and LOS to vegetation production. From the spatial pattern of RI, we found that LOS played a more dominant role than MAG in regulating either the long-term trends (Figure 7a,b) or the inter-annual variability (Figure 7c,d) of vegetation productivity, indicating that summer productivity in each vegetation type was still mainly controlled by LOS, not by MAG. The high RI of MAG was found in the western steppe grassland, followed by central croplands with a relatively smaller RI. Daxing’anling coniferous forests and eastern broad-leaved forests had the lowest RI for MAG. The RI of MAG in the steppe in the southwestern part is not significant. In contrast, the RI of LOS was relatively lower in the western steppe area, but very conspicuous in the Daxing’anling coniferous forests and in the eastern broad-leaved forests. Parts of croplands were slightly controlled by LOS. A visual comparison of Figure 7a,b (row one) with Figure 7c,d (row 2), suggests that the RI of LOS and MAG was likely stronger for long-terms trends than for inter-annual variability across the various vegetation types.

For each vegetation type, Figure 8 shows that LOS was still the main factor controlling the long-term trend of vegetation productivity (mean magnitude of RI around 75%). LOS in the needle- and broad-leaved mixed forest (MF) had the largest contribution to VIsum trends. The relative contribution ratio of MAG to VIsum in each vegetation area was approximately 23%. Although LOS in the steppe area still mainly controlled the long-term trend of productivity, the influencing degree of MAG was significantly greater than that of other vegetation types, which is consistent with the spatial distribution (Figure 7).

In addition to the long-term trends, we also calculated the IAV of LOS, MAG and VIsum for different vegetation types. The effects of IAV of LOS and MAG on IAV of regional plant productivity (VIsum) were also examined (Figure 9). IAV of VIsum in each region had a significant positive correlation with IAV of LOS, but had an insignificant negative correlation with MAG. In Figure 9, except for grasslands, all other linear regressions were statistically significant (p < 0.01). Among them, the IAV of LOS in the evergreen needleleaf forest (ENF) had the greatest impact on IAV of VIsum (R² = 0.85), followed in order by needle-leaved, broad-leaved mixed forests (MF) and deciduous broad-leaved forests (DBF). The smallest contribution of IAV of LOS occurred in the steppe area (53%) (Figure 9m). The IAV correlation between MAG and VIsum was most pronounced in DBF forest (R² = 0.39), followed by meadows, ENF and croplands. IAV of VIsum in grass exhibited the weakest relationship with IAV of MAG. Vegetation growth in croplands is generally intensively regulated by anthropogenic activities (e.g., irrigation, fertilization), and LOS had a greater impact on the variability of productivity (33%). Correlation analysis indicated that the variability of LOS was still the main factor controlling the inter-annual variability of productivity. The impacts of LOS and MAG in forests on productivity variability were the most significant, and the control on productivity variability was the weakest in grassland areas.

4. Discussion

4.1. Plant Phenological Changes over NEC

In the Northern Hemisphere, the impacts of warming-related shifts in seasonality on ecosystem productivity have been extensively investigated [2,16]. Generally, vegetation growing season has progressively lengthened due to both earlier spring and delayed autumn [8] and consequently the photosynthesis period for land vegetation has increased, leading to an increase in vegetation productivity. At the same time, widespread browning and reversal of or stalled greening have been reported at high latitudes [45,46,47]. However, the adaptation of vegetation to climate change in different regions is distinct and limited, so vegetation phenology is location-specific and shows diverse patterns across various biomes. This study identified prominently delayed spring onset over 73% of the vegetated areas and advanced autumn senescence across 82% of the vegetated areas of NEC, which is the opposite of what has been found by some previous studies of the middle and high latitudes of the Northern Hemisphere [14]. The delayed arrival of SOS and the earlier autumn senescence jointly control the shortening of the overall vegetation growing season over NEC. Recent studies have also suggested that the inherently short spring greenup period of northern vegetation will limit its ability to advance further [48,49]. On the other hand, the warming-induced drought stress in autumn will adversely affect the persistence of a delayed end of the growing season [48,49]. The slowing of the SOS advancement since the early 2000s has been well documented over some regions in the Northern Hemisphere, such as Siberia and the Tibetan Plateau [50,51]. These findings could support our results of delayed spring onset and earlier autumn offset.

For the shortening of the cropland growing season, our finding is consistent with the results of a former NPP-based study focusing on NEC [24]. However, this study reveals that the summer peak growth of cultivated vegetation shows an increasing trend over NEC and this may be directly attributed to human farming activities. In addition, warming and CO₂ fertilization also indirectly promote peak growth of crops.

According to relative contribution analysis, LOS and MAG jointly control the IAV of VIsum, but present an opposite effect. It is hard for vegetation to maintain a longer LOS and higher MAG concurrently. Thus, the negative temporal correlation between MAG and VIsum may be attributed to the strong positive control of LOS on VIsum. The ten-year periodicity in LOS, MAG and VIsum variability is an interesting finding in this study, which may be attributed to various climatic and gene type factors. As prior study identified, temperature and precipitation have a consistent step with vegetation growth rhythm [52]. From a planet perspective, this periodicity may be caused by the cosmic power. For example, according to our exploratory analysis, these phenological indicators (1982–2015) have a strong correlation with sunspot activity that has been reported an eleven-year periodicity [53,54,55]. This finding highlights the drivers of plant phenological shifts beyond Earth climate.

4.2. Attributing Regulation of Growing Amplitude and Length on Productivity

Across the entire NEC, the insignificant trend in vegetation productivity may be attributed to both shortened LOS and increasing MAG of vegetation greenness. In other words, the shortened duration of the growing season may partly offset the effects of increased greenness amplitude, resulting in the insignificant trend in total vegetation productivity. From the assessment of the relative contribution over different biomes, the impacts of LOS and MAG on the long-term trend and IAV of vegetation productivity are more significant in forests than in various grasslands (Figure 7 and Figure 9). The underlying mechanism may be the much higher productivity and longer growing season of forests compared to grasslands and croplands. However, the contributions of increasing MAG to productivity in croplands and grasslands (mainly steppe) are relatively more pronounced than in forests (Figure 8). Although in this study in both forests and grasslands growing season length has a significant advantage over peak magnitude in controlling IAV of vegetation productivity, this may be altered in the future due to climatic changes and vegetation growth ability. Previous studies have also shown that forest ecosystems are major contributors to land surface productivity, but the maximum growing season amplitude of other vegetation types is likely to continue to increase. Photosynthesic capacity of leaves varies greatly among different vegetation functional types. Studies under the warming condition also showed that the relationship between vegetation growth activity and temperature variability was weakened, and the growing season lengthening was also slowdown [22,56]. Therefore, if the MAG continues to increase, the main controlling effect of the MAG on the carbon sequestration capacity of the vegetation may occur in the future. Recent study also demonstrated that the maximum LAI during the peak growing season accounts for one half of the annual mean LAI, suggesting that growing amplitude played a key in the carbon cycle [12].

The effect of growing season amplitude on productivity trends and variability is greater in herbaceous vegetation areas such as crops and grasslands than in woody vegetation areas such as forests. This is not only related to the ecological functioning of vegetation, but also the impact of human management activities. Recent studies have shown that China and India dominate the increase in global leaf area, mainly due to the contribution of cropland areas [57]. The low herbaceous vegetation has a short growth period, which also highlights the effect of summer growing peaks on productivity development.

4.3. Validation from MODIS EVI

MODIS EVI at 500 m spatial resolution was analyzed to validate the trends and interannual variabilities of VIsum, LOS and MAG (Figure 10). For the entire NEC, all three indices demonstrated a positive trend, in which VIsum and MAG had a significant trend slope. In contrast to the monotonic trend of MODIS EVI-based VIsum, NDVI3g-based VIsum exhibited fluctuations in its temporal variability, which first rised during 2001–2006, then declined until 2012, after which it had an increasing trend. The obvious differences in VIsum trends derived from EVI and NDVI may be attributed to their ability in capturing vegetation greenness. NDVI exhibits a well-known saturation effect over dense vegetated biomes. EVI is enhanced over NDVI through de-coupling of the canopy background signal (through the use of the soil line) and a reduction in atmospheric influences (through the use of the blue band). Thus, cumulative NDVI over the growing season may affect its inter-annual variability and long-term trend. MODIS EVI-based LOS tended to be prolonged slightly, but insignificantly. Although NDVI3g-based LOS showed a slightly shorter trend, it was not statistically significant. In other words, neither of these findings can confirm significant trends for the developing dynamics of LOS. The trend for MODIS EVI-based MAG was in agreement with that derived from NDVI3g, which generally showed an increase during 2001–2015, but with fluctuations. Due to the brevity of the MODIS EVI time series, there no periodic cycle was found in these indices.

The relative importance of LOS and MAG to VIsum was also examined using the time series dataset of MODIS EVI. Due to the weak impact of IAV, the EVI time series was not decomposed into long-term trend and IAV. Similarly to the NDVI3g derived results (Figure 6), the contributions of LOS and MAG to VIsum exhibited a contrasting spatial pattern (Figure 11). LOS contributed significantly to VIsum in a small extent within the eastern broad-leaved forest and northern coniferous forest, but had a slight impact on VIsum in the central cropland and grassland in the west, whereas MAG played a relatively important role in the eastern grassland and central cropland. It is noteworthy that the contributing magnitude of LOS (mean = 0.51) to VIsum was slightly larger than that of MAG (mean = 0.49). The results that LOS mainly affected forests while MAG was more important in croplands and grasslands were generally in agreement with the results detected with NDVI3g, but with fewer significant pixels.

4.4. Uncertainty Analysis

During the past thirty-four years, anthropogenic activities such as land use and land cover change have also impacted the shift in phenological change patterns due to the conversion of natural vegetated land (forestlands or grassland) into agricultural land over NEC. Further, this phenomenon will create uncertainties in the quantification of the relative contribution of LOS and MAG to growing season productivity. Yet, land cover changes information has to be contained and reflected in the NDVI time series. In other words, the phenological events derived from each NDVI sequence using mature algorithms are generally true and reliable, but they may be the phenological attributions belonging to different vegetation types. To address this issue, we analyzed the land cover changes over NEC. From 1980 to 2000, the shifts in land cover types over NEC were mainly attributed to the conversion of forestlands and grasslands into croplands in DaXing’anling and XiaoXing’anling Mountain areas and of grasslands into croplands in the eastern NEC. During 2000-2006, land cover changes mainly emerged in the Songnen Plain (central part of NEC) with the conversion of dry farmland to paddy farmland and in some degree grasslands to croplands. In general, the more pronounced land cover changes consisted in the conversion of grasslands to croplands and shifts in various agricultural [58]. Due to the 8 km spatial resolution of NDVI and land cover changes occurring mainly in the earlier years of the study period in the NEC, land cover changes probably did not play a dominant role in impacting vegetation phenological dynamics. The long-term trends of SOS and EOS exhibited a similar spatiotemporal pattern in grasslands and croplands that may imply that grasses and crops did not have a significant difference in genetically controlled phenology. Additionally, from Figure 7, we can see that LOS still had the much greater dominant impact factor for vegetation productivity among all vegetation types, so the limited land cover changes would not lead to the reversal of roles LOS and MAG playing in controlling vegetation productivity.

5. Summary and Conclusions

In this paper, the long-term satellite-retrieved vegetation greenness dataset was used to quantify the regulating effects of growing duration length and amplitude on long-term trends and inter-annual variability (IAV) of vegetation productivity over Northeast China (NEC). The conclusions drawn are the following:

• During the past three decades, there was no remarkable change in NDVI-derived vegetation productivity (VIsum) over NEC. However, over vast regions, growing season length (LOS) and amplitude (MAG) presented decreasing and increasing trends, respectively. The temporal profiles for the three phenological variables studied here showed an obvious fluctuation with an approximate cycle of ten years. In the case of persistent trends in LOS and MAG, the regulating effect of MAG on vegetation productivity may be enhanced or even reversed as the main driving force in the future.
• Compared to MAG, LOS was the main factor in controlling the long-term trends and variability of vegetation productivity (VIsum) in NEC. There was a clear spatial cluster characteristic of LOS and MAG in influencing VIsum. The impacts of LOS on VIsum was significant in the northern needleleaf forest and eastern broad-leaved forest, but relatively weak in croplands and steppes. The controlling magnitude of MAG on VIsum was generally small over the whole NEC, with a relatively more pronounced contribution on VIsum in the central-southern croplands and in the western steppe grasslands than in forested areas.
• The key vegetation phenological parameters SOS and EOS indicated a significant shift profile over time. Spatially, SOS was delayed in the majority of NEC (about 72% of vegetated area), except in the northern coniferous forest, and EOS was advanced in the central croplands and western steppe grasslands (about 82% of vegetated area). According to the temporal correlation analysis, EOS had a greater impact on vegetation productivity than SOS in coniferous forest area, while SOS played a more pivotal role than EOS in croplands and in some broadleaf forests.

Overall, this study reveals the spatiotemporal shifts in growing season length and maximum amplitude and their contribution to long-term trends and IAV of vegetation productivity. The approach used here could be utilized in various datasets, such as ground carbon flux and near-ground remote sensing data. Although this is a case study over a temperate zone, it provides broad insights for understanding the mechanisms underlying the trend and IAV of land carbon uptake.