Spatio-Temporal Evolution, Future Trend and Phenology Regularity of Net Primary Productivity of Forests in Northeast China

Abstract

Net Primary Productivity (NPP) is one of the significant indicators to measure environmental changes; thus, the relevant study of NPP in Northeast China, Asia, is essential to climate changes and ecological sustainable development. Based on the Global Production Efficiency (GLO-PEM) model, this study firstly estimated the NPP in Northeast China, from 2001 to 2019, and then analyzed its spatio-temporal evolution, future changing trend and phenology regularity. Over the years, the NPP of different forests type in Northeast China showed a gradual increasing trend. Compared with other different time stages, the high-value NPP (700–1300 gC·m⁻²·a⁻¹) in Changbai Mountain, from 2017 to 2019, is more widely distributed. For instance, the NPP has an increasing rate of 6.92% compared to the stage of 2011–2015. Additionally, there was a significant advance at the start of the vegetation growth season (SOS), and a lag at the end of the vegetation growth season (EOS), from 2001 to 2019. Thus, the whole growth period of forests in Northeast China became prolonged with the change of phenology. Moreover, analysis on the sustainability of NPP in the future indicates that the reverse direction feature of NPP change will be slightly stronger than the co-directional feature, meaning that about 30.68% of the study area will switch from improvement to degradation. To conclude, these above studies could provide an important reference for the sustainable development of forests in Northeast China.

1. Introduction

In order to restore the fragile ecosystem, reduce environmental pollution, and reduce restrictions of natural resources on socioeconomic development, the 18th National Congress of the Communist Party of China proposed a strategic plan that actively promotes the construction of an ecological civilization by integrating ecological aspect into the entire process of economic, political, cultural and social construction. However, the realization of an ecological civilization is a complex and comprehensive project [1,2], which not only involves politics, economy, culture and pollution control and other aspects, but also is closely related to the sustainable development of forest ecosystem. In the past few decades, many scholars have demonstrated that the forest ecosystems in the middle and high latitudes of the Northern Hemisphere play an important role in the global carbon sink, using the matrix model [3] and TRIPLEX1.0 model [4], based on national forest resource inventories and Landsat time-series data, especially in North America and European [3,4,5,6]. However, as an important player in the harmony between human and nature, forest ecosystem in East Asia also plays an irreplaceable role in maintaining global biodiversity, regulating climate, conserving water, and providing biological resources and other economic goods [7,8,9,10]. With the acceleration of economic development and industrialization, the environmental pollution and ecological damage have become increasingly serious and even the stability of forest ecosystem has been severely damaged. Therefore, the sustainable development of forest ecosystem has gradually and globally become of concern [11,12,13].

As an important factor of forest ecosystem, vegetation Net Primary Productivity (NPP) is one of the indispensable core parameters for estimating terrestrial ecosystem function to explore the impact, response and countermeasures of global changes on ecosystem [14]. Moreover, the NPP, which is the main indicator of ecosystem ecological balance and ecosystem health, is closely related to major issues of concern such as carbon cycle and carbon disturbance, land use change, climate change and natural resource management. In addition, changes in regional phenology regularity, especially the data from the start of the vegetation growth season (SOS) and end of the vegetation growth season (EOS), not only influence terrestrial carbon and water cycle, but also exert significant impact to climate by regulating the land-atmosphere energy exchange [15,16]. Therefore, studies on how accurately simulate and predict the NPP and explore the phenology regularity in the regional scale ecosystem are of great significance [17,18]. In particular, they can provide an important basis for the dynamic carbon cycle in the atmosphere, the vegetation and soil continuum as well as the sustainable management of forest ecosystem [19,20,21].

In recent years, the rapid development of remote-sensing techniques has promoted its wide application in vegetation dynamic monitoring and prediction. Plenty of studies worldwide have conducted extensive studies on estimating forest biomass [22,23,24], the analysis on spatio-temporal distribution of NPP [25], the responses of NPP to climate changes in the past decades [26,27]. The Light Utilization Energy (LUE) models are commonly used to estimate and evaluate vegetation NPP on a regional or global scale. On the one hand, this kind of model simplifies the physiological and ecological process of vegetation, and the calculation is relatively simple. On the other hand, this model can be used to estimate the temporal and spatial variation of vegetation NPP on a large regional or global scale with high accuracy. These advantages promote the development and research of such models. Currently, there are four LUE models of NPP estimation models, such as the Vegetation Photosynthesis Model (VPM) [28,29], Carnegie–Ames–Stanford Approach (CASA) [30], Global Production Efficiency Model (GLO-PEM) [31] and Moderate Resolution Imaging Spectroradiometer (MODIS) GPP algorithm [32].

Among these LUE models, the CASA model has been widely used, whose principle of that is to estimate NPP by vegetation photosynthetic effective radiation and actual utilization of light energy. For example, Bao et al. [33] proposed an improved CASA model to achieve a comparable result of NPP with the results from the early CASA model and ground measured NPP. Matsushita et al. [26] analyzed quantitatively the NPP of different vegetation types on a monthly scale by inputting different climate datasets into the ecological process model, and then further explored the influence of different climate variables on NPP. Most of the abovementioned studies focused on the quantitative analysis and on the changes of the previous vegetation cover at inter-annual and monthly scales, lacking the analysis on a more accurate timescales, which was not conducive to obtain the information of phenology change. In addition, several recent studies have highlighted that there are many uncertainties and systematic errors when implementing all these models for estimating forest NPP in some specific ecological regions [34,35], only considering the effect of water stress in soil or air on vegetation. Zhang et al. [36] also found that the GLO-PEM performed better by comparing the four LUE models to estimate terrestrial gross primary production. The GLO-PEM developed by Prince and Goward [31] is driven by the variables derived almost entirely from satellite products (please see Section 3.2.1 for specific explanations). To sum up, the past studies most concentrate on the past analysis of NPP and comparison of different models to estimate NPP [35], not considering the phenology regularity of different crops and future sustainability analysis.

As a typical region in the Northeast Asia, the forest in Northeast China is not only a significant ecological barrier for urban agglomerations, but it is also the part of northeast forest belt of the national ecological pattern named “two screens and three belts” [37], as they accounted for over 30% of the total forest area in China. Although the Chinese government has implemented a series of ecological construction projects, such as “Returning Farmland to Forest” or “Grassland and Three North Revetment Forests”, several ecological and environmental issues still exist [38], such as the declining forest productivity and deteriorating water quality and so on. In order to better solve the above problems, this study takes Northeast China as an example to explore the ecological environment changes based on the past limitations. Our objectives were to do the following: (1) improve the inversion estimation of forest NPP on the timescale of eight days, by using the improved GLO-PEM and MODIS remote-sensing data; (2) further explore the phenology regularity and its sustainability trend in the future, combined with the evolution characteristics of its spatio-temporal pattern; and (3) reveal the ecological environment quality in the newest year, which is from 2016 to 2019, and its future challenge in Northeast China. Here ecological environment quality reflects the suitability of ecological environment for human survival and sustainable social and economic development at the ecosystem level. These will provide a reference for decision-making about the rational exploitation and utilization of regional forest resources, as well as the slowing down of the global warming trend and cleaner production energy for the surrounding area.

2. Study Area

In this work, Changbai Mountain and Zhangguangcai Mountain in Northeast China, which are located at the edge of the Eurasian continent, are selected as study areas (Figure 1). The whole area spans the warm temperate zone, the middle temperate zone and the cold temperate zone from south to north, respectively, with a typical temperate monsoon climate [39,40,41]. It is one of the coldest regions in China because of its relatively high latitude, and its annual average temperature is −7 to 3 °C, especially in Changbai Mountain. Moreover, its annual temperature difference is larger than the region with the same latitude. The latter is characterized by with high temperature and abundant precipitation during the same period, whereas the annual precipitation decrease as the distance to the shoreline in the southeast area increases [42]. Moreover, the whole study area is also an important forests production base in China, and it is dominated by cold temperate and temperate mountain coniferous forest, mixed forest of temperate coniferous and deciduous broadleaf and temperate deciduous broad-leaved forests in the south of Changbai Mountain. In the central area of Zhangguangcai Mountain, it mainly consists of temperate deciduous shrub, temperate deciduous broad-leaved forest, cold temperate and temperate mountain coniferous forest. In short, it has become an important area for the study of global changes and terrestrial ecosystems due to its typical variation of heat and humidity as well as its unique distribution pattern of vegetation [43,44].

3. Data and Methodology

3.1. Data Source and Processing

3.1.1. Modis Data

In this study, considering the difference of temporal and spatial resolution, we apply the NDVI data in our study that were derived from the MODND1D product (the temporal resolution is one day and the spatial resolution is 500 m) acquired by the Terra satellites in 2012, which is processed by daily MOD09GA (Collection 6 and the spatial resolution is 500 m) through splitting, cutting, projection conversion, unit conversion and other processes. The NDVI data (MODND1D product) from 2001 to 2019 were derived from the Geospatial Data Cloud website (http://www.gscloud.cn/). They were used to obtain the monthly maximum NDVI dataset, at a spatial resolution of 500 m, in the selected years. By contrast, the NDVI product MOD13A1 (Collection 6, the spatial resolution is 500 m), which is acquired by the Terra satellites, has a temporal of 16 days, and the Global MYD13A1 product is also provided every 16 days, at 500-m spatial resolution, as a gridded level-3 product in the Sinusoidal projection, which is acquired by the Aqua satellites (https://ladsweb.modaps.eosdis.nasa.gov/). In addition, this study used other MODIS products with a spatial resolution of 500 m and a time resolution of 8 days, such as fraction of absorbed photosynthetically active radiation (FPAR) from Terra Product (MOD15A2H) and surface reflectance data from Terra Product (MOD09A1). Moreover, surface reflectance data were used to calculate the Land Surface Water Index (LSWI) to require the limiting coefficient of water on photosynthesis of vegetation. These data were obtained from Nasa Earth Science Data (https://ladsweb.modaps.eosdis.nasa.gov/) and were spliced and projected using the MRT tool.

In order to verify the accuracy of GLO-PEM in simulating NPP, the NPP at the spatial resolution of 500 m in 2003 and in 2008 (MOD17A3HGF Version 6 product) was derived from the NASA website (https://cmr.earthdata.nasa.gov/). Here annual NPP is derived from the sum of all 8-day Net Photosynthesis (PSN) products ([MOD17A2H] (https://0-doi-org.brum.beds.ac.uk/10.5067/MODIS/MOD17A2H.006)) from the given year. The PSN value is the difference of the Gross Primary Productivity (GPP) and the Maintenance Respiration (MR).

3.1.2. Meteorological Data

Meteorological data, including daily mean temperature, sunshine time, solar radiation, etc., from 2001 to 2019, were acquired from the China Meteorological Administration (http://data.cma.cn/). In order to maintain a consistent resolution between meteorology data and remote-sensing data from space, the Anuspline interpolation method was used to interpolated point data into raster data with the same spatial resolution as NDVI, FPAR and other data based on 25 meteorological stations in Northeast China.

3.1.3. Land-Use and Land-Cover Data

Land-use and land-cover data, encompassing 2000, 2005, 2010 and 2015 in Heilongjiang and Jilin provinces, were provided by the Resources and Environmental Science Data Center of the Chinese Academy of sciences (http://www.resdc.cn/). These data were obtained based by manual visual interpretation, and they were classified in more than 43 categories. In this study, the land-use and land-cover data were regrouped into five categories, according to a new classification logic, resulting in (i) temperate coniferous forests, (ii) cold temperate and temperate mountain coniferous forests, (iii) temperate deciduous broad-leaved forest, (iv) mixed forest of temperate coniferous and deciduous broadleaf, and (v) temperate deciduous shrub.

3.2. Methodology

3.2.1. GLO-PEM

In this study, the Global Production Efficiency Model (GLO-PEM), which was based on light use efficiency, was applied to simulate the NPP of forests in Northeast China. This model fully considered the characteristics of plant growth and abiotic environmental factors that could regulate plant photosynthesis, such as temperature, soil and water conditions. The NPP estimated by the GLO-PEM (Figure 2) was mainly determined by considering the absorbed photosynthetic active radiation (APAR) and the light use efficiency (ε), using Equation (1) [36,45,46,47,48]:

$$\mathrm{NPP}\left(\mathrm{x}, \mathrm{t}\right) = \mathrm{APAR}(\mathrm{x}, \mathrm{t}) \times \mathsf{\epsilon }(\mathrm{x}, \mathrm{t})$$

(1)

where NPP(x, t) represents net primary productivity of pixel x at time t, APAR(x, t) refers to the absorbed photosynthetic active radiation of that and ε(x, t) indicates light use efficiency under the same conditions. The APAR mainly depends on the following two variables:

$$\mathrm{APAR}(\mathrm{x}, \mathrm{t}) = \mathrm{FPAR}(\mathrm{x}, \mathrm{t}) \times \mathrm{PAR}(\mathrm{x}, \mathrm{t}).$$

(2)

where PAR(x, t) refers to photosynthetically active radiation of pixel x at time t, FPAR(x, t) is fraction of absorbed photosynthetically active radiation of that.

Photosynthetically Active Radiation (PAR)

In this study, PAR was estimated by considering the geographic location of typical sites in the study area, including longitude, latitude, elevation and the sunshine hours as the input data for the crop mode [36,49,50].

(1)

Solar radiation top of atmosphere (H₀)

$${\mathrm{H}}_0{ = (1/\mathsf{\pi }) \times \mathrm{G}}_{\mathrm{SC}}{ \times \mathrm{E}}_0 \times (\cos \mathsf{\varphi }{\cos \mathsf{\delta }\sin \mathrm{W}}_{\mathrm{S}}{ + (\mathsf{\pi }/180) \times \sin \mathsf{\Phi }\sin \mathsf{\delta } \times \mathrm{W}}_{\mathrm{S}}).$$

(3)

where G_sc is the solar constant with value equal to 1367 W·m⁻² that is equivalent to 118.108 MJ· m⁻²·d⁻¹. E₀ represents eccentricity correction factor of the earth’s orbit; Φ, δ, W_s indicate latitude, declination of the sun and time Angle, respectively.

(2)

Solar radiation in clear sky (H_L)

According to existing studies [51], the transparency coefficient of total radiation in the atmosphere is usually about 0.8 under ideal conditions, i.e., clear sky, and it is different under specific environmental conditions.

$${\mathrm{H}}_{\mathrm{L}} = 0{.8 \times \mathrm{H}}_0$$

(4)

(3)

Daily solar radiation (H)

$${\mathrm{H} = \mathrm{H}}_{\mathrm{L}}{ \times (\mathrm{a} + \mathrm{b} \times \mathrm{S}/\mathrm{S}}_{\mathrm{L}})$$

(5)

where H is the daily solar radiation; H_L is solar radiation in clear sky; S and S_L are sunshine hours and sunshine days, respectively; and a and b are empirical parameters generally obtained by regression simulation based on measured values of solar radiation. According to the average value of total daily solar radiation, the sunshine percentage and the solar radiation in clear sky in different regions of China. Podesta et al. [52] and Wu et al. [53] calculated that a was equal to 0.248, and b was equal to 0.752. PAR is calculated based on the daily solar radiation; the formula is expressed as Equation (6).

$$\mathrm{PAR} = 0.5 \times \mathrm{H}$$

(6)

Light-Use Efficiency (ε)

The light-use efficiency under realistic conditions (ε) refers to the light use efficiency under the stress of environmental driving factors such as temperature and vapor pressure, which can be extracted from meteorological indicators [28,29,30,45,46]. ε is calculated by using Equation (7):

$${\mathsf{\epsilon }(\mathrm{x}, \mathrm{t}) = \mathrm{T}}_{\mathsf{\epsilon }1}{(\mathrm{x}, \mathrm{t}) \times \mathrm{T}}_{\mathsf{\epsilon }2}{(\mathrm{x}, \mathrm{t}) \times \mathrm{W}}_{\mathsf{\epsilon }}(\mathrm{x}, \mathrm{t}) \times \mathsf{\epsilon }*$$

(7)

where ${\mathrm{T}}_{\mathsf{\epsilon }1}\left(\mathrm{x}, \mathrm{t}\right)$ and ${\mathrm{T}}_{\mathsf{\epsilon }2}(\mathrm{x}, \mathrm{t})$ represent the stress coefficients of high temperature and low temperature on light utilization efficiency, respectively.${ \mathrm{W}}_{\mathsf{\epsilon }}\left(\mathrm{x}, \mathrm{t}\right)$ represents the influence coefficient of water on the utilization of light energy. ε* is the maximum light use efficiency under ideal conditions (please see References [28,29,30,45,46] for specific calculations).

3.2.2. Sen–Mann–Kendall Trend

In order to quantitatively study the variation trend of forests NPP, Sen–Mann–Kendall trend analysis was adopted in this study to comprehensively reflect the evolution characteristics of regional pattern of NPP during the period from 2001 to 2019 by exploring the spatial variation trend and variation slope of each grid [54,55]. The calculation formula is as follows:

$$\mathrm{Slope}=\frac{\mathrm{n}\times {{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{n}}\mathrm{i}\times {\mathrm{NPP}}_{\mathrm{i}}-{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{n}}\mathrm{i}{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{n}}{\mathrm{NPP}}_{\mathrm{i}}}{\mathrm{n}\times {{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{n}}{\mathrm{i}}^2-{\left({{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{n}}\mathrm{i}\right)}^2}$$

(8)

where Slope is the linear fitting slope of the equation, the positive and negative value of the slope represents the increase and decrease of forests NPP in Changbai Mountain area; NPP_i is the total annual NPP in the i year; and N is the span of research period of 2001–2019.

3.2.3. Hurst Index

Hurst index can not only be used to describe the sustainability of the NPP in Northeast China but also predict the future development trend from time series. Therefore, this study attempts to explore the future ecological environment changes of forests in Northeast China by using Hurst index to provide reference for the sustainable management of forests in the future. Hurst index (H), which is based on the re-scale range (R/S) analysis, is an effective method to describe quantitatively the long-term dependence of time series information [56,57]. The following sequences are all established at any moment $\left\{\mathrm{i}=1,2,\dots \dots , \mathrm{n}\right\}$after given time series variable of NPP. For any positive integer m, define the above time series. Here this study mainly predicts the future changes based on the NPP from 2001 to 2019, so n is 19.

Difference sequence:

$$\mathsf{\Delta }{\mathrm{NPP}}_{\mathrm{i}}{ = \mathrm{NPP}}_{\mathrm{i}}{ - \mathrm{NPP}}_{\mathrm{i}-1}.$$

(9)

Mean sequence:

$$\overline{\mathsf{\Delta }\mathrm{NPP}\left(\mathrm{m}\right) } = \frac{1}{\mathrm{m}}{\displaystyle \sum }_{\mathrm{i} = 1}^{\mathrm{m}}\mathsf{\Delta }{\mathrm{NPP}}_{\mathrm{i}}\hspace{1em}\hspace{1em}\hspace{1em} \mathrm{m} = 1,2 \dots \mathrm{n}$$

(10)

Cumulative dispersion:

$$\mathrm{X}\left(\mathrm{t}\right) = {\displaystyle \sum }_{\mathrm{i} = 1}^{\mathrm{m}}(\mathsf{\Delta }{\mathrm{NPP}}_{\mathrm{i}}-\overline{\mathsf{\Delta }\mathrm{NPP}\left(\mathrm{m}\right)})\hspace{1em}\hspace{1em}\hspace{1em} 1 \le \mathrm{t} \le \mathrm{m}$$

(11)

Range:

$$\mathrm{R}\left(\mathrm{m}\right)=\mathrm{X}{\left(\mathrm{t}\right)}_{\max }-\mathrm{X}{\left(\mathrm{t}\right)}_{\min }\hspace{1em}\hspace{1em}\mathrm{m}=1,2\dots \mathrm{n}$$

(12)

Standard deviation:

$$\mathrm{S}\left(\mathrm{m}\right)={\left[\frac{1}{\mathrm{m}}{\displaystyle \sum }_{\mathrm{i}=1}^{\mathrm{m}}{(\mathsf{\Delta }{\mathrm{NPP}}_{\mathrm{i}}-\overline{\mathsf{\Delta }\mathrm{NPP}\left(\mathrm{m}\right)})}^2\right]}^{\frac{1}{2}} \mathrm{m}=1 \dots \mathrm{n}$$

(13)

This formula $\mathrm{R}/\mathrm{S} = \mathrm{R}\left(\mathrm{m}\right)/\mathrm{S}\left(\mathrm{m}\right)$ is defined based on computing $\mathrm{R}\left(\mathrm{m}\right)$. and $\mathrm{R}\left(\mathrm{m}\right)$. When $\mathrm{R}/\mathrm{S}\propto {\mathrm{T}}^{\mathrm{H}}$, Hurst phenomenon exists in the considered time series, and the representational significance of Hurst index is evident. When 0.5 < H < 1, the times series is persistent, meaning that the future trend is consistent with the past trend. The closer the value is to 1, the stronger the continuity will be. However, when 0 < H < 0.5, the time series has an anti-continuity, indicating that the future change trend is opposite to the past change trend. The closer the value is to 0, the stronger the anti-continuity is. When H = 0.5, the event sequence is a random sequence, and the future change trend is unpredictable [41].

4. Model Validation

As the GLO-PEM is a light use efficiency model, the estimation of NPP mainly depends on light use efficiency and photosynthetically active radiation. Therefore, this study applied the method of phased calculation and comparison for verification.

4.1. Comparisons of Solar Radiation and PAR Variation of Different Research Sites

In this study, to evaluate the accuracy of the PAR, the estimated PAR was compared with the ground observation data obtained from meteorological radiation stations in 2008. Firstly, the solar radiation top of atmosphere, solar radiation in clear sky and daily solar radiation were estimated based on the model. As shown by the analysis on the solar radiation of six typical stations (Figure 3), it was evident that the distribution results of the stations basically conformed to the typical theoretical distribution model, and the simulated values of solar radiation in clear sky were consistent with the actual daily maximum values observed over many years. However, due to differences in atmospheric composition, cloud cover, atmospheric moisture content and atmospheric suspended matter content, the absorption and scattering of daily solar radiation top of the atmosphere were different, and this affected the total radiation reaching the ground. The correlation coefficient between the estimated PAR and the observed data surpassed 0.90 in Harbin and Yanji stations (Figure 4). The results indicated that the simulation accuracy of PAR generally met the requirements.

4.2. NPP Validation between Simulated NPP and MOD17 NPP

To further assess the model accuracy, this study also performed a comparative analysis on the estimated NPP based on GLO-PEM and MOD17 NPP products in 2003 and 2008 (Figure 5). The results indicated that the NPP values estimated by the model were generally consistent with the ones of MOD17 NPP, and the correlation coefficients at the pixel scale in 2003 and 2008 were 0.91 and 0.89, respectively. Therefore, the accuracy of NPP estimated by this model, which could represent the distribution of forests NPP of Northeast China, was relatively high.

5. Results

5.1. Spatio-Temporal Evolution of NPP in Northeast China

5.1.1. Variation Characteristics of Forests NPP from 2001 to 2019

In order to reveal the evolution characteristics of NPP in northeast forests over time, this study selected the average NPP as an indicator of ecological environment quality from 2001 to 2019. Overall, this study found an increasing trend with an average increase of 3.95 gC·m⁻²·a⁻¹ (Figure 6) in the analyzed years. Specifically, there was a downward trend from 2001 and 2005, with a decrease of 19.85%. Then, due to the large-scale implementation of natural forests protection and other ecological restoration projects in Northeast China, the ecological environment quality rapidly increased during the period of 2005–2008 and the NPP reached 859.16 gC·m⁻²·a⁻¹ in 2008. However, the average NPP did not show a trend of continuous increase after 2008. Instead, it showed a severe decline with a reduction of −17.22% from 2008 to 2010, which may be related to large-scale drought events. Subsequently, the NPP gradually increased again from 2010 to 2015, mainly due to the climate and the national attention to the protection of natural forests. During the period of 2016 to 2019, the mean NPP had a slight drop in 2016, 2017. Subsequently, the mean NPP reached the maximum 976.17 gC·m⁻²·a⁻¹ in 2019. It is likely to be affected by a range of climate changes, such as rising regional temperatures, permafrost degradation, shorter freezing–thawing times and so on [58,59,60].

5.1.2. Inter-Annual Variation of NPP of Different Forests Types

Next, this study estimated the NPP of different kinds of forests to reveal the variation characteristics of these forests. The average value of five forests NPP in Northeast China from 2001 to 2019 were ranked from high to low as follows: temperate coniferous forest > mixed forest of temperate coniferous forest and deciduous broad-leaved mixed forest > temperate deciduous broad-leaved forest > cold temperate and temperate mountain coniferous forest > temperate deciduous shrub. Compared with other four forests, the NPP of cold temperate and temperate mountain coniferous forest showed a higher increasing trend, with a rate of 8.89 gC·m⁻²·a⁻¹. The reason is likely that global warming has become an unstoppable global climate problem in the last 20 years, leading to a significant increase in NPP [59]. Similarly, the NPP of other forest types all showed an increasing trend of different amplitude. The remaining forest types with different increasing rate from high to low is Mixed forest of temperate coniferous and deciduous broadleaf, temperate coniferous forests, temperate deciduous broadleaved forest and temperate deciduous shrub, whose rate is 8.08, 7.41, 6.34 and 6.05 gC·m⁻²·a⁻¹, respectively.

From the NPP changes of different forests in different time periods (Figure 7), it was evident that the NPP of various forests showed a decreasing trend from 2001 to 2005, and then went up gradually and reached the highest point in 2007 and 2008. However, the NPP declined significantly again from 2009, mainly because of the severe summer drought. In recent years, that is from 2017 to 2019, the NPP of these forests increased significantly. Specifically, the temperate deciduous broad-leaved forest, which accounted for 49.14% of the study area, its average NPP reached 1067.91 gC·m⁻²·a⁻¹. During the period of 2001 to 2019, there are two peaks for this forest type, which peaked at 942.75 gC·m⁻²·a⁻¹ in 2008 and 2067.91 gC·m⁻²·a⁻¹ in 2018, respectively. As for the cold and temperate mountain coniferous forest, mixed forest of temperate coniferous and deciduous broadleaf, their change trends of NPP were overall consistent with the one of temperate deciduous broad-leaved forest. Among them, the NPPs of the two types of forests decreased significantly during the period of 2004–2005 and 2008–2009. In the first stage, they varied from 803.66 gC·m⁻²·a⁻¹ to 753.08 gC·m⁻²·a⁻¹ and from 902.31 gC·m⁻²·a⁻¹ to 769.40 gC·m⁻²·a⁻¹_, respectively. Next, the NPP in the second phrase declined more seriously than previously, with a rate of 18.32% and 16.06%, respectively. The highest NPP values reached 951.19 gC·m⁻²·a⁻¹ and 1013.29 gC·m⁻²·a⁻¹ in 2008, respectively. However, the NPP inter-annual changes in temperate deciduous shrub were relatively stable compared with the other types of forests.

5.1.3. Spatial Evolution of NPP, from 2001 to 2019, in Northeast China, in Different Periods

Based on the GLO-PEM, this study firstly estimated the NPP of Northeast China from 2001 to 2019, and then analyzed its spatio-temporal evaluation in different stages, namely the phases 2001–2005, 2006–2010, 2011–2015 and 2016–2019 (Figure 8). The dynamic evolution of NPP in the study area showed a clear spatial heterogeneity, with a decreasing trend from southeast to northwest. Overall, the average NPP of the study area was 778.54 gC·m⁻²·a⁻¹ and the change of NPP in the high-value areas was from 700 to 1300 gC·m⁻²·a⁻¹, with a small fluctuation every five years. The regions with high NPP were mainly located in the east of Changbai Mountain and the south of Zhangguangcai Mountain, which were mainly covered with deciduous broad-leaved forest, mixed forest of temperate coniferous and deciduous broad-leaved, as well as temperate deciduous shrub. Because this region has favorable hydrothermal conditions affected by the monsoon climate and it is located in the humid zone of the mid-temperate zone, thus resulting in better vegetation growth. For the early stage, the median-area NPP (700–900 gC·m⁻²·a⁻¹) increased by 10.36% from the first stage (2001–2005) to the second stage (2006–2010), mainly in the northwest of Bin County in Zhangguangcai Mountain, the west of Wuchang City, Acheng City, Ningan City, Jilin City and the west of Shulan City. The reason is likely that the lower-value area (0–500 gC·m⁻²·a⁻¹) gradually moved to the middle value area over time, only accounting for 9.05% of the study area during the period of 2006–2010. Compared with other different time stages, the high-value NPP (700–1300 gC·m⁻²·a⁻¹) in Changbai Mountain from 2017 to 2019 occupied a larger area. For instance, the NPP had an increasing rate of 6.92% compared to the stage of 2011–2015. This indicated that the ecological environment in Changbai Mountains had recovered quickly due to the implementation of natural forest protection project in recent years.

5.2. Future Trend of Forests NPP in Northeast China

5.2.1. Analysis on Spatial Changing Trend of NPP, from 2001 to 2019, in Northeast China

In order to explore the changing trend variation of NPP in Northeast China, this study estimated the Sen trend degree of NPP from 2001 to 2019 (z = 0.0990, slope = 0.2982). The overall NPP in the study area presented a not-significant increasing trend from 2001 to 2019, and the area with an increase and decrease trend accounted for 52.58% and 47.42%, respectively (Figure 9a). In terms of varying rates, the area with the NPP change rate between −10 and 0 gC·m⁻²·a⁻¹ accounted for 44.36% of the whole area and it was mainly distributed in the south of Zhangguangcai Mountain and the north of Changbai Mountain. To a certain degree, the ecological environment in many regions showed a weak trend of decline. Additionally, there was a large proportion of the area with a change rate between 0–5 gC·m⁻²·a⁻¹, mainly located in the western part of Wuchang City, the north of Yushu City, and the eastern part of Ningan City in Zhangguangcai Mountain. On the contrary, the regions with the change range of NPP between 5–10 gC·m⁻²·a⁻¹ and over 10 gC·m⁻²·a⁻¹ accounted for 14.37% and 2.19% of the whole area, respectively, and they mainly distributed around the edges of the study area within the spatial pattern of belts. The regions with a change rate under −10 gC·m⁻²·a⁻¹ only made up 3.07% and they were mainly distributed in the central part of Wuchang City, the western part of Shangzhi City, the eastern part of Jiahe City and the western part of Dunhua City in the central part of Zhangguangcai Mountain. Thus, only few regions showed very serious ecological environment degradation.

In order to further determine whether the above regions showed a significant trend, we firstly performed the Mann–Kendall (M–K) test, and divided the significance test results into significant changes (Z > 1.96 or Z < −1.96) and insignificant changes (−1.96 ≤ Z ≤ 1.96) at 0.05 confidence level. Then, the NPP change trend at the pixel scale was calculated by superposing the changing trend results of Theil–Sen Median and significance results by M–K test, and by further dividing the results into four types (Figure 9b). The outcome demonstrated that the NPP change trend in most Northeast China was not significant from 2001 to 2019, and the areas with significant change only accounted for 3.99% of the whole study area. In terms of the upward trend, 50.78% of the areas showed a slight improvement, and only 1.79% showed a significant improvement. The areas with extremely significant improvement in ecological environment were mainly present in Bin County, the middle of Yanshou County, Acheng City, the west of Wuchang City, the east of Shuangcheng City, and distributed in the west of Shulan City in Zhangguangcai Mountain. The areas with slight improvements were mainly distributed in the periphery of the whole area. In terms of declining trend, 45.22% of the regions in Northeast China showed a slight degradation trend. They were mainly distributed in the south of Shulan City, east of Wuchang City, west of Dunhua City, Shangzhi City, and the intersections of Acheng City and Yanshou County in Zhangguangcai Mountain. The severely degraded areas accounted for only 2.20%, mainly distributed in the south of Longjing City, Helong City in the east of the Changbai Mountains, and the Jiahe City in the middle region.

5.2.2. Future Trend of NPP in Northeast China

In order to represent the change trend of ecological environment in the future, this study estimated the Hurst index of NPP in Northeast China in the past 19 years. The Hurst index in Northeast China will fluctuate between 0.17 and 0.93, and the average Hurst index will be 0.49 (Figure 10). The regions with Hurst index less than 0.50 will account for 55.79%, indicating that the area with a reverse-direction NPP change will be bigger than the one with a co-directional feature in Northeast China. Particularly, the pixel number of Hurst index between 0.40 and 0.60 will account for 73.98%, while the pixel number with an index lower than 0.40 and greater than 0.60 will account for 15.78% and 10.24%, respectively. This illustrates that the regions with strong reverse direction feature and strong co-directional feature in Northeast China will be relatively low abundant. The region with high Hurst index (0.60–0.93) will mainly distribute in the northern part of the mixed forest and of the coniferous and deciduous broad-leaved belt in the south of the temperate zone, whereas the region with a low value index (0.16–0.40) will mainly distribute in the northern part of Zhangguangcai Mountain and in the southern part of Changbai mountains.

To better reveal the future trend and sustainability of NPP in Northeast China, this study overlaid the current NPP change trend and Hurst index. About 30.68% of the study area is predicted to switch from improvement to degradation, especially in Yanshou County, Shangzhi City, Fangzheng County, Ningan City and Bin County in the north of Zhangguangcai Mountain, as well as in Fusong County and Changbai Korean Autonomous County in Changbai Mountain range. The areas with continuous improvement will account for 21.99% of the total area and they will mainly distribute in Shuangcheng and Acheng City in the north of Zhangguangcai Mountain and the south of Changbai Mountain. Moreover, about 25.11% of regions will exhibit a descending trend of the NPP that will be reversed in the future. These areas will be clustered in the west of Shangzhi City, the east of Acheng City and the north of Tonghe County in the northwest of Zhangguangcai Mountain, and in strips in Antu County of Changbai Mountain. However, the regions with a continuous degradation will account for 22.22% and they comprise Jiahe City and Dunhua City in the central part of Northeast China.

5.3. Phenology Regularity of Forests over Time in Northeast China

In this study, we mainly explored the changing trend of forest NPP every eight days, from 2001 to 2019, to reveal the phenology regularity in Northeast China (Figure 11). Specifically, the first NPP represents the average accumulated NPP from 2001 to 2019, from the first to the eight day. The results indicated that the period of NPP accumulation in northeast forests mainly ranged from the 145th to 257th days, i.e., from May to September, when good hydrothermal conditions were occurring. The total amount of NPP in this growth period accounted for 81.19% of the annual average. In addition, the inter-annual fluctuation range of NPP in forests during growing season was larger than the one during non-growing season. In terms of seasonal changes, the NPP in northeast forests accumulation mainly took place from the 161st to 241st days, during the summer time, where the cumulative amount of NPP was 489.13 gC·m⁻²·(8d)⁻¹, accounting for 63.88% of the total annual NPP. In the spring (65th–153rd days), the cumulative amount of NPP was relatively high, accounting for 19.26% of the total annual NPP. However, as for the autumn (249th–329th days) and the winter (1st–57th and 337th–353rd days), their amounts of the accumulation NPP were 108.01 and 21.04 gC·m⁻²·(8d)⁻¹, only accounted for 14.11% and 2.75%, respectively. During this period, the plant grew extremely slowly due to the influence of Mongolian high pressure and low temperature in the majority of the Northeast China areas.

According to the average NPP, every eight days, from 2001 to 2019, the phenomenon of phenology advance or lag occurred in all the vegetation growth seasons in different periods. During the period of 2001–2005, there was a remarkable increase of the growth of forests on the 105th–153rd days of 2004 and 2005. On the contrary, the growth season of forests in Northeast China advanced in 2007 and 2008. This was probably due to the extremely warm winter leading to the early increase of NPP during that period. Furthermore, the NPP growth in 2008 had an evident advance, up to 2.5 times the average on the 105th day. During the period of 2010–2015, a significant advance has occurred in the growth of forests. However, by the 153rd day, the NPP accumulation had a significant decrease from 55.80 to 22.97 gC·m⁻²·(8d)⁻¹. From 2016 to 2019, due to the aggravation of global warming and the implementation of ecological engineering, the accumulation of forest NPP increased significantly from day 1 to day 353.

5.3.1. Phenology Regularity of Forest at the Start of the Vegetation Growth Season (SOS)

In order to reveal the phenological regularity of forest, this study intends to explore the changing trend of NPP every eight days in the initial growth period. According to the relevant research [61], the starting growth date of most forests in Changbai Mountain from 2001 to 2019 were mainly 97th days (representing the accumulated NPP of 97th–105th days) to 137th days (representing the accumulated NPP of 137th–145th days), that is, from mid-late April to early May (Figure 12). Firstly, what we need to explain is that 97th NPP represents the average accumulated NPP from day 97th to day 105th from 2001 to 2019. As can be seen Figure 12, the forests entered the rapid growth stage from the 97th day. Because of the rapid rise of the temperature, the NPP also increased significantly. From 2001 to 2019, the forest NPP all showed an apparent increasing trend, especially on the 129th and 137th days. On the whole, the starting date of growing season of forests in Changbai Mountain from 2001 to 2019 has been advanced.

5.3.2. Phenology Regularity of Forest at the End of the Vegetation Growth Season (EOS)

Since the Changbai Mountain, located in the high latitude area in the northeast, is influenced by climate changes, the ending date of forests is relatively late compared to most areas in China. The ending date of forests in Northeast China were concentrated on the 265th–289th days (Figure 13). It can be seen that the accumulated NPP in the ending date of forests growth all presented a slight increase, with an increasing of 0.26, 0.27, 0.11 and 0.11 gC·m⁻²·a⁻¹, respectively. The above results showed that the phenological regularity of forest from 2001 to 2019 had an apparent lag in the ending date in Northeast China.

6. Discussion

6.1. Optimization about Station-Temporal Regularity Exploration

In recent 19 years, although natural forests protection and other ecological restoration projects had been largely implemented in Northeast China, the ecological environment did not all improved persistently during the whole study period. Hence, the construction and ecological restoration of forests in Northeast China still require continuous attention of the government. This study displayed that the NPP in Northeast China from 2001 to 2019 generally exhibited a slow rising trend with certain fluctuations, especially in the period of 2005–2008 and 2017–2019. Similarly, a latest study in Northeast China based on climate simulation scenarios of RCP 4.5 and RCP 8.5, shows an increase from 2011 to 2100, with corresponding reference years (1981–2010) under the future climate change due to the drastic climate change [60].

As for the inter-annual NPP changes of different forests types, on the whole, the NPP of all forest types all showed an increasing trend, especially cold temperate and temperate mountain coniferous forest and mixed forest of temperate coniferous and deciduous broadleaf. Lindholm et al. [61] also performed the NPP in China has a significant increase influenced by the temperature, precipitation and CO₂ concentration. Moreover, there is a larger increase in NPP in broad-leaf and needle-leaf mixed forest in Northeast China [62]. During the period of 2001 to 2019, they experienced three major fluctuations indicating its poor stability. Two reasons may be accounted for this change. On the one hand, large-scale drought events had likely an influence on the NPP. On the other hand, this sensitive area was susceptible to human activities and the transformation of different land use types [4,37,63]. Therefore, the policy of under-forest economy, that is the government firstly could take some intercropping measures between forests in vertical structures, should be introduced in areas with a good forests ecological environment to promote the common development of under-forest economy and forest environment protection. In addition, the government should then establish ecological restoration and comprehensive treatment plan to enhance the ecological service function.

6.2. Summaries of and Sustainability Exploration of Forests and the Outlook for Future Ecological Environment

As a traditional forest region of China and an ecological barrier to revitalize the old industrial base, the decrease of forests NPP in Northeast China will seriously affect the ecological environment and national ecological security pattern [64]. According to our data, the NPP in Northeast China generally showed a minimal increase trend from 2001 to 2019. However, we also found that 47.42% of the area still exhibited a decreasing trend in this study. Notably, the ecological environment still showed a trend of continuous degradation in a large proportion of regions. Therefore, the management of the continuously degraded area should be accentuated in the future ecological environment construction of northeast forests, in order to achieve the sustainable development of forests ecosystem.

The study on the sustainability of NPP in Northeast China in the future indicates that the reverse direction feature of NPP change will be slightly stronger than the co-directional feature. It means that the degraded area will be slightly bigger than the improved areas. On the contrary, another study about Northeast China reported that the combined effects of climate change and CO₂ fertilization on the increase of NPP were estimated to be 10%–20% for the 2030s and 28%–37% in the 2090s, using the improved TRIPLEX model [4]. Differences in the future prediction are mainly attributed to the following reasons. In this study, we focus on estimating NPP every eight days, using GLO-PEM and predicting the future ecological conditions using Hurst index [57,65,66] based on the estimated NPP from 2001 to 2019 instead of estimating the forest NPP in future different scenarios. However, Peng et al. [4] concentrated on estimating the NPP changes in 2030s to 2090s based on the reference baseline year (1999), using the TRIPLEX model, not considering the potential impacts of land-use changes.

In addition, about 30.68% of the study area will switch from improvement to degradation, especially in the Changbai Mountain. Since most of the areas in Yanshou County, Shangzhi City and Fangzheng County in this area are temperate broad-leaved forests, it is advisable to intercrop the cultivated vegetation and the shrubs with higher ecological benefits. In addition, a single forest type may lead to more frequent outbreaks of pests and diseases, which could destroy the forest health, and then seriously limit the sustainable improvement of forest in the future [67]. Hence, the government also should restore the local precious broad-leaved tree resources so as to increase the diversity of forests types. Remarkably, the regions with a trend of continuous degradation will also account for 22.22% of the total area. As the layout of the degraded area is relatively dispersed, natural environmental resources are extremely vulnerable to regional damage. Therefore, more attention should be paid to the pattern optimization of forests types in Northeast China, to improve the construction of ecological civilization.

6.3. Uncertainty Discussions for Phenology Regularity in Northeast China

Phenology, as a comprehensive indicator of climate and natural environment changes, has been widely used worldwide in many assessment projects about environmental change impact. In this study, we determined the SOS, that is to say, the date of green-up date, in the spring phenology, and EOS in the winter phenology of forests based on the statistics references of the NPP time series every eight days, from 2001 to 2019, in Northeast China. Our study found that the whole growth period of forest in Northeast China became extended with the change of phenology. Concretely, our results performed that there is an advanced SOS and a significant delayed EOS in forest ecosystem in Northeast China, from 2001 to 2019. In addition, other scholars have developed many phenology retrieving methods based on remoting-sensing data, such as machine learning and vegetation indices from time series based on software TIMESAT, but the strength and weakness between different retrieving methods are still unclear [68,69,70,71]. A latest study in Northeast China also found that there is an advanced at the start of the growing season (SOS) with in forest ecosystems based on statistics and analysis [72]. However, there are some uncertainties in exploring the phenological period. For example, this study only explored the EOS and SOS of forests roughly based on the eight-day NPP instead of some photo-temperature driven phenological model. Moreover, since we have not yet found the more precise way to study the growth period of each forest type in detail, this also creates uncertainty.

The study of phenology regularity not only improves the understanding of vegetation response to climate change, but also influences the forest construction [16,73,74]. Of course, there are some limitations in the exploration of vegetation phenology. Firstly, we did not explore the spatial phenology since we pay more attention to the analysis from multiple time series. In addition, we recommend that the multi-model retrieving methods based on satellite data should be compared to explore the vegetation phenology in the different ecosystems in the future. To conclude, to realize the sustainable development of forest ecosystem, it is suggested that the forest management plan should be scientifically formulated based on understanding and mastering the growth regularity and phenological periods of different forests in Northeast China.

6.4. Limitations in the NPP Estimation of GLO-PEM and Verification

In summary, this study revealed the phenological regularity by exploring the NPP trend every eight days in the early, middle and ending growth in Northeast China, from 2001 to 2019, based on the GLO-PEM. It could provide new useful knowledge for promoting symbiotic socioeconomic–ecological interactions and regional sustainable development. However, the used model showed few limitations that needed to be further addressed. Firstly, a certain error in the accuracy verification may occur, due to the limited amount of measured data in the northeast forest. However, the simulated results have been verified by using the other remoting-sensing data (MOD17 NPP) in 2003 and 2008, proving that the accuracy of our model is very high (please see Section 4.2 for specific validation). Furthermore, the scale difference between the remote-sensing data and the ground-measured data may also lead to some errors. There were also some uncertainties when running GLO-PEM, mainly considering the temperature, precipitation, land surface water index, different forest types and so on. In the future, the NPP precision could be well enhanced if more useful data were taken into the improved model, such as human distribution factor, CO₂ concentration and other factors.

7. Conclusions

In this study, firstly, we applied the verified GLO-PEM to estimate the NPP every eight days, from 2001 to 2019, based on remote-sensing data and climate data. Studies showed that an increasing trend occurred with an average of 3.95 gC·m⁻²·a⁻¹, especially during the period of 2016–2019. In addition, we explored the vegetation phenology, using the estimated NPP data from the time series at the early and ending date of the growth season. We found that there is a significant advance in SOS and a delayed EOS from 2001 to 2019, causing the prolong of the whole growing season of forests in Northeast China. The above studies have provided a basis for further phenology regularity. In future research, we also need to consider the driving mechanism of vegetation phenology and pay more attention to the inversing methods.

Finally, methods including Sen–Mann–Kendall trend and Hurst index were combined to explore the sustainability analysis and future ecological environment in Northeast China. Generally, more researchers explore the future ecological environment using different climate scenarios based on the older reference year. We proposed a different method to explore the NPP changes from the time series of 2001–2019 in the future. Studies demonstrated that the overall NPP presented a not-significant increasing trend from 2001 to 2019, but 47.42% of the area still exhibited a decreasing trend, mainly in the central part of Zhangguangcai Mountain and the northern part of Changbai Mountain. In the future, the regions with a reverse-direction NPP change will be bigger than the one with a co-directional feature in Northeast China. Moreover, the regions with a continuous degradation will account for 22.22%, which distributed in Jiahe City and Dunhua City in the central part of Northeast China. These conclusions could provide a strong potential template for tropical managers and sustainable management in Asia, even in the global climate change.