Mapping Spatiotemporal Changes in Vegetation Growth Peak and the Response to Climate and Spring Phenology over Northeast China

: Global climate change has led to significant changes in seasonal rhythm events of vegetation growth, such as spring onset and autumn senescence. Spatiotemporal shifts in these vegetation phenological metrics have been widely reported over the globe. Vegetation growth peak represents plant photosynthesis capacity and responds to climate change. At present, spatiotemporal changes in vegetation growth peak characteristics (timing and maximum growth magnitude) and their underlying governing mechanisms remain unclear at regional scales. In this study, the spatiotemporal dynamics of vegetation growth peak in northeast China (NEC) was investigated using long-term NDVI time series. Then, the effects of climatic factors and spring phenology on vegetation growth peak were examined. Finally, the contribution of growth peak to vegetation production variability was estimated. The results of the phenological analysis indicate that the date of vegetation green up in spring and growth peak in summer generally present a delayed trend, while the amplitude of growth peak shows an increasing trend. There is an underlying cycle of 11 years in the vegetation growth peak of the entire study area. Air temperature and precipitation before the growing season have a small impact on vegetation growth peak amplitude both in its spatial extent and magnitude (mainly over grasslands) but have a significant influence on the date of the growth peak in the forests of the northern area. Spring green-up onset has a more significant impact on growth peak than air temperature and precipitation. Although green-up date plays a more pronounced role in controlling the amplitude of the growth peak in forests and grasslands, it also affects the date of growth peak in croplands. The amplitude of the growth peak has a significant effect on the inter-annual variability of vegetation production. The discrepant patterns of growth peak response to climate and phenology reflect the distinct adaptability of the vegetation growth peak to climate change, and result in different carbon sink patterns over the study area. The study of growth peak could improve our understanding of vegetation photosynthesis activity over various land covers and its contribution to carbon uptake.


Introduction
Vegetation phenology is a sensitive indicator of global climate change. The primary vegetation phenological metrics, including dates for green-up onset and senescence, as well as their interactions with climate change, have been widely examined at regional to global scales [1][2][3][4]. Previous studies, based on data from satellite remote sensing and in situ observations, have reported that vegetation in the mid-high latitudes of the northern hemisphere generally shows a trend of earlier green-up onset and delayed senescence but with complex spatiotemporal patterns [5,6]. As a result, a phenologically induced extension of the growing season has been considered a key factor in the long-term variability of vegetation production [7,8]. Climate attribution analyses confirm that changes in vegetation phenology in the northern mid-high latitudes are mainly controlled by the interannual variation of air temperature, but, in arid and semi-arid areas, vegetation phenology is more significantly affected by temporal and spatial patterns of precipitate [9,10]. Additionally, vegetation growth peak, a primary phenological parameter, also plays a key role in regulating land carbon uptake and carbon cycling [11]. The characteristics of vegetation growth peak, including its date and magnitude, remain poorly understood regarding its spatiotemporal dynamics and the corresponding influence on vegetation production [12]. Thereby, it is of great importance to further explore changes in vegetation growth peak at regional and global scales.
Vegetation growth peak can be described in terms of the date of peak growth and the maximum magnitude of growth. The date of peak growth is usually defined as the position of peak (POP) in Julian day units, and refers to the change point where the vegetation growth process shifts from green up to leaf fall. Due to its position in the plant photosynthetic process, POP is closely linked with the transition stage in annual carbon uptake. Vegetation growth peak (PEAK) is a proxy of the maximum photosynthesis capacity and has substantial impact on carbon uptake and atmospheric CO 2 concentration [13]. POP and PEAK generally imply that a plant has the most abundant resources (e.g., water, temperature, photoperiod) for growth [14]. Changes in POP and PEAK will, in turn, affect the energy balance between land surface and atmosphere, and water and material cycling. The change of vegetation growth peak reflects the adaptability of vegetation in obtaining the best growth status with the least cost, under the background of environmental change [15]. It has been reported that POP in mid-high latitudes of North America tends to occur earlier in spring, which leads to an increase of annual gross primary productivity of vegetation [14,16]. However, the effect of POP is very region specific. Vegetation in the arid and semi-arid areas in the northern mid-high latitudes is mainly controlled by water availability, and the advance of vegetation POP may lead to a decline of PEAK and annual total productivity [16,17].
Northeast China (NEC) is located at the mid latitudes of the northern hemisphere, with a typical temperate continental monsoon climate. This area has a rich composition of plant species and is sensitive to climate change. It has become a typical transect zone for global change research in the International Geosphere Biosphere Programme [18,19]. Previous studies in this region mainly focused on the long-term trends of key phenological parameters and climate response, such as the dynamics of spring green-up onset and changes in growth season length. Among these, an earlier start of the growing season and delayed senescence date have been frequently reported and mainly attributed to warming [2,20,21]. Taking northeast China as a case study will improve our understanding of the spatiotemporal variation characteristics of vegetation growth peak and its impact on the carbon cycle in temperate monsoon regions.
Overall, the changing mechanism and characteristics of vegetation growth peak remain unclear and present a complex spatiotemporal pattern at regional scales. This paper aimed to use a long-term normalized difference vegetation index (NDVI) dataset to explore the spatiotemporal dynamics of vegetation growth peak, and sought to respond to two questions: (1) what is the impact of climatic factors and phenology on vegetation growth peak? and (2) how does vegetation growth peak affect vegetation annual production?

Study Area and NDVI Data Processing
Due to the variety in vegetation types and high plant coverage, northeast China is an ideal place for exploring how vegetation responds to climate based on remote sensing technology. It covers an area of 124 × 10 4 km 2 , and is composed of the provinces of Liaoning, Jilin, Heilongjiang, and the eastern part of Inner Mongolia. The majority of this area is located in the warm temperate and temperate zones, with a small area in the north in the cold temperate zone. The annual average temperature ranges from 11.5 to −4 • C [22]. From east to west, it spans across humid, semi-humid, and arid-semi-arid zones, with annual precipitation ranging from 1100 to 250 mm [22]. Vegetation types mainly consist of broadleaf forest, needleleaf forest, mixed forest, meadow, steppe, grassland, shrub, and cropland ( Figure 1) [23]. Vegetation greenness (NDVI) derived from multi-year satellite images demonstrate that croplands in the center and grasslands in the west have a lower vegetation coverage status than forests in northern and eastern parts.

Study Area and NDVI Data Processing
Due to the variety in vegetation types and high plant coverage, northeast China is an ideal place for exploring how vegetation responds to climate based on remote sensing technology. It covers an area of 124 × 10 4 km 2 , and is composed of the provinces of Liaoning, Jilin, Heilongjiang, and the eastern part of Inner Mongolia. The majority of this area is located in the warm temperate and temperate zones, with a small area in the north in the cold temperate zone. The annual average temperature ranges from 11.5 to −4 °C [22]. From east to west, it spans across humid, semi-humid, and arid-semi-arid zones, with annual precipitation ranging from 1100 to 250 mm [22]. Vegetation types mainly consist of broadleaf forest, needleleaf forest, mixed forest, meadow, steppe, grassland, shrub, and cropland ( Figure 1) [23]. Vegetation greenness (NDVI) derived from multi-year satellite images demonstrate that croplands in the center and grasslands in the west have a lower vegetation coverage status than forests in northern and eastern parts.

Vegetation Greenness Dataset
Vegetation status over the study area is measured in terms of NDVI, which is a widely used proxy of vegetation canopy greenness. Here, we applied the long-term NDVI3g V1.0 dataset (https://ecocast.arc.nasa.gov/data/pub/gimms/) derived from AVHRR satellites, released by the group of Global Inventory Monitoring and Modeling System (GIMMS), in the USA. GIMMS NDVI3g V1.0 covers the period 1982-2015, with a spatial and temporal resolution of 1/12 degree and 1/24 a year, respectively [24]. In the literature of biosphere studies, this dataset has been extensively used to study the carbon cycle, vegetation phenology dynamics, and plant response to climate change [9,25,26]. In GIMMS NDVI3g V1.0, data quality issues caused by AVHRR satellites orbital drifts, inter-sensor calibration, and aerosols impacts from volcanic eruption have been corrected [24,27]. Long-term time series of NDVI was reconstructed using the Savitzky-Golay (S-G) filter to reduce the potential noise in the raw data. The moving average window in the S-G method was set to seven points.

Meteorological Dataset
We selected gridded air temperature and precipitation datasets to examine their relationship with vegetation growth peak. The meteorological dataset was derived from surface observation data obtained from the National Meteorological Information Center (NMIC) of the China Meteorological

Vegetation Greenness Dataset
Vegetation status over the study area is measured in terms of NDVI, which is a widely used proxy of vegetation canopy greenness. Here, we applied the long-term NDVI3g V1.0 dataset (https: //ecocast.arc.nasa.gov/data/pub/gimms/) derived from AVHRR satellites, released by the group of Global Inventory Monitoring and Modeling System (GIMMS), in the USA. GIMMS NDVI3g V1.0 covers the period 1982-2015, with a spatial and temporal resolution of 1/12 degree and 1/24 a year, respectively [24]. In the literature of biosphere studies, this dataset has been extensively used to study the carbon cycle, vegetation phenology dynamics, and plant response to climate change [9,25,26]. In GIMMS NDVI3g V1.0, data quality issues caused by AVHRR satellites orbital drifts, inter-sensor calibration, and aerosols impacts from volcanic eruption have been corrected [24,27]. Long-term time series of NDVI was reconstructed using the Savitzky-Golay (S-G) filter to reduce the potential noise in the raw data. The moving average window in the S-G method was set to seven points.

Meteorological Dataset
We selected gridded air temperature and precipitation datasets to examine their relationship with vegetation growth peak. The meteorological dataset was derived from surface observation data obtained from the National Meteorological Information Center (NMIC) of the China Meteorological Administration (http://www.csdata.org/en/p/80/) [28,29]. The ANUSPLIN software was applied to interpolate observed meteorological data to a grid format file with high resolution (1 km), and the results were validated with temperature and precipitation records from AsiaFlux stations (http: //asiaflux.net) [28,29]. These gridded meteorological data were resampled into the same spatial resolution Remote Sens. 2020, 12, 3977 4 of 16 of 8 km of NDVI3g, using the nearest neighbor method in ArcGIS 10.3. Daily air temperature and precipitation were composited to a 16-day dataset using the average value. Averaged air temperature and cumulative precipitation during one month before the position of peak (POP) date were calculated to evaluate their impact on vegetation growth peak. This one-month interval was determined by several regression experiments. We tested the correlation of air temperature and precipitation one, two, and three months before POP with growth peak, and found that one month was the most significant period.

Calculation of Phenological Metrics
According to the single seasonal growth trajectory of vegetation over the study area, the double sigmoid logistics approach was used to fit the growth curve and derive key phenological metrics. This method is based on the assumption that vegetation canopy greenness observed by satellite is a function of chlorophyll content and leaf transmittance [30]. In the early growing season, the photosynthetic capacity of leaves reached the maximum, and the canopy greenness was only a function of the leaf transmittance in the developing stage. In the late growing season, the photosynthetic capacity of leaves gradually decreased until the dormancy stage. To build a model of the vegetation growth process, the sigmoid logistics curve is used for delineating vegetation green up and the senescence phase. The key technique is based on the extraction of the six parameters in Equation (1) [30]: where v (t) is vegetation greenness (NDVI) at time t (day of year), and v min and v amp represent the background greenness and maximum greenness, respectively. m 1 , m 2 , n 1 , and n 2 are fitted parameters. Specifically, m 1 and n 1 are expected to control the phase and slope of the growth curve during the green-up period, while m 2 and n 1 control these during the senescence period. The sigmoid logistics method is suitable for depicting changes in vegetation structure and is not sensitive to noise [30].
In the fitted sigmoid logistics curve, the date for the midpoint (50%) of maximum NDVI is assigned to the start of the growing season (SOS). This SOS point has several unique characteristics as the maximum slope in the growth curve and the peak point in the first derivative. In the actual growth process, the SOS date could appropriately represent the vegetation leaf unfolding status. The midpoint is defined as a percentage, so even in the case the absolute value of vegetation greenness changes, it can maintain stationarity in the time domain. Likewise, the date of the end of the growing season (EOS) is defined by the midpoint in the dormancy logistic curve. The position of the peak (POP) is determined as the date with the maximum NDVI value (PEAK in the growth curve). Based on these key time points, the green-up period and growing season length can be defined ( Figure 2).

Trend and Correlation Analysis
The Mann-Kendall (M-K) test, which is a non-parametric method, was used to detect the changing trends of vegetation POP and PEAK at the pixel scale. The M-K test does not assume that the samples follow a specific statistical distribution, because it is not sensitive to a few outliers. Thus,

Trend and Correlation Analysis
The Mann-Kendall (M-K) test, which is a non-parametric method, was used to detect the changing trends of vegetation POP and PEAK at the pixel scale. The M-K test does not assume that the samples follow a specific statistical distribution, because it is not sensitive to a few outliers. Thus, it has a wide detection range and high skill for trend detection [31]. The tau value in the M-K test is used to depict the direction and magnitude of changes in vegetation POP and PEAK. The process for obtaining tau is as follows: where the sign (Equation (4)) is the symbolic function, n = 34 (for years 1982-2015). A positive (negative) tau indicates an increasing (decreasing) trend in vegetation growth, and the absolute value of tau represents the magnitude of the growth trend.
The relationship between vegetation peak phenology and environmental factors are examined using the Spearman method, which is a non-parametric rank correlation analysis method. For samples X and Y, both with n observations, X and Y are first sorted, and then the Spearman correlation coefficient is calculated according to Equation (6). In Equation (6), ρ is the Spearman correlation coefficient. A greater ρ indicates a stronger relationship between two samples. cov rg x , rg y is the covariance between sorted X (rg x ) and Y (rg y ). σ rgx and σ rgy are the standard deviation for rg x and rg y , respectively: ρ rgx,rgy = cov rg x , rg y /(σ rgx σ rgy ).

Spatiotemporal Pattern of Vegetation Peak Phenology
The interannual variability of POP and PEAK for the entire area indicates an obvious fluctuation with an insignificant increasing trend ( Figure 3). From visual observation, there is a potential cycling period of about 10 years for POP and PEAK (Figure 3a,b), such as the epoch 1982-1992 in the POP curve. We used the spectral analysis method in the R package 'forecast' to quantify this phenomenon and found an 11-year and 10-year cycling period for POP and PEAK (insignificant for PEAK), respectively. The coefficient of variation for the time series of POP and PEAK are 3.5 and 4.3, and they demonstrate a weak correlation (Spearman r < 0.1). The temporal dynamics indicates that POP and PEAK do not consistently grow towards an increasing trend, especially in the context of warming in the northern mid-high latitudes, indicating that climatic changes are heterogeneous over distinct regions and spatial scales.
The date of the start of the growing season (start of season, SOS), a primary time point representing the beginning of vegetation green up, can affect vegetation growth status during the growing season [32]. Thus, the temporal dynamics of SOS should be explored here. The result of the linear regression indicates that SOS has a slightly increasing trend (insignificant) during 1982-2015, which implies a delayed vegetation green-up date (Figure 3c). There is an analogous time series between SOS and PEAK, of which the Spearman correlation r is 0.86. SOS presents an insignificant relationship with POP (Spearman r = 0.18). The SOS time series is identified with a 10-year cycling period, and has a coefficient of variation of 3.8.
Remote Sens. 2020, 12, 3977 6 of 16 respectively. The coefficient of variation for the time series of POP and PEAK are 3.5 and 4.3, and they demonstrate a weak correlation (Spearman r < 0.1). The temporal dynamics indicates that POP and PEAK do not consistently grow towards an increasing trend, especially in the context of warming in the northern mid-high latitudes, indicating that climatic changes are heterogeneous over distinct regions and spatial scales. The date of the start of the growing season (start of season, SOS), a primary time point representing the beginning of vegetation green up, can affect vegetation growth status during the growing season [32]. Thus, the temporal dynamics of SOS should be explored here. The result of the linear regression indicates that SOS has a slightly increasing trend (insignificant) during 1982-2015, which implies a delayed vegetation green-up date (Figure 3c). There is an analogous time series between SOS and PEAK, of which the Spearman correlation r is 0.86. SOS presents an insignificant The MK test at the pixel scale indicates that the magnitude of PEAK in croplands, southern grasslands, and northern needleleaf forests shows a significant enhanced trend (accounting for 81% of study area), while there is a decreasing trend in the southwestern steppe area and southeastern broadleaf forests (Figure 4a). The date of POP is obviously delayed in eastern and northern broadleaf forests (positive MK tau value), accounting for about 69% of the study area, but is advanced in the steppe and croplands (Figure 4b). There is a visible opposite spatial pattern in Figure 4 that could be fit in the logic of the impact of climate change; a place with earlier POP (negative trend) usually corresponds to an enhanced PEAK (positive trend). The spatial pattern for the trend of SOS presents greater spatial heterogeneity with a discrete distribution of SOS with a negative trend. SOS shows a substantial delayed trend in parts of eastern and central croplands and broadleaf forest (90% of NEC) (Figure 4c). In contrast, only some parts in the northern needleleaf forest and western and southern steppe present an SOS advanced trend (negative trend). The area with a significant SOS trend is smaller than that of POP and PEAK. Trends for POP and SOS have a slightly similar spatial pattern indicating that, in these regions, SOS and POP have a consistent developing trend. Thus, an earlier spring onset will commonly come with an advanced peak growth date, and vice versa. substantial delayed trend in parts of eastern and central croplands and broadleaf forest (90% of NEC) (Figure 4c). In contrast, only some parts in the northern needleleaf forest and western and southern steppe present an SOS advanced trend (negative trend). The area with a significant SOS trend is smaller than that of POP and PEAK. Trends for POP and SOS have a slightly similar spatial pattern indicating that, in these regions, SOS and POP have a consistent developing trend. Thus, an earlier spring onset will commonly come with an advanced peak growth date, and vice versa.

Impacts of pre-POP Temperature and Precipitation on Vegetation Peak
According to the lagged effect of climatic factors on vegetation changes, we selected precipitation and temperature one month prior to the POP date to examine their relationship with POP and PEAK. Actually, we tested the relationship of precipitation and temperature one, two, and three months prior to the POP date with POP and PEAK, and found the most significant relationship

Impacts of Pre-POP Temperature and Precipitation on Vegetation Peak
According to the lagged effect of climatic factors on vegetation changes, we selected precipitation and temperature one month prior to the POP date to examine their relationship with POP and PEAK. Actually, we tested the relationship of precipitation and temperature one, two, and three months prior to the POP date with POP and PEAK, and found the most significant relationship one month before the vegetation growth peak date. A previous study also demonstrated the important impacts of pre-season precipitation and temperature on vegetation peak date and magnitude [32]. In Figure 5, we identified that pre-POP precipitation mainly correlates with POP in northern needleleaf forests and in some patches in the southern and eastern parts (Figure 5a). Meanwhile, pre-POP temperature has a strong correlation with POP in most parts of study area but is not significant in the western steppe and central croplands (Figure 5b). Vegetation PEAK is only correlated with pre-POP precipitation and temperature in small parts in western and southeastern steppe but not correlated in most parts (Figure 5c,d). These findings imply that environmental factors play a more important role in controlling the timing of the peak date than the magnitude of the peak. In addition, precipitation and temperature have very small impacts on croplands that are highly managed by human management activities. Even these environmental factors could affect vegetation peak growth that is just working on the phenological date of the peak. correlated with pre-POP precipitation and temperature in small parts in western and southeastern steppe but not correlated in most parts (Figure 5c,d). These findings imply that environmental factors play a more important role in controlling the timing of the peak date than the magnitude of the peak. In addition, precipitation and temperature have very small impacts on croplands that are highly managed by human management activities. Even these environmental factors could affect vegetation peak growth that is just working on the phenological date of the peak.

Controlling Effect of Spring Phenology on Vegetation Peak
Green-up date has been a widely studied phenological metric that closely relates to global changes. For example, in North America's boreal forests, an earlier green-up date in spring can lead to a reduction of the maximum magnitude of vegetation peak growth in summer [32]. Here, we also explored the relationship of SOS with POP and PEAK using Spearman correlation analysis ( Figure  6). There is a strong positive correlation between SOS and POP in various forests, excluding the central croplands and the western steppe (Figure 6a), indicating a synchronous development of SOS and POP in natural vegetation. Due to anthropogenic management practices, the dates of SOS and POP have some uncertainties that caused an insignificant relationship between SOS and POP. Only in a small part of broadleaf forests, croplands, and steppe, SOS is significantly correlated with PEAK

Controlling Effect of Spring Phenology on Vegetation Peak
Green-up date has been a widely studied phenological metric that closely relates to global changes. For example, in North America's boreal forests, an earlier green-up date in spring can lead to a reduction of the maximum magnitude of vegetation peak growth in summer [32]. Here, we also explored the relationship of SOS with POP and PEAK using Spearman correlation analysis ( Figure 6). There is a strong positive correlation between SOS and POP in various forests, excluding the central croplands and the western steppe (Figure 6a), indicating a synchronous development of SOS and POP in natural vegetation. Due to anthropogenic management practices, the dates of SOS and POP have some uncertainties that caused an insignificant relationship between SOS and POP. Only in a small part of broadleaf forests, croplands, and steppe, SOS is significantly correlated with PEAK ( Figure 6b) but is not significant in most parts of the study area. POP date is synchronized with vegetation growth PEAK in northern needleleaf forest and parts of the eastern broadleaf forests (Figure 6c), indicating that vegetation greenness will grow higher in the summer peak in the case of delayed POP. However, POP is not tightly correlated with PEAK in a majority of croplands.
Remote Sens. 2020, 12, x FOR PEER REVIEW 9 of 16 ( Figure 6b) but is not significant in most parts of the study area. POP date is synchronized with vegetation growth PEAK in northern needleleaf forest and parts of the eastern broadleaf forests (Figure 6c), indicating that vegetation greenness will grow higher in the summer peak in the case of delayed POP. However, POP is not tightly correlated with PEAK in a majority of croplands.

Relative Importance of Spring Phenology and Climate Factors on Vegetation Peak
According to previous analyses, climatic factors and spring phenological parameters have various effects on vegetation peak growth (e.g., POP, PEAK). Here, we used the relative importance method to detect the diverse roles that SOS, pre-POP temperature, and pre-POP precipitation play in POP and PEAK (Figure 7). In the relative importance analysis, SOS, pre-POP temperature, and pre-POP precipitation were set as explanation variables and regressed with POP (PEAK). The outputs are the relative contributions of multiple explanatory variables to the explained variable. Results indicate that pre-POP temperature and pre-POP precipitation have small impacts on vegetation PEAK in small parts of western and southwestern steppe (blue pixels in Figure 7a,c). However, SOS has an inverse effect on PEAK in the spatial distribution in which the significant regions cover a majority of the forests (Figure 7e). Besides the spatial extent, the correlated magnitude between SOS and PEAK is substantially larger than between temperature (precipitation) and PEAK. SOS has a greater impact on POP in croplands and steppe, showing an opposite spatial pattern with PEAK (Figure 7e vs. Figure 7f). Pre-POP temperature only affects the POP date in the northern needleleaf forests and southeastern broadleaf forests. These findings imply that under the context of global and local climatic changes, meteorological factors may not play such a pivotal role in affecting the vegetation summer growth peak as spring phenological events do (e.g., SOS). In other words, vegetation genotype may be a key driver controlling the variability of vegetation growth patterns.

Relative Importance of Spring Phenology and Climate Factors on Vegetation Peak
According to previous analyses, climatic factors and spring phenological parameters have various effects on vegetation peak growth (e.g., POP, PEAK). Here, we used the relative importance method to detect the diverse roles that SOS, pre-POP temperature, and pre-POP precipitation play in POP and PEAK (Figure 7). In the relative importance analysis, SOS, pre-POP temperature, and pre-POP precipitation were set as explanation variables and regressed with POP (PEAK). The outputs are the relative contributions of multiple explanatory variables to the explained variable. Results indicate that pre-POP temperature and pre-POP precipitation have small impacts on vegetation PEAK in small parts of western and southwestern steppe (blue pixels in Figure 7a,c). However, SOS has an inverse effect on PEAK in the spatial distribution in which the significant regions cover a majority of the forests (Figure 7e). Besides the spatial extent, the correlated magnitude between SOS and PEAK is substantially larger than between temperature (precipitation) and PEAK. SOS has a greater impact on POP in croplands and steppe, showing an opposite spatial pattern with PEAK (Figure 7e vs. Figure 7f). Pre-POP temperature only affects the POP date in the northern needleleaf forests and southeastern broadleaf forests. These findings imply that under the context of global and local climatic changes, meteorological factors may not play such a pivotal role in affecting the vegetation summer growth peak as spring phenological events do (e.g., SOS). In other words, vegetation genotype may be a key driver controlling the variability of vegetation growth patterns.

Impacts of Growth Peak on Vegetation Production
The date and magnitude of vegetation growth peak are key parameters that can shape the vegetation seasonal growth profile and will further affect the accumulation of production during the growing season. We used the Spearman correlation method to test the link between vegetation growing season production with POP and PEAK. Here, averaged NDVI (mean NDVI during the growing season, MGS) was used as the proxy of vegetation production, according to the strong relationship between NDVI and GPP (gross primary productivity) [33,34]. The results indicate that over the entire region, POP is slightly and positively associated with MGS in the majority parts, with northern needleleaf forest and the northeastern part showing a relatively high correlation (Figure 8a). This finding of a positive relationship may reflect that a delayed POP date may improve the accumulation of vegetation production during the growing season. In the other case, an earlier POP Figure 7. Relative Importance of SOS, temperature, and precipitation to PEAK (left column) and POP (right column) (black dots represent significant pixels) (a-f).

Impacts of Growth Peak on Vegetation Production
The date and magnitude of vegetation growth peak are key parameters that can shape the vegetation seasonal growth profile and will further affect the accumulation of production during the growing season. We used the Spearman correlation method to test the link between vegetation growing season production with POP and PEAK. Here, averaged NDVI (mean NDVI during the growing season, MGS) was used as the proxy of vegetation production, according to the strong relationship between NDVI and GPP (gross primary productivity) [33,34]. The results indicate that over the entire region, POP is slightly and positively associated with MGS in the majority parts, with northern needleleaf forest and the northeastern part showing a relatively high correlation (Figure 8a). This finding of a positive relationship may reflect that a delayed POP date may improve the accumulation of vegetation production during the growing season. In the other case, an earlier POP date may reduce the accumulation of production. In contrast, in most parts, the variability of MGS can be explained by PEAK except in the croplands of the central area and some steppe areas in the southwestern part (Figure 8b). The tight relationship between PEAK and MGS demonstrates that high PEAK in summer vegetation growth can enhance vegetation production. Cropland growth is highly affected by artificial management practices, including irrigation, fertilization, pesticides application, and seed selection, thus its PEAK growth varies from year to year and further weakens the impact of phenology on cropland production.
Remote Sens. 2020, 12, x FOR PEER REVIEW 11 of 16 date may reduce the accumulation of production. In contrast, in most parts, the variability of MGS can be explained by PEAK except in the croplands of the central area and some steppe areas in the southwestern part (Figure 8b). The tight relationship between PEAK and MGS demonstrates that high PEAK in summer vegetation growth can enhance vegetation production. Cropland growth is highly affected by artificial management practices, including irrigation, fertilization, pesticides application, and seed selection, thus its PEAK growth varies from year to year and further weakens the impact of phenology on cropland production.

Mechanism Analysis of Impact Factors
According to the weak relationship and small spatial distribution, we can confirm that temperature and precipitation have a relatively small influence on the peak amplitude of vegetation. On the other hand, SOS had a significant correlation with the date and amplitude of growth peak. This phenomenon indicates that climatic factors may not be the decisive factors in controlling vegetation growth peak changes. Vegetation physiological factors, such as vegetation genotype, may play a key role in regulating vegetation growth peak [35][36][37]. Temperature and precipitation before the growing season can generally result in a delaying of POP, because good hydrothermal conditions can promote the sustainable growth of vegetation during the green-up period. A previous study also confirms the positive effect of pre-season temperature and precipitation on POP [13].
In this study, SOS over the study area is characterized generally by a delayed trend, with a small part of advanced SOS in the western grasslands and northern needle forests. This finding is not consistent with some previous studies reporting substantial advancements of spring SOS [20,21]. Large-scale studies focusing on northern mid-high latitudes also discovered the complex trend and spatial difference in the SOS trend in various biomes, including delayed trend, earlier trend, and no trend [14,38]. Although the northern hemisphere is generally undergoing a change in greening patterns, there are a considerable number of regions that are browning. As a case study at a regional scale, this study uncovers the complexity of SOS dynamics.

Inferring Vegetation Dynamic in the Future
Under the driving forces from global climate changes and local environmental shifts, the seasonal growth pattern and maximum growth potential for northeast China vegetation may adjust to adapt to these changes. On the basis of the long-term trend and interannual variability of POP and PEAK, we inferred the future possible profiles for vegetation growth peak and key phenological metrics (Figure 9. There are several findings supporting the potential growth curves shown in Figure   Figure 8. Impacts of POP and PEAK on vegetation productivity (correlation coefficient) (a,b).

Mechanism Analysis of Impact Factors
According to the weak relationship and small spatial distribution, we can confirm that temperature and precipitation have a relatively small influence on the peak amplitude of vegetation. On the other hand, SOS had a significant correlation with the date and amplitude of growth peak. This phenomenon indicates that climatic factors may not be the decisive factors in controlling vegetation growth peak changes. Vegetation physiological factors, such as vegetation genotype, may play a key role in regulating vegetation growth peak [35][36][37]. Temperature and precipitation before the growing season can generally result in a delaying of POP, because good hydrothermal conditions can promote the sustainable growth of vegetation during the green-up period. A previous study also confirms the positive effect of pre-season temperature and precipitation on POP [13].
In this study, SOS over the study area is characterized generally by a delayed trend, with a small part of advanced SOS in the western grasslands and northern needle forests. This finding is not consistent with some previous studies reporting substantial advancements of spring SOS [20,21]. Large-scale studies focusing on northern mid-high latitudes also discovered the complex trend and spatial difference in the SOS trend in various biomes, including delayed trend, earlier trend, and no trend [14,38]. Although the northern hemisphere is generally undergoing a change in greening patterns, there are a considerable number of regions that are browning. As a case study at a regional scale, this study uncovers the complexity of SOS dynamics.

Inferring Vegetation Dynamic in the Future
Under the driving forces from global climate changes and local environmental shifts, the seasonal growth pattern and maximum growth potential for northeast China vegetation may adjust to adapt to these changes. On the basis of the long-term trend and interannual variability of POP and PEAK, we inferred the future possible profiles for vegetation growth peak and key phenological metrics (Figure 9. There are several findings supporting the potential growth curves shown in Figure 9. Firstly, POP and SOS in northern and eastern vegetated regions show a delayed trend, but in the western parts, they show an advancing trend. Secondly, SOS and POP generally maintain a synchronously developing trend, and PEAK over the whole region shows a rising trend. In this case, future vegetation growth patterns may follow type 1 and 2 in Figure 9. Growth type 1 may occur in northern and eastern forests, while type 2 may occur in the western steppe and southwestern grasslands. Compared with the current growth profile, vegetation in the future may have more production during the growing season that is caused by a higher growth peak.
Remote Sens. 2020, 12, x FOR PEER REVIEW 12 of 16 9. Firstly, POP and SOS in northern and eastern vegetated regions show a delayed trend, but in the western parts, they show an advancing trend. Secondly, SOS and POP generally maintain a synchronously developing trend, and PEAK over the whole region shows a rising trend. In this case, future vegetation growth patterns may follow type 1 and 2 in Figure 9. Growth type 1 may occur in northern and eastern forests, while type 2 may occur in the western steppe and southwestern grasslands. Compared with the current growth profile, vegetation in the future may have more production during the growing season that is caused by a higher growth peak.

Uncertainty Analysis
In the process of calculating long-term vegetation phenology, land use/land cover (LULC) change must be mentioned due to its impact on the accuracy of phenological analysis [39]. Land cover change, including human-induced and natural forcing changes, will introduce heterogeneity in the time series of vegetation greenness (e.g., NDVI) or structure [1]. Additionally, current phenological metrics extraction methods cannot identify information regarding vegetation type. Thus, the multiyear time series of NDVI may represent mixes of different ecosystem types, such as forest and grassland in different periods. During the period 1982-2015, there were some changes of LULC in the study area. In the period 1980-2000, these LULC changes were primarily changes of forest and grassland to pasture in the eastern area, the change of farming-pastoral land to farmland, and the change of forest to grassland in the northern area. After the year 2000, there were mainly land cover changes between different farmland planting systems, such as the transition from non-irrigated farmland to paddy field in the Northeast Plain, with a small amount of grassland converted to farmland in the central area [40].
These LULC changes not only result in the complexity of vegetation index time series but also further induce uncertainties in computing vegetation phenological metrics. However, the results in this paper can be considered reliable. First, during the study period, there were limited regions with significant LULC changes over northeast China [40,41], which could not dominate the general trend of vegetation phenology. Second, land cover change from natural vegetation to cultivated vegetation mainly occurred in the early stages of the study period, which promises that most samples in the phenological parameters time series were derived from the same vegetation type. Thus, the time series trend could be basically recognized as the tendency for one vegetation. Third, the study of vegetation phenology mainly focused on the trends of phenological metrics and not on their absolute dates, so vegetation type is not the key issue. Fourth, due to the relatively coarse resolution of satellite images (8 km in GIMMS NDVI3g data), remote sensing-based vegetation monitoring can only

Uncertainty Analysis
In the process of calculating long-term vegetation phenology, land use/land cover (LULC) change must be mentioned due to its impact on the accuracy of phenological analysis [39]. Land cover change, including human-induced and natural forcing changes, will introduce heterogeneity in the time series of vegetation greenness (e.g., NDVI) or structure [1]. Additionally, current phenological metrics extraction methods cannot identify information regarding vegetation type. Thus, the multi-year time series of NDVI may represent mixes of different ecosystem types, such as forest and grassland in different periods. During the period 1982-2015, there were some changes of LULC in the study area. In the period 1980-2000, these LULC changes were primarily changes of forest and grassland to pasture in the eastern area, the change of farming-pastoral land to farmland, and the change of forest to grassland in the northern area. After the year 2000, there were mainly land cover changes between different farmland planting systems, such as the transition from non-irrigated farmland to paddy field in the Northeast Plain, with a small amount of grassland converted to farmland in the central area [40].
These LULC changes not only result in the complexity of vegetation index time series but also further induce uncertainties in computing vegetation phenological metrics. However, the results in this paper can be considered reliable. First, during the study period, there were limited regions with significant LULC changes over northeast China [40,41], which could not dominate the general trend of vegetation phenology. Second, land cover change from natural vegetation to cultivated vegetation mainly occurred in the early stages of the study period, which promises that most samples in the phenological parameters time series were derived from the same vegetation type. Thus, the time series trend could be basically recognized as the tendency for one vegetation. Third, the study of vegetation phenology mainly focused on the trends of phenological metrics and not on their absolute dates, so vegetation type is not the key issue. Fourth, due to the relatively coarse resolution of satellite images (8 km in GIMMS NDVI3g data), remote sensing-based vegetation monitoring can only represent a comprehensive status for various vegetation types in a pixel, which is suitable for exploring phenological change at the landscape scale. In this case, the phenological information derived from this dataset should be a comprehensive proxy for different vegetation types. Thus, the coarse spatial resolution may contribute to mitigate the impact of land cover change on landscape phenology. Fifth, the spatial distribution of the temporal trend of SOS and PEAK has the prominent feature of regional aggregation, indicating that the impact of land cover change on the phenological trend is not dominant. In addition, the majority of previous studies do not mention the effect of LULC on the phenological dynamics due to the complexity of their patterns [20,[42][43][44]. Based on the analysis above, we can confirm that the findings in this study are rational.
Additionally, there may be potential issues on the AVHRR NDVI3g data source that may induce uncertainties for spatiotemporal patterns of phenological indices. A number of aspects potentially introduce noise in the NDVI data set due to the change of AVHRR sensors and data processing [45,46]. For instance, GIMMS NDVI shows a tendency towards higher positive regression slope values than Terra MODIS in more humid areas [47]. This may be attributed to the longer time range of GIMMS NDVI than MODIS NDVI, as well as the degeneration of the AVHRR sensor. In this study, the data source could have some impacts on the phenological trend but would not reject the significantly statistical trends for key phenological parameters, such as the trends of PEAK and POP (Figure 4). Discrepancy in the absolute values for various NDVI data sources will not change the pattern of the seasonal growth curve that is the focus point in this study.

Conclusions
This study employed a long-term NDVI time series to explore the spatiotemporal dynamics of vegetation growth peak and the corresponding drivers over northeast China. For the entire region, both the date and magnitude of the vegetation growth peak present a slightly increasing trend and an underlying cyclic period of 11 years. The date of the growth peak shows a negative (earlier) trend in croplands, while the peak magnitude has a positive (increasing) trend. The spring green-up date is a key element impacting the date and magnitude of the vegetation growth peak. Generally, an increasing growth peak could improve the accumulation of vegetation productivity. In the future, the vegetation growth pattern may be classified into two types: (1) POP and SOS advancing synchronously, and (2) POP and SOS delaying together; both types have an amplifying PEAK. These findings and concepts could represent a pathway to improve our understanding of the effect of maximum photosynthesis on plant production and its response to climate change.