Three Decades of Gross Primary Production (GPP) in China: Variations, Trends, Attributions, and Prediction Inferred from Multiple Datasets and Time Series Modeling

Abstract

The accurate estimation of gross primary production (GPP) is crucial to understanding plant carbon sequestration and grasping the quality of the ecological environment. Nevertheless, due to the inconsistencies of current GPP products, the variations, trends and short-term predictions of GPP have not been sufficiently well studied. In this study, we explore the spatiotemporal variability and trends of GPP and its associated climatic and anthropogenic factors in China from 1982 to 2015, mainly based on the optimum light use efficiency (LUEopt) product. We also employ an autoregressive integrated moving average (ARIMA) model to forecast the monthly GPP for a one-year lead time. The results show that GPP experienced an upward trend of 2.268 g C/m² per year during the studied period, that is, an increasing rate of 3.9% per decade since 1982. However, these trend changes revealed distinct heterogeneity across space and time. The positive trends were mainly distributed in the Yellow River and Huaihe River out of the nine major river basins in China. We found that the dynamics of GPP were concurrently affected by climate factors and human activities. While air temperature and leaf area index (LAI) played dominant roles at a national level, the effects of precipitation, downward shortwave radiation (SRAD), carbon dioxide (CO₂) and aerosol optical depth (AOD) exhibited discrepancies in terms of degree and scope. The ARIMA model achieved satisfactory prediction performance in most areas, though the accuracy was influenced by both data values and data quality. The model can potentially be generalized for other biophysical parameters with distinct seasonality. Our findings are further verified and corroborated by four widely used GPP products, demonstrating a good consistency of GPP trends and prediction. Our analysis provides a robust framework for characterizing long-term GPP dynamics that shed light on the improved assessment of the environmental quality of terrestrial ecosystems.

1. Introduction

Gross primary production (GPP), the summation of carbon fixed by plants during photosynthesis, is of great significance in a variety of ecosystem functions such as growth and respiration [1]. As the main driving force of terrestrial carbon sequestration, GPP has a profound impact on global carbon balance and can partially regulate the atmospheric carbon dioxide (CO₂) emissions from fossil fuel consumption [2,3]. Meanwhile, GPP is a key element for ecological regulatory behaviors and an important indicator of ecosystem health [4]. It affects food, fiber, energy, and provides ecosystem services that alleviate the increasingly serious greenhouse effects [5,6].

Characterizing the spatiotemporal variations and trends of GPP is vital for a deeper understanding of the carbon cycle between terrestrial ecosystems and the atmosphere [7]. Due to the profound impact of GPP on terrestrial vegetation, its variations and trends reflect ecological environmental quality [8]. More importantly, GPP dynamics integrate geochemical, ecological and human impacts into terrestrial ecosystems. As a result, the continuous monitoring of GPP will enhance our ability to explain the complex mechanisms of the biosphere [9]. However, since current GPP products rely on assumptions and are based on different estimation models, there is not a set of GPP products that can perform better consistently in different ecosystems and external conditions. As a result, the variations and trends obtained from different GPP products are often differential [10], and thus the research scope of GPP dynamics is limited to the local scale. In previous studies, Xie et al. [11] analyzed the GPP dynamics in China’s Three-North Region over the last forty years by comparing the trend differences and the underlying causes in four sub-regions. Ge et al. [12] evaluated the spatiotemporal pattern of the phenological index and the GPP from 2001 to 2017 in the Yungui Plateau, and explored its response to droughts. The aforementioned research has improved the understanding of GPP dynamics in specific areas to varying degrees. However, the external conditions of these study areas were relatively monotonic, i.e., they were at a local scale. The climatic and hydrothermal environments of a study area which can be called large scale are generally not unified, and the research results obtained from single GPP products are also unconvincing [13]. With that being said, local-scale analysis is not conducive to understanding the macro patterns of vegetation productivity, and it is therefore urgent to investigate the variations and trends of GPP at a large scale. At the same time, in order to ensure the accuracy of the analysis results, multiple sets of GPP products should be used to verify each other.

Climate change affects ecosystem productivity by changing the absorption of light, heat, and water during plant growth and development [14,15]. The effects of climate factors on vegetation productivity have significant spatiotemporal differences [16], and the response of the productivity of different vegetation types to climate factors is also inconsistent [17]. In the past few decades, the global climate has experienced major changes represented by frequent extreme weather catastrophes and continuous warming [18]. Ecosystems and terrestrial carbon sinks will be significantly affected by climatic variation. For example, the plant productivity of Europe has decreased by about 30% because of droughts and heat waves [19]. On the other hand, GPP is substantially impacted by human activities [20]. Although a massive implementation of afforestation is conducive to the restoration of vegetation, which could increase GPP, population growth and rapid economic development have intensified urbanization and undermined the balance of terrestrial ecosystems. These human interventions have a significant impact on the formation of GPP dynamics [21,22]. Nevertheless, a thorough impact analysis of climatic variation and human interventions on GPP is necessary and still worth studying, owing to the complex interactions between the carbon cycle and environmental factors.

Since GPP provides important ecological services for human survival [23,24], grounded insights into the short-term prediction of GPP lays the basis for untangling how the terrestrial carbon cycle will evolve in the near future. The autoregressive integrated moving average (ARIMA) model is recognized as a functional statistical tool for analyzing and predicting time series data. The model is valued for its powerful performance [25], with its strength mainly reflected in the recognition and consideration of seasonality and sequence correlation in data [26]. At present, the ARIMA model has been applied in many disciplines such as aerosols [27,28], energy [29], climatology [30], and hydrology [31]. However, the currently available GPP products are inconsistent in terms of periodicity change, space distribution, and data availability [32], and such inconsistencies makes it a challenge to obtain a coherent prediction analysis. Accordingly, whether the ARIMA model can reasonably delineate and predict spatial-temporal patterns of GPP remains unknown.

In recent decades, China’s ecological construction has been increasingly conducted because of the rapid growth of the economy and population [33,34]. In view of the important role of GPP in the ecosystem, the study of GPP in China is of great significance in deepening the understanding of the terrestrial ecosystem and providing basic support for China’s ecological construction and policy-making. In light of this, our study aims to analyze the variations, trends, and impact factors of GPP dynamics over more than three decades (from 1982 to 2015) in China, and further predict short-term GPP dynamics. The specific objectives of our study are: (1) to characterize the GPP dynamics in nine major river basins in China at both seasonal and annual scales; (2) to assess the response of GPP tendencies to climatic variations and human activities; and (3) to examine the feasibility of the ARIMA model in simulating and forecasting GPP profiles and to pinpoint the potential sources of uncertainty associated with model performance.

2. Method

2.1. Study Area

China is located in the middle and low latitudes of the eastern hemisphere, with a huge land span and a vast territory. In eastern China, the continental monsoon that is dry and cold prevails in winter, while the maritime monsoon that is rainy and humid prevails in summer. As for western China, the climate is drought due to the distance from the sea [35].

Due to the complex geographical condition of China, the distribution of GPP has substantial spatiotemporal variability. To better examine such variability, we divide China into nine sub-regions based on its river systems [36] and watershed ranges. They are the Songhua and Liaohe River Basin (SRB), the Continental Basin (CB), the Haihe River Basin (HRB), the Huaihe River Basin (HuRB), the Yellow River Basin (YERB), the Yangtze River Basin (YRB), the Southwest Basin (SWB), the Southeast Basin (SEB) and the Pearl River Basin (PRB). Figure 1a shows the geographical distribution of each basin and the corresponding vegetation types.

2.2. GPP Product

As mentioned in the introduction, since large inconsistencies exist in current GPP products, five sets of GPP products were collected to ensure the robustness of our research. The available data included the optimum light use efficiency (LUEopt) product, the GPPNIRv product, the Global Land Surface Satellite (GLASS) product, the FluxCom product, and the global Orbiting Carbon Observatory-2 solar-induced chlorophyll fluorescence (GOSIF) product. The individual product information for each product used in this study can be found in Table 1, as detailed below.

The LUEopt product is obtained from The Oak Ridge National Laboratory (https://earthdata.nasa.gov/, accessed on 17 October 2021). Based on the light use efficiency (LUE) models, this product is improved by the optimum LUE originating from the global FLUXNET network [37] and is proven to obtain a better understanding of the response of primary productivity to climatic variation [38]. Zhang et al. [10] analyzed and compared 45 global GPP products and found that the estimated value of LUEopt on the seasonal scale was closer to the average value of the 45 sets of products compared with the other products we collected. Considering that analysis on a seasonal scale is an important part of this study, and LUEopt itself has relatively high spatio-temporal resolution, it was mainly used for analysis. The spatial distribution of the annual average GPP of this product during the period 1982–2015 is shown in Figure 1b.

Among the remaining GPP products for verifying the reliability of the analysis results, the GPPNIRv product was obtained from the National Tibetan Plateau Data Center (http://data.tpdc.ac.cn/en/, accessed on 17 October 2021). This dataset provides global long-term GPP based on the strong correlation between near-infrared reflectance (NIRv) and GPP, demonstrating the ability of NIRv in capturing the long-term trends of GPP [39].

The GLASS product was provided by the National Earth System Science Data Center (http://www.geodata.cn/, accessed on 25 October 2021). With the characteristics of satisfactory precision, high resolution, and a wide range of spatial-temporal coverage, this product provides powerful assistance for the study of global ecological and environmental changes [40].

The FluxCom product was downloaded from the FLUXCOM initiative (http://fluxcom.org/, accessed on 29 October 2021). It combines remote sensing and meteorological data based on a variety of machine learning algorithms, such as artificial neural network and random forest, and generates global flux tower GPP products through eddy covariance technology [41].

The GOSIF product is available from the Global Ecology Group (https://globalecology.unh.edu/data/GOSIF.html, accessed on 2 November 2021). This product is generated based on the solar-induced chlorophyll fluorescence (SIF) retrieved from the Orbiting Carbon Observatory-2 (OCO-2). SIF is regarded as an indicator of the functional state of actual plant photosynthesis [42]. The GOSIF product is suitable for studying carbon cycle and ecosystem function at various temporal and spatial scales with high spatial resolution, a good time range, and globally continuous coverage [43].

2.3. Climate and Anthropogenic Variables

The impact factors considered in this study included climatic and anthropogenic factors. Because of the impact of the hydrothermal environment and solar radiation on plant growth, the influencing factors of climate variation were selected as air temperature, precipitation, and downward shortwave radiation (SRAD). In terms of human activities, the influencing factors included carbon dioxide (CO₂), aerosol optical depth (AOD), and leaf area index (LAI). Among them, CO₂ and AOD can partly reflect industrial production activities and LAI can partly reflect the degree of vegetation restoration and protection.

Air temperature, precipitation and SRAD data were provided by the Terraclimate dataset (https://www.climatologylab.org/terraclimate.html, accessed on 4 November 2021). Using climatically aided interpolation, the Terraclimate dataset combines the WorldClim dataset with multiple sources of coarse resolution data, generating monthly climate data at a 4 km spatial resolution. These datasets provide an important basis for large-scale ecological environment research [44].

LAI data are available from the Advanced Very-High-Resolution Radiometer. The data are based on the National Oceanic and Atmospheric Administration’s climate data record (CDR) program and are produced by an artificial neural network algorithm, providing high-quality time series data for climate research [45].

CO₂ data are available from the greenhouse gases observing satellite (GOSAT) “IBUKI” (https://www.gosat.nies.go.jp/en/, accessed on 7 November 2021). Launched in 2019, the satellite is specially designed for the observation of greenhouse gases. The data obtained by the satellite were compared with the ground station data and proven to be reliable [46].

AOD data were provided by Modern-Era Retrospective analysis for Research and Applications, Version 2 (MERRA-2) (https://gmao.gsfc.nasa.gov/reanalysis/MERRA-2/, accessed on 11 November 2021). MERRA-2 has been improved on the basis of its predecessor, MERRA. It assimilates AOD through multi-source observational records, including the Moderate Resolution Imaging Spectroradiometer (MODIS), the Aerosol Robotics Network (AERONET), the Multi-angle Imaging SpectroRadiometer (MISR), and the Advanced Very-High-Resolution Radiometer (AVHRR) [47].

2.4. Trend Analysis

The significance of time series trends can be detected by parametric or non-parametric methods. Since GPP generally does not follow normal distribution, commonly used parametric methods, such as the t-test, were not applicable [48]. In light of this, we used non-parametric methods that only required the data to be independent to detect the trends of GPP.

2.4.1. Mann–Kendall Test

The Mann–Kendall test reflects the variation of the latter part of the observed value compared with the previously observed value [49]. It can be described as follows:

(1)$\begin{array}{c}S={\displaystyle {\displaystyle \sum }_{i=1}^{n-1}{\displaystyle \sum }_{j=i+1}^n}sgn\left(x_j-x_i\right) \end{array}$

where $n$ is the number of research years, $x_i$ and $x_j$ are the GPP values in year $i$ and j ($j$ > $i$), respectively, and the definition of $sgn\left(x_j-x_i\right)$ is:

(2)$\begin{array}{c}sgn\left(x_j-x_i\right)=\left\{\begin{array}{c}+1, \mathit{if} x_j-x_i>0\\ 0, \mathit{if} x_j-x_i=0\\ -1, \mathit{if} x_j-x_i<0\end{array}\right. \end{array}$

The variance $Var\left(S\right)$ is calculated as:

(3)$\begin{array}{c}Var\left(S\right)=\frac{n\left(n-1\right)\left(2n+5\right)-{{\displaystyle \sum }}_{i=1}^m{k}_i\left(k_i-1\right)\left(2k_i+5\right)}{18}\end{array}$

where $n$ is the length of years, $m$ represents the number of tied groups, and $k_i$ is the number of ties in extent $i$. The standard normal deviate ($Z$ statistics) is calculated as:

(4)$\begin{array}{c}{Z}_S=\left\{\begin{array}{c}\frac{S-1}{\sqrt{Var\left(S\right)}},\mathit{if} S>0\\ 0,\mathit{if} S=0\\ \frac{S+1}{\sqrt{Var\left(S\right)}},\mathit{if} S<0\end{array}\right.\end{array}$

A positive value of $Z_S$ means an upward trend while a negative $Z_S$ value indicates a downward trend. We tested the trend at the significance α level. If |$Z_S$| > $Z_{1-\mathsf{\alpha }/2}$, it means rejecting the null hypothesis at the α level and that there is a significant trend in the GPP time series.

2.4.2. Sen’s Slope Estimator

Sen proposed a non-parametric method for calculating the slope of trend in time series data with the sample of $N$ pairs [50]:

(5)$\begin{array}{c}{Q}_i=\frac{x_j-x_k}{j-k} for i=1,\dots ,N \end{array}$

where $x_j$ and $x_j$ represent the GPP values at year $j$ and $k$ ($j$ > $k$), respectively. If each period has only one datum, then $N=n\left(n-1\right)/2$, where $n$ represents the number of data points. If there are multiple values in periods, then $N<n\left(n-1\right)/2$, where $n$ represents the total number of values.

Values of $Q_i$ are ranked from large to small, and the median of the slope estimator is calculated as:

(6)$\begin{array}{c}{Q}_{med}=\left\{\begin{array}{c}{Q}_{\left[\frac{N+1}{2}\right]} ,\mathit{if} N \mathit{is} odd\\ \frac{Q_{\left[\frac{N}{2}\right]}+Q_{\left[\frac{N+2}{2}\right]}}{2},\mathit{if} N \mathit{is} even\end{array}\right.\end{array}$

The $Q_{med}$ sign reflects the variation of data, and its value shows the gradient of the trend. We need to get the confidence interval of $Q_{med}$ at the α level to judge whether the values of $Q_{med}$ are statistically different from zero.

The confidence interval of the slope of time series can be calculated as:

(7)$\begin{array}{c}{C}_{\alpha }=Z_{1-\frac{\alpha }{2}}\sqrt{Var\left(S\right)}\end{array}$

We define that $M_1=\left(N-C_{\alpha }\right)/2$ and $M_2=\left(N+C_{\alpha }\right)/2$. $Q_{max}$ and $Q_{min}$ represent the range of the confidence interval corresponding to the ($M_2+1$)th largest and the $M_1$th largest values in the $N$ order slope estimation, respectively. The slope $Q_{med}$ is different than zero statistically if $Q_{min}$ and $Q_{max}$ have the same sign.

2.5. Correlation Analysis

The GPP is a coupling process of multiple influencing factors in its growing environment. To analyze the relationship between GPP and these influencing factors, the correlation coefficients between GPP and air temperature, precipitation, SRAD, LAI, CO₂, and AOD are calculated separately.

The correlation value is [–1, 1], and the calculation formula is as follows:

(8)$\begin{array}{c}{R}_{xy}=\frac{{{\displaystyle \sum }}_{i=1}^n\left(x_i-\overline{x}\right)\left(y_i-\overline{y}\right)}{\sqrt{{{\displaystyle \sum }}_{i=1}^n{\left(x_i-\overline{x}\right)}^2}\sqrt{{{\displaystyle \sum }}_{i=1}^n{\left(y_i-\overline{y}\right)}^2}}\end{array}$

where $R_{xy}$ is the correlation coefficient between factor $x$ and $y$, indicating the correlation degree between the two factors; $n$ is the number of research years; $x_i$ is the GPP of vegetation in year $i$; and $y_i$ is the mean value of correlation factor in year $i$. $R_{xy}$ > 0 indicates a positive correlation whereas $R_{xy}$ < 0 indicates a negative correlation, and $\overline{x}$ and $\overline{y}$ are the average values of the sample values of the two elements, respectively. The significance of the correlation is further tested afterwards.

2.6. Autoregressive Integrated Moving Average (ARIMA) Model

The time series ARIMA model is used for analysis and prediction. The time series data formed by samples changing with time are regarded as a random series. Its profile is simulated by fitting a best-fit ARIMA model according to the past and current values of the time series, and the built model is used to forecast future GPP values. Broadly, the establishment of the ARIMA model mainly has three processes: model identification, parameter and diagnostic checking, and prediction. After completing the three processes, we carefully selected the most suitable model from the above stages to predict the future values. Figure 2 shows the schematic workflow of the ARIMA model.

2.6.1. Identification

The first step of the ARIMA model is to estimate the stationarity of the time series, which can be examined by the autocorrelation function and partial autocorrelation function. The attenuation rate of the autocorrelation function plot reflects whether the data are stable or not. Non-stationary time series data sets require a different process to make them stationary. The ARIMA model can be defined by three operators (p, d, q), where p, d, q denotes the orders of the auto regressive (AR) process, the differencing process, and the moving average (MA) process, respectively, and is given by:

(9)$\begin{array}{c}{\varnothing }_p\left(B\right){\left(1-B\right)}^d{X}_t={\theta }_q\left(B\right){\epsilon }_{t }\end{array}$

where $X_t$ is the GPP at time point $t$; $X_t$ is the backshift operator that follows $B{X}_t=X_{t-1}$, $B^2{X}_t=BB{X}_t=B{X}_{t-1}=X_{t-2}$, …, $B^k{X}_t=X_{t-k}$; ${\mathsf{\varnothing }}_p\left(B\right)$ is the AR operator and ${\mathsf{\varnothing }}_p\left(B\right)=1-{\mathsf{\varnothing }}_{1}B-{\mathsf{\varnothing }}_2{B}^2-\dots -{\mathsf{\varnothing }}_p{B}^p$; ${\theta }_q\left(B\right)$ is the MA operator and ${\theta }_q\left(B\right)=1+{\theta }_{1}B+{\theta }_2{B}^2+\dots +{\theta }_q{B}^q$; and ${\epsilon }_t$ represents the random errors that resemble a white noise process ${\epsilon }_t~WN\left(0,{\sigma }^2\right)$.

2.6.2. Parameter and Diagnostic Checking

After passing the stationarity test of time series data, we fit the ARIMA model by estimating the values of p and q. For a model that can adequately describe the trends of the GPP sequence, the distribution of the residuals should be similar to white noise. Meanwhile, the goodness of fit statistics of the model should be the most accurate among the candidate models. Based on this, we considered the following test statistics: mean absolute error (MAE), root mean squared error (RMSE), relative mean absolute error (rMAE), and relative root mean square error (rRMSE), the specific definitions of these statistics can be found in the Appendix A. Taken together, we used these metrics to assess the fitting effect of the ARIMA model on the GPP series.

2.6.3. Prediction and Evaluation

The last procedure is to establish a prediction model. Specifically, the model that provides the best diagnostic statistics (the values of MAE, RMSE, rMAE, and rRMSE are the lowest) is used to predict the GPP values of the next 12 months. In our study, we divided China into 410 × 684 pixels, and each pixel of GPP data was regarded as a separate model. The above three stages were used to iteratively fit the model to the individual grid cells and finally generate the national prediction map.

2.7. Analysis Framework

In this study, the trends of GPP at annual and seasonal scales were first analyzed in China during 1982–2015 by using the Mann–Kendall test and Sen’s method. The Mann–Kendall test was used to evaluate the reliability of the results, and Sen’s slope method was used to quantify the trends of GPP. We reorganized the original GPP products into seasonal time series data and annual time series data, which were used as the input of the Mann–Kendall test and Sen’s method, respectively. We divided the twelve months into four seasons according to the following definitions: March–April–May representing spring; June–July–August representing summer; September–October–November representing autumn; and December–January–February representing winter.

The Pearson correlation coefficient was calculated to ascertain the correlation degree between GPP and the impacting factor. We calculated the annual mean of each impact factor and then computed its correlation with the annual mean of GPP. Potential causes that accounted for the trend results were proposed from the perspectives of climatic variation and human activities.

The ARIMA model was applied to forecast the monthly GPP in 2016. Taking RSME, MAE, rRMSE and rMAE as indicators, the predicted values were compared against the reference values. We also assessed the model’s uncertainty from the perspectives of data range and data availability.

To investigate the robustness of our framework, we further compared the variations, trends, and prediction results from LUEopt GPP against those from the other three products, i.e., the FluxCom product, the GPPNIRv product, and the GLASS product. The GOSIF product was not involved in the trend analysis because it did not contain the data before 2001. However, all sets of the products were used for the input of the ARIMA model. Similarly, we took RMSE and MAE as indicators to assess the reliability of the ARIMA prediction results.

3. Result

3.1. Annual and Seasonal Trends of GPP Calculated from LUEopt Products

Figure 3 exhibits the spatial distribution of the annual and seasonal trends of GPP nationwide while Figure 4 shows the GPP trends at the basin scale. During the last thirty-four years, the GPP of terrestrial vegetation in China has shown an increasing trend, consistent with the studies of Anav et al. [32] that focused on multiple sets of GPP data. This can be explained by the effect of CO₂ concentration and air temperature rise on vegetation production in recent years [51]. On the other hand, in the 21st century, the Chinese government has issued a series of policies to restore the ecological environment, such as the Grain to Green Program [52,53]. These policies may facilitate an increase in GPP [54,55]. However, specific to basin scale, certain differences are detected in the trends of varying river basins at different time scales.

Annually, except for the SEB which showed a decreasing trend of −0.19 g C/m²/year), the remaining basins all exhibited an increasing trend. Among them, the HRB showed the largest average ascending trend (3.43 g C/m²/year), and the YRB ranked the second highest (2.67 g C/m²/year), followed by the SRB, the HuRB, the YERB, and the SWB, with the four being nearly equal (2.19, 2.07, 2.06 and 1.82 g C/m²/year). Because of the large variability within pixels (the standard deviation is 5.48 g C/m²/year), the PRB showed an increasing trend of only 1.05 g C/m²/year. With an inherently small mean GPP value, the CB only showed a slight increase (0.24 g C/m²/year). Relatively significant upward trends were mainly distributed in the HRB, the YERB, and the border areas of the YRB, the SWB, and the PRB. Meanwhile, the significant downward trends were mainly distributed in the eastern region of the PRB. Our results are consistent with the conclusions of Ge et al. [12].

Seasonally, the four northern basins (SRB, CB, HRB and YERB) revealed an increasing trend in spring, summer and autumn, despite a stagnant period in winter, which was probably because of the low air temperature and SRAD in winter that leads to the stagnation of plant photosynthesis [56]. In the HuRB and the YRB, which are located in the central part of China, a significant increasing trend was detected in spring. Although there was a downward trend in summer, the ascending trends in the other seasons rendered the annual trend of the two basins to be rising. In the SEB and the PRB, statistically significant decreasing trends were found in summer. This can be attributed to the fact that high air temperature is more likely to occur in summer, which inhibits vegetation growth and is not conducive to increases in GPP [57]. Due to the increasing trend of GPP in spring and winter, the annual trend of the PRB exhibited an increasing state. The SEB was the only basin with a downward annual trend because of the decreasing trend in summer (−2.56 g C/m²/year). Nevertheless, the SWB maintained an upward trend consistent with the other three seasons in summer when the other basins in the south of China generally showed a downward trend. This can be explained by the fact that the SWB area has a high altitude and the terrain is mainly mountainous, which is not favorable to having extremely high air temperatures in summer [58]. On the whole, the noticeable discrepancies of the GPP trends in the nine river basins mainly occurred in summer, which mostly dominated the annual trend pattern of the basin.

3.2. Attributions of Long-Term Changes in GPP

Figure 5 displays the spatial pattern of primary impacting variables and Figure 6 shows the Pearson correlation coefficient between the individual factor and GPP annual variation. At a national scale, the correlation coefficients (R) were 0.318 (LAI), 0.225 (air temperature), 0.16 (CO₂), 0.053 (precipitation), −0.012 (SRAD) and −0.056 (AOD), respectively. At the basin scale, different factors had obvious regional discrepancies in the impacts of the annual variation of GPP [9,59]; these regional discrepancies were mainly reflected in the differences in hydrothermal conditions and the strength of human activities. The average value and standard deviation of the Pearson correlation coefficient in each basin are summarized in Table 2.

The hydrothermal conditions of water and air temperature and their interaction are generally considered the main climate factors affecting vegetation growth and distribution [60]. In our study, GPP and air temperature did show a relatively positive correlation in most parts of China, except the PRB and the SEB, which can be related to the fact that high air temperature in summer is beyond the optimal growth range of plants, leading to an increase in evapotranspiration and a decrease in soil water availability for vegetation, and thus inhibiting the growth of plants [57,59]. Notably, the positive correlation between GPP and air temperature in the SRB was not as substantial as in the other basins, probably owing to the reduction in precipitation that weakens the positive effect of warming on GPP [61]. A discrepancy in the GPP–precipitation correlation between the north and south regions was observed (Figure 6c), i.e., the correlation decreased from north to south. Areas with positive correlation were mainly distributed in the CB, the HRB and the YERB, which possess relatively little precipitation. This implies that, in arid areas, the increase in precipitation had a positive effect on GPP, while the sensitivity of GPP to precipitation was less prominent in humid areas. Unlike the driving force of photosynthesis [62], our study did not show a strong positive correlation between SRAD and GPP, and some basins in the north, such as the CB and the HRB, even showed a negative correlation (R < −0.2). One explanation for this phenomenon is that strong solar radiation often reflects a lack of cloud cover, resulting in a reduction in precipitation and the enhancement of soil water evaporation [63]. As mentioned above, in arid regions such as northern China, vegetation is very sensitive to water changes, and the promotion of photosynthesis by radiation is not sufficient to offset the limitations brought by drought.

Human activities can have substantial effects on vegetation productivity. Influenced by the policies of the Chinese government, for example, during 1999–2010, 16,000 square kilometers of farmland were turned into forest or grass by the Grain to Green Program [64]. LAI in most parts of China have shown an upward trend, which has a significant correlation with the trends of GPP. In the YERB, the northern HRB, the western SRB, and the eastern CB, the correlation was particularly evident (R > 0.3), implying that the government’s efforts in improving the environmental quality of ecosystems have been effective.

Another factor regarding human activities is the increase in CO₂ emissions. For example, in 2011, the use of fossil fuels related to human activities resulted in 33.2 billion tons of CO₂ emissions worldwide [65]. The increase in atmospheric CO₂ concentration can improve vegetation productivity. Devaraju et al. [66] found that global primary productivity increases by about 2.3 Pg C each year. Schimel et al. [6] reported that about 60% of the carbon absorbed by land comes from CO₂ emissions. Our results show that the rise in CO₂ emissions in most parts of China can promote the growth of GPP. However, the correlation between CO₂ and GPP is not significant in the SWB, the southern CB, the western YRB, or the northern SRB. Among them, the SWB was even negatively correlated (R = −0.101). This may have been due to the air temperature, which is low in these areas, and the GPP dynamic is more affected by hydrothermal conditions rather than CO₂.

Furthermore, due to the human activities caused by rapid urbanization and industrial expansion, the air pollution level of more than 90% of China’s cities exceeds the standards of the World Health Organization [67]. Since the amplitude of aerosols reflects the degree of air pollution to a certain extent, AOD is used to quantify aerosols [68,69]. In this study, we analyzed the correlation between AOD and GPP. A negative correlation between AOD and GPP was observed on a national scale. This is acceptable, since the absorption of atmospheric aerosols reduces the photosynthetically active radiation in plants [70,71], which affects vegetation productivity. However, this correlation was weak, regardless of whether it was positive or negative. For example, even where there was a positive correlation in some areas, such as the HRB, the correlation was still very weak. This is probably because factors such as air temperature and precipitation often have a more obvious impact on vegetation growth, which dampens the response of vegetation productivity to AOD. Another possible explanation is the complex climatic effects of aerosol species. The net effect might be small due to their counteract responses to GPP. In future work, the multivariate linear regression that can be used to control other factors will be considered for critical analysis.

3.3. Short-Term Prediction of GPP by ARIMA Model

Time series analysis provides an effective tool for GPP modeling. In this study, the ARIMA model was used to simulate and predict GPP based on the monthly dataset during 1982–2015. Figure 7 shows the comparison between the annual average of reference GPP and the ARIMA-predicted GPP in 2016. The spatial pattern of the predicted value is basically consistent with that of the reference values, especially in the northern regions such as the SRB, the HRB and the YERB. Areas with a large degree of uncertainty were mainly located in the YRB and the SEB. Overall, the model has a strong capacity for predicting short-term GPP time series.

Figure 8 exhibits the accuracy of the prediction results. In general, the ARIMA model achieved good results in most parts of China. The national averaged values of RMSE, MAE, rRMSE, and rMAE were 0.571 g C/m²/d, 0.401 g C/m²/d, 24.8%, and 16.5%, respectively. Regionally, the ARIMA model showed relatively stable results for the northern basins such as the SRB, the HRB, the CB and the YERB (the averaged values of RMSE and MAE in these four basins were 0.383 g C/m²/d and 0.248 g C/m²/d, respectively), whereas its performance for southern basins with high GPP values such as the YRB and the PRB tended to be more volatile (the average values of RMSE and MAE in these two basins were 0.822 g C/m²/d and 0.606 g C/m²/d, respectively). This is probably related to the high correlation between low accuracies and high GPP values.

Earlier studies have demonstrated that the availability of data greatly affects the prediction ability of the ARIMA model [72]. The uncertainty of the ARIMA model regarding the data value range and data availability is shown in Figure 9. It is observed that pixels with poor availability are mainly distributed in the CB and the SWB, and the GPP values of these two regions are relatively low. The results indeed show that the model capability is closely related to data availability. The more available the data are, the better the forecast result will be, and vice versa. Moreover, the ability of the ARIMA model is noticeably affected by GPP values (Figure 9a,b). The curves of Figure 9c,d both show an upward trend in the availability of 100%. This is because the pixels with availability below 100% have lower GPP values than those with almost 100% data coverage, leading to their comparatively high accuracy. In other words, despite some places with poor data quality, the precision of GPP prediction is relatively high because of the low GPP values. All in all, the precision of the model is affected by the value and quality of the data together [73].

4. Discussion

4.1. Comparison of Multi-Source Data

Considering the inconsistencies in the values of current GPP products among different regions, seasons and vegetation types [10], multiple sets of GPP products were selected to further corroborate our results. It should be mentioned that our study merely focuses on consistencies and discrepancies in GPP variation, trends and predication, rather the intrinsic difference between the GPP products.

Figure 10 shows the temporal and spatial comparison of the monthly mean values of different GPP products. It is observed that five sets of GPP products are consistent in seasonality, i.e., the value is smaller in winter while the value reaches the peak in summer. Spatially, the coincident patterns occur in the SRB, the HRB and the YRB regions. Remarkable differences are found in the northwest inland region represented by the CB and the southeast coastal region represented by the SEB. The differences in the CB may be related to the lesser vegetation coverage and the inconsistencies of the GPP anomalies among different products. As for the SEB, the discrepancy in GPP values can be explained by the fact that different products are discordant in the quantification of the impact of the hydrothermal environment on GPP.

The trends of GPP in each basin obtained by different products are shown in Figure 11. The overall trends calculated by each set of products are consistent, which demonstrates the growing trend of GPP in China. However, there are certain differences among basins. Such differences mainly appear in the SEB and the PRB where the GPPNIRv product shows a positive trend (0.112 g C/m²/year and 0.447 g C/m²/year), and the GLASS product a negative one (−0.338 g C/m²/year and −0.028 g C/m²/year). The LUEopt product shows a negative trend in the SEB (−0.186 g C/m²/year), and a positive trend in the PRB (1.051 g C/m²/year). Compared with the above products, the results generated by FluxCom products were more conservative, i.e., the maximum value occurred in the SEB (0.07 g C/m²/year) and the minimum value happened in the HuRB (−0.014 g C/m²/year).

Figure 12 shows the accuracy of the predicted GPP using ARIMA obtained by the four products. The prediction results of each set of products have good accuracy. It is observed that the FluxCom product showed better accuracy, which may be related to its low GPP value. The GPPNIRv product had the most variable accuracy. Consistent with the above conclusions that the performance of the ARIMA model is affected by the GPP values, the predicted value of the model had high uncertainty in areas with high GPP values, such as the PRB (the averaged RMSE and MAE of the four sets of products were 0.36 g C/m²/d and 0.289 g C/m²/d, respectively). In contrast, in areas with low GPP values such as the CB, the predicted value of the model had better accuracy (the averaged RMSE and MAE of the four sets of products were 0.03 g C/m²/d and 0.025 g C/m²/d, respectively).

4.2. Implications for Other Biophysical Studies

Our previous work proposed the use of the ARIMA model in estimating atmospheric parameters such as aerosols at local and global scales [27,28,72]. In a broader sense, we are motivated to apply the time series model to biophysical parameters such as GPP. The reasonable forecasts of the present study demonstrate the generalization of the model, and that it was able to capture the trends and periodicity inherent in the time series. We hope that the stochastic ARIMA model will benefit the simulation and prediction of biophysical parameters that are characterized by strong seasonality, such as those related to the vegetation patterns and dynamics. To make it more appropriate for planning and regulation purposes, as well as suitable for reducing uncertainties in global vegetation models [32], our future work plans to capitalize on cutting-edge computing platforms such as the Google Earth Engine (GEE), which enables automatic and mass production.

On the other hand, it should also be noted that this trend analysis of GPP was mainly based on long-term variations of time series, and the ability to capture short-term mutations caused by extreme phenomena was insufficient. Gampe et al. [74] identified local extremes at the grid level and joined them with a flood-filling method to capture extreme GPP events, which provides a good strategy for our further work. In addition, the correlation analysis between impacting factors and GPP was conducted only at an annual scale, but not a seasonal scale, and there was no specific quantification of the relative contribution of each factor. Another deficiency in the correlation analysis was that, subject to the problem of spatio-temporal resolution, we selected LAI, AOD, and CO₂ as the influencing factors to reflect human activities. However, these data can only partially reflect human activities. How to quantify the importance and contribution of each impact factor seasonally and how to extract information that can fully reflect human activities will be another focus of our future work. Moreover, in terms of the discrepancies in the results obtained by different GPP products, an exploration of the reasons for these inconsistencies could be given more attention. The production principles and input parameters for each set of data should be further considered to clarify the source of these discrepancies.

5. Conclusions

This study used the Mann–Kendall and Sen’s methods to analyze the trends of GPP dynamics in China from 1982 to 2015 and explore the impacting factors behind these variations. Our aim was to contribute to understanding the response of vegetation growth to climatic variation and human activities at a national scale. On this foundation, we further used the ARIMA model to forecast the one-year lead GPP values in China and assess the potential and uncertainty of the model. To corroborate the reliability of our results, we also extended our trend and prediction analyses into multiple sets of GPP products. Our main findings were that:

(1) Over the whole of China, an upward GPP trend was observed with an average trend of 2.268 g C/m²/year during 1982–2015. However, the temporal pattern of GPP trends showed certain heterogeneity, and a more obvious downward trend occurs in the summer. Spatially, the GPP trends in most regions were mainly increasing, but in the SEB and the PRB the trends calculated by five sets of products were not consistent;

(2) Climate factors and human interventions affect long-term GPP jointly. On a national scale, air temperature is a major climatic factor impacting GPP and the rise in LAI is a major human factor. The effect of precipitation on GPP is mainly reflected in the arid areas in the north and the effect of CO₂ on GPP is reflected in areas with better hydrothermal conditions in the east. AOD and SRAD also affect GPP, but not as obviously as the above factors;

(3) The ARIMA model performs well in most parts of China, as evidenced by the relatively low values of MAE, RMSE, rMAE and rRMSE, as well as the good consistency with the reference GPP from 2016. Moreover, the performance of the ARIMA model was affected by the values and quality of the GPP data, indicating that the low value and high quality will produce more accurate predictions;

(4) Further multiple GPP analysis showed consistent seasonality and trends among four sets of products that contained the data from 1982 to 2015. However, there were certain differences in spatial distribution, especially in the CB, the SEB, and the PRB. Using the ARIMA model, five GPP products all achieved satisfactory prediction results which were consistent in space. The comparison of multi-source products showed certain discrepancies in GPP trends and predictions at the basin scale but, on the whole, the conclusions were consistent.

Author Contributions

Conceptualization, Y.B. and K.L.; methodology, Y.B., S.W. and X.L.; formal analysis, X.L. and X.Z.; data curation, K.L., X.Z. and H.Z.; writing—original draft preparation, Y.B.; writing—review and editing, Y.B., K.L., X.L. and X.G.; visualization, X.Z. and H.Z.; supervision, X.L., S.W. and X.G. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

The data presented in this study are available on request from the corresponding author.

Acknowledgments

We acknowledge the providers of all the data used in the study.

Conflicts of Interest

We declare that they have no conflict of interest to this work. Additionally, we have no financial and personal relationships with other people or organizations that can inappropriately influence our work.

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures and Tables

Figure 1. The spatial distribution of (a) vegetation type and (b) annual mean of GPP from 1982 to 2015 in each river basin. Vegetation-type data are available from NASA LP DAAC (https://lpdaac.usgs.gov/, accessed on 11 October 2021). The box-plot in (b) represents the interquartile range (box), median (line within the box), and whiskers of the annual mean GPP in different basins. The nan value is displayed in white.

Figure 2. Schematic map of the autoregressive integrated moving average (ARIMA) model. Each pixel corresponds to an independent modeling process. The diagram depicts the actual time series and simulated time series obtained from the ARIMA model.

Figure 3. Spatial distribution of GPP trends at annual and seasonal scales from 1982 to 2015. The values in the figure represent the change in GPP cumulative value in a year (or a season) compared with the last year (or the last same season) in the unit area. Values with statistical significance at a 95% confidence level are displayed with cross symbols. The nan value is displayed in white.

Figure 4. Annual and seasonal GPP trends from 1982 to 2015 in each basin. The error bars reflect the standard deviation of the trends.

Figure 5. Annual mean of air temperature, downward shortwave radiation (SRAD), precipitation, leaf area index (LAI), and aerosol optical depth (AOD) from 1982 to 2015, and CO2 emissions from 2010 to 2015. The definition of box and whiskers on the lower right of each panel is the same as Figure 1.

Figure 6. Spatial distribution of the relevance of (a) air temperature, (b) SRAD, (c) precipitation, (d) LAI, (e) AOD, and (f) CO2 to GPP. The correlation map is based on the pixel-wise annual mean of GPP and each factor. Values with statistical significance at the 95% confidence level are displayed with cross symbols. The nan value is displayed in white.

Figure 7. The annual mean of the (a) reference GPP value and (b) predicted GPP value from the ARIMA model in 2016. The nan value is displayed in white.

Figure 8. Spatial patterns of RMSE, MAE, rRMSE, and rMAE computed by the 12-month reference GPP values and predicted values using the ARIMA model in 2016. The nan value is displayed in white.

Figure 9. Uncertainty analysis of the ARIMA model regarding (a,b) data value range and (c,d) data availability, assessed with RMSE and MAE. The shaded areas represent one standard deviation around the mean RMSE or MAE.

Figure 10. Comparison of monthly mean values from different GPP products.

Figure 11. The trend of GPP in each basin obtained by different products from 1982 to 2015. The error bars reflect the standard deviation.

Figure 12. Accuracy of predicted GPP using ARIMA obtained by the FluxCom product, the GPPNIRv product, and the GLASS product from 1982 to 2015 and the GOSIF product from 2001 to 2015.

Table 1

Descriptions of multi-source GPP products and impact factors used in the study.

ID	Variables	Source	Time Period	TemporalResolution	SpatialResolution	Units
1	GPP	Oak Ridge National Laboratory	1982–2016	Monthly	8 × 8 km	gC/m2/day
2	GPP	National Tibetan Plateau Data Center	1982–2016	Monthly	0.05 × 0.05°	gC/m2/day
3	GPP	National Earth System Science Data Center	1982–2016	8 days	0.05 × 0.05°	gC/m2/8 days
4	GPP	FluxCom	1982–2016	Monthly	0.5 × 0.5°	gC/m2/day
5	GPP	Global Ecology Group	2001–2016	8 days	0.05 × 0.05°	gC/m2/month
6	Airtemperature	Terraclimate	1982–2015	Monthly	4 × 4 km	Degree (°)/month
7	Precipitation	Terraclimate	1982–2015	Monthly	4 × 4 km	mm/month
8	SRAD	Terraclimate	1980–2015	Monthly	4 × 4 km	W/m2/month
9	LAI	AVHRR	1982–2015	8 days	0.05 × 0.05°	m2/m2
10	CO2	GOSAT	2010–2015	Monthly	2.5 × 2.5°	ppm/month
11	AOD	MERRA-2	1982–2015	Monthly	0.625 × 0.50°	\

Table 2

Area-averaged (mean) correlation coefficient and its standard deviation (STD) in each basin.

Region	Air Temperature		SRAD		Precipitation		LAI		AOD		CO2
	Mean	STD	Mean	STD	Mean	STD	Mean	STD	Mean	STD	Mean	STD
SRB	0.145	0.212	0.065	0.305	−0.03	0.343	0.4	0.217	−0.126	0.254	0.229	0.479
CB	0.265	0.266	−0.327	0.208	0.263	0.222	0.307	0.249	−0.132	0.198	0.013	0.484
HRB	0.206	0.155	−0.241	0.135	0.372	0.172	0.501	0.198	0.101	0.186	0.204	0.385
YERB	0.252	0.165	−0.201	0.208	0.224	0.189	0.427	0.239	0.008	0.211	0.268	0.431
HuRB	0.225	0.164	−0.029	0.218	0.143	0.215	0.364	0.187	0.052	0.216	0.251	0.34
YRB	0.291	0.214	0.178	0.269	−0.03	0.228	0.271	0.225	0.021	0.243	0.21	0.483
SWB	0.359	0.194	−0.003	0.24	−0.038	0.22	0.228	0.214	−0.204	0.155	−0.101	0.438
SEB	−0.002	0.196	0.231	0.168	−0.194	0.159	0.057	0.195	−0.072	0.212	0.174	0.359
PRB	0.039	0.208	0.134	0.146	−0.186	0.16	0.152	0.272	−0.025	0.311	0.176	0.432
China	0.225	0.229	−0.012	0.307	0.053	0.29	0.318	0.248	−0.056	0.246	0.16	0.47

Appendix A

Statistical metrics used in our study are:(A1)$$\begin{array}{c}RMSE=\sqrt{\frac{1}{n}{\displaystyle {\displaystyle \sum }_{i=1}^n}{\left(GP{P}_{\left(pre\right)i}-GP{P}_{\left(ref\right)i}\right)}^2}\end{array}$$ (A2)$$\begin{array}{c}rRMSE=\frac{RMSE}{GP{P}_{\left(ref\right)max}-GP{P}_{\left(ref\right)min}}\end{array}$$ (A3)$$\begin{array}{c}MAE=\frac{1}{n}{\displaystyle {\displaystyle \sum }_{i=1}^n}\left|GP{P}_{\left(pre\right)i}-GP{P}_{\left(ref\right)i}\right|\end{array}$$ (A4)$$\begin{array}{c}rMAE=\frac{MAE}{GP{P}_{\left(ref\right)max}-GP{P}_{\left(ref\right)min}}\end{array}$$ where $GP{P}_{\left(pre\right)}$ is predicted GPP values from the ARIMA model, $GP{P}_{\left(ref\right)}$ is reference GPP values from different products; $GP{P}_{\left(ref\right)max}$ represents the maximum values of reference GPP, while $GP{P}_{\left(ref\right)min}$ represents the minimum; $n$ is the number of research months.

References

1. Beer, C.; Reichstein, M.; Tomelleri, E.; Ciais, P.; Jung, M.; Carvalhais, N.; Rodenbeck, C.; Arain, M.A.; Baldocchi, D.; Bonan, G.B. et al. Terrestrial gross carbon dioxide uptake: Global distribution and covariation with climate. Science; 2010; 329, pp. 834-838. [DOI: https://dx.doi.org/10.1126/science.1184984] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/20603496]

2. Janssens, I.A.; Freibauer, A.; Ciais, P.; Smith, P.; Nabuurs, G.J.; Folberth, G.; Schlamadinger, B.; Hutjes, R.W.; Ceulemans, R.; Schulze, E.D. et al. Europe’s terrestrial biosphere absorbs 7 to 12% of European anthropogenic CO₂ emissions. Science; 2003; 300, pp. 1538-1542. [DOI: https://dx.doi.org/10.1126/science.1083592] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/12764201]

3. Battin, T.J.; Luyssaert, S.; Kaplan, L.A.; Aufdenkampe, A.K.; Richter, A.; Tranvik, L.J. The boundless carbon cycle. Nat. Geosci.; 2009; 2, pp. 598-600. [DOI: https://dx.doi.org/10.1038/ngeo618]

4. Gao, Q.; Li, Y.; Wan, Y.; Qin, X.; Jiangcun, W.; Liu, Y. Dynamics of alpine grassland NPP and its response to climate change in Northern Tibet. Clim. Change; 2009; 97, pp. 515-528. [DOI: https://dx.doi.org/10.1007/s10584-009-9617-z]

5. Norby, R.J.; Warren, J.M.; Iversen, C.M.; Medlyn, B.E.; McMurtrie, R.E. CO₂ enhancement of forest productivity constrained by limited nitrogen availability. Proc. Natl. Acad. Sci. USA; 2010; 107, pp. 19368-19373. [DOI: https://dx.doi.org/10.1073/pnas.1006463107]

6. Schimel, D.; Stephens, B.B.; Fisher, J.B. Effect of increasing CO₂ on the terrestrial carbon cycle. Proc. Natl. Acad. Sci. USA; 2015; 112, pp. 436-441. [DOI: https://dx.doi.org/10.1073/pnas.1407302112]

7. Cheng, Y.-B.; Zhang, Q.; Lyapustin, A.I.; Wang, Y.; Middleton, E.M. Impacts of light use efficiency and fPAR parameterization on gross primary production modeling. Agric. For. Meteorol.; 2014; 189–190, pp. 187-197. [DOI: https://dx.doi.org/10.1016/j.agrformet.2014.01.006]

8. Zhao, M.; Running, S.W. Drought-induced reduction in global terrestrial net primary production from 2000 through 2009. Science; 2010; 329, pp. 940-943. [DOI: https://dx.doi.org/10.1126/science.1192666]

9. Nemani, R.R.; Keeling, C.D.; Hashimoto, H.; Jolly, W.M.; Piper, S.C.; Tucker, C.J.; Myneni, R.B.; Running, S.W. Climate-driven increases in global terrestrial net primary production from 1982 to 1999. Science; 2003; 300, pp. 1560-1563. [DOI: https://dx.doi.org/10.1126/science.1082750]

10. Zhang, Y.; Ye, A. Would the obtainable gross primary productivity (GPP) products stand up? A critical assessment of 45 global GPP products. Sci. Total Environ.; 2021; 783, 146965. [DOI: https://dx.doi.org/10.1016/j.scitotenv.2021.146965]

11. Xie, S.; Mo, X.; Hu, S.; Liu, S. Contributions of climate change, elevated atmospheric CO₂ and human activities to ET and GPP trends in the Three-North Region of China. Agric. For. Meteorol.; 2020; 295, 108183. [DOI: https://dx.doi.org/10.1016/j.agrformet.2020.108183]

12. Ge, W.; Han, J.; Zhang, D.; Wang, F. Divergent impacts of droughts on vegetation phenology and productivity in the Yungui Plateau, southwest China. Ecol. Indic.; 2021; 127, 107743. [DOI: https://dx.doi.org/10.1016/j.ecolind.2021.107743]

13. Song, L.; Li, Y.; Ren, Y.; Wu, X.; Guo, B.; Tang, X.; Shi, W.; Ma, M.; Han, X.; Zhao, L. Divergent vegetation responses to extreme spring and summer droughts in Southwestern China. Agric. For. Meteorol.; 2019; 279, 107703. [DOI: https://dx.doi.org/10.1016/j.agrformet.2019.107703]

14. Piao, S.; Ciais, P.; Lomas, M.; Beer, C.; Liu, H.; Fang, J.; Friedlingstein, P.; Huang, Y.; Muraoka, H.; Son, Y. et al. Contribution of climate change and rising CO₂ to terrestrial carbon balance in East Asia: A multi-model analysis. Glob. Planet. Chang.; 2011; 75, pp. 133-142. [DOI: https://dx.doi.org/10.1016/j.gloplacha.2010.10.014]

15. Liu, K.; Li, X.; Long, X. Trends in groundwater changes driven by precipitation and anthropogenic activities on the southeast side of the Hu Line. Environ. Res. Lett.; 2021; 16, 094032. [DOI: https://dx.doi.org/10.1088/1748-9326/ac1ed8]

16. Zhu, W.; Pan, Y.; Yang, X.; Song, G. Comprehensive analysis of the impact of climatic changes on Chinese terrestrial net primary productivity. Chin. Sci. Bull.; 2007; 52, pp. 3253-3260. [DOI: https://dx.doi.org/10.1007/s11434-007-0521-5]

17. Khalifa, M.; Elagib, N.A.; Ribbe, L.; Schneider, K. Spatio-temporal variations in climate, primary productivity and efficiency of water and carbon use of the land cover types in Sudan and Ethiopia. Sci. Total Environ.; 2018; 624, pp. 790-806. [DOI: https://dx.doi.org/10.1016/j.scitotenv.2017.12.090]

18. Mirza, M. Climate change and extreme weather events: Can developing countries adapt?. Clim. Policy; 2003; 3, pp. 233-248. [DOI: https://dx.doi.org/10.3763/cpol.2003.0330]

19. Ciais, P.; Reichstein, M.; Viovy, N.; Granier, A.; Ogee, J.; Allard, V.; Aubinet, M.; Buchmann, N.; Bernhofer, C.; Carrara, A. et al. Europe-wide reduction in primary productivity caused by the heat and drought in 2003. Nature; 2005; 437, pp. 529-533. [DOI: https://dx.doi.org/10.1038/nature03972]

20. Zhang, Y.; Zhang, C.; Wang, Z.; Chen, Y.; Gang, C.; An, R.; Li, J. Vegetation dynamics and its driving forces from climate change and human activities in the Three-River Source Region, China from 1982 to 2012. Sci. Total Environ.; 2016; 563–564, pp. 210-220. [DOI: https://dx.doi.org/10.1016/j.scitotenv.2016.03.223]

21. Luo, L.; Ma, W.; Zhuang, Y.; Zhang, Y.; Yi, S.; Xu, J.; Long, Y.; Ma, D.; Zhang, Z. The impacts of climate change and human activities on alpine vegetation and permafrost in the Qinghai-Tibet Engineering Corridor. Ecol. Indic.; 2018; 93, pp. 24-35. [DOI: https://dx.doi.org/10.1016/j.ecolind.2018.04.067]

22. Jiang, Y.; Guo, J.; Peng, Q.; Guan, Y.; Zhang, Y.; Zhang, R. The effects of climate factors and human activities on net primary productivity in Xinjiang. Int. J. Biometeorol.; 2020; 64, pp. 765-777. [DOI: https://dx.doi.org/10.1007/s00484-020-01866-4]

23. Woodward, F.I.; Lomas, M.R. Vegetation dynamics--simulating responses to climatic change. Biol. Rev. Camb. Philos. Soc.; 2004; 79, pp. 643-670. [DOI: https://dx.doi.org/10.1017/S1464793103006419] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/15366766]

24. Ardo, J. Comparison between remote sensing and a dynamic vegetation model for estimating terrestrial primary production of Africa. Carbon Balance Manag.; 2015; 10, 8. [DOI: https://dx.doi.org/10.1186/s13021-015-0018-5] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/25960765]

25. Romilly, P. Time series modelling of global mean temperature for managerial decision-making. J. Environ. Manag.; 2005; 76, pp. 61-70. [DOI: https://dx.doi.org/10.1016/j.jenvman.2005.01.008]

26. Yürekli, K.; Simsek, H.; Cemek, B.; Karaman, S. Simulating climatic variables by using stochastic approach. Build. Environ.; 2007; 42, pp. 3493-3499. [DOI: https://dx.doi.org/10.1016/j.buildenv.2006.10.046]

27. Li, X.; Zhang, C.; Zhang, B.; Liu, K. A comparative time series analysis and modeling of aerosols in the contiguous United States and China. Sci. Total Environ.; 2019; 690, pp. 799-811. [DOI: https://dx.doi.org/10.1016/j.scitotenv.2019.07.072]

28. Li, X.; Liu, K.; Tian, J. Variability, predictability, and uncertainty in global aerosols inferred from gap-filled satellite observations and an econometric modeling approach. Remote Sens. Environ.; 2021; 261, 112501. [DOI: https://dx.doi.org/10.1016/j.rse.2021.112501]

29. Yuan, C.; Liu, S.; Fang, Z. Comparison of China’s primary energy consumption forecasting by using ARIMA (the autoregressive integrated moving average) model and GM(1,1) model. Energy; 2016; 100, pp. 384-390. [DOI: https://dx.doi.org/10.1016/j.energy.2016.02.001]

30. Zhang, L.; Lin, J.; Qiu, R.; Hu, X.; Zhang, H.; Chen, Q.; Tan, H.; Lin, D.; Wang, J. Trend analysis and forecast of PM2.5 in Fuzhou, China using the ARIMA model. Ecol. Indic.; 2018; 95, pp. 702-710. [DOI: https://dx.doi.org/10.1016/j.ecolind.2018.08.032]

31. Valipour, M.; Banihabib, M.E.; Behbahani, S.M.R. Comparison of the ARMA, ARIMA, and the autoregressive artificial neural network models in forecasting the monthly inflow of Dez dam reservoir. J. Hydrol.; 2013; 476, pp. 433-441. [DOI: https://dx.doi.org/10.1016/j.jhydrol.2012.11.017]

32. Anav, A.; Friedlingstein, P.; Beer, C.; Ciais, P.; Harper, A.; Jones, C.; Murray-Tortarolo, G.; Papale, D.; Parazoo, N.C.; Peylin, P. et al. Spatiotemporal patterns of terrestrial gross primary production: A review. Rev. Geophys.; 2015; 53, pp. 785-818. [DOI: https://dx.doi.org/10.1002/2015RG000483]

33. Wu, M.; Liu, Y.; Xu, Z.; Yan, G.; Ma, M.; Zhou, S.; Qian, Y. Spatio-temporal dynamics of China’s ecological civilization progress after implementing national conservation strategy. J. Clean. Prod.; 2021; 285, 124886. [DOI: https://dx.doi.org/10.1016/j.jclepro.2020.124886]

34. Jiang, L.; Liu, Y.; Wu, S.; Yang, C. Analyzing ecological environment change and associated driving factors in China based on NDVI time series data. Ecol. Indic.; 2021; 129, 107933. [DOI: https://dx.doi.org/10.1016/j.ecolind.2021.107933]

35. He, Y.; Wang, M. China’s geographical regionalization in Chinese secondary school curriculum (1902–2012). J. Geogr. Sci.; 2013; 23, pp. 370-383. [DOI: https://dx.doi.org/10.1007/s11442-013-1016-8]

36. Zhang, Z.; Chen, X.; Xu, C.-Y.; Yuan, L.; Yong, B.; Yan, S. Evaluating the non-stationary relationship between precipitation and streamflow in nine major basins of China during the past 50 years. J. Hydrol.; 2011; 409, pp. 81-93. [DOI: https://dx.doi.org/10.1016/j.jhydrol.2011.07.041]

37. Madani, N.; Kimball, J.S.; Running, S.W. Improving Global Gross Primary Productivity Estimates by Computing Optimum Light Use Efficiencies Using Flux Tower Data. J. Geophys. Res. Biogeosci.; 2017; 122, pp. 2939-2951. [DOI: https://dx.doi.org/10.1002/2017JG004142]

38. Madani, N.; Parazoo, N.C.; Kimball, J.S.; Ballantyne, A.P.; Reichle, R.H.; Maneta, M.; Saatchi, S.; Palmer, P.I.; Liu, Z.; Tagesson, T. Recent Amplified Global Gross Primary Productivity due to Temperature Increase Is Offset by Reduced Productivity due to Water Constraints. AGU Adv.; 2020; 1, e2020AV000180. [DOI: https://dx.doi.org/10.1029/2020AV000180]

39. Wang, S.; Zhang, Y.; Ju, W.; Qiu, B.; Zhang, Z. Tracking the seasonal and inter-annual variations of global gross primary production during last four decades using satellite near-infrared reflectance data. Sci. Total Environ.; 2021; 755, 142569. [DOI: https://dx.doi.org/10.1016/j.scitotenv.2020.142569]

40. Liang, S.; Zhao, X.; Liu, S.; Yuan, W.; Cheng, X.; Xiao, Z.; Zhang, X.; Liu, Q.; Cheng, J.; Tang, H. et al. A long-term Global LAnd Surface Satellite (GLASS) data-set for environmental studies. Int. J. Digit. Earth; 2013; 6, pp. 5-33. [DOI: https://dx.doi.org/10.1080/17538947.2013.805262]

41. Tramontana, G.; Jung, M.; Schwalm, C.R.; Ichii, K.; Camps-Valls, G.; Ráduly, B.; Reichstein, M.; Arain, M.A.; Cescatti, A.; Kiely, G. et al. Predicting carbon dioxide and energy fluxes across global FLUXNET sites with regression algorithms. Biogeosciences; 2016; 13, pp. 4291-4313. [DOI: https://dx.doi.org/10.5194/bg-13-4291-2016]

42. Meroni, M.; Rossini, M.; Guanter, L.; Alonso, L.; Rascher, U.; Colombo, R.; Moreno, J. Remote sensing of solar-induced chlorophyll fluorescence: Review of methods and applications. Remote Sens. Environ.; 2009; 113, pp. 2037-2051. [DOI: https://dx.doi.org/10.1016/j.rse.2009.05.003]

43. Li, X.; Xiao, J. A Global, 0.05-Degree Product of Solar-Induced Chlorophyll Fluorescence Derived from OCO-2, MODIS, and Reanalysis Data. Remote Sens.; 2019; 11, 517. [DOI: https://dx.doi.org/10.3390/rs11050517]

44. Abatzoglou, J.T.; Dobrowski, S.Z.; Parks, S.A.; Hegewisch, K.C. TerraClimate, a high-resolution global dataset of monthly climate and climatic water balance from 1958–2015. Sci. Data; 2018; 5, 170191. [DOI: https://dx.doi.org/10.1038/sdata.2017.191]

45. Claverie, M.; Matthews, J.; Vermote, E.; Justice, C. A 30+ Year AVHRR LAI and FAPAR Climate Data Record: Algorithm Description and Validation. Remote Sens.; 2016; 8, 263. [DOI: https://dx.doi.org/10.3390/rs8030263]

46. Qu, Y.; Zhang, C.; Wang, D.; Tian, P.; Bai, W.; Zhang, X.; Zhang, P.; Dai, H.; Wu, Q. Comparison of atmospheric CO₂ observed by GOSAT and two ground stations in China. Int. J. Remote Sens.; 2013; 34, pp. 3938-3946. [DOI: https://dx.doi.org/10.1080/01431161.2013.768362]

47. Gelaro, R.; McCarty, W.; Suárez, M.J.; Todling, R.; Molod, A.; Takacs, L.; Randles, C.A.; Darmenov, A.; Bosilovich, M.G.; Reichle, R. et al. The Modern-Era Retrospective Analysis for Research and Applications, Version 2 (MERRA-2). J. Clim.; 2017; 30, pp. 5419-5454. [DOI: https://dx.doi.org/10.1175/JCLI-D-16-0758.1]

48. Wu, W.Y.; Lo, M.H.; Wada, Y.; Famiglietti, J.S.; Reager, J.T.; Yeh, P.J.; Ducharne, A.; Yang, Z.L. Divergent effects of climate change on future groundwater availability in key mid-latitude aquifers. Nat. Commun.; 2020; 11, 3710. [DOI: https://dx.doi.org/10.1038/s41467-020-17581-y]

49. Hamed, K.H. Exact distribution of the Mann–Kendall trend test statistic for persistent data. J. Hydrol.; 2009; 365, pp. 86-94. [DOI: https://dx.doi.org/10.1016/j.jhydrol.2008.11.024]

50. Gocic, M.; Trajkovic, S. Analysis of changes in meteorological variables using Mann-Kendall and Sen’s slope estimator statistical tests in Serbia. Glob. Planet. Chang.; 2013; 100, pp. 172-182. [DOI: https://dx.doi.org/10.1016/j.gloplacha.2012.10.014]

51. Zscheischler, J.; Michalak, A.M.; Schwalm, C.; Mahecha, M.D.; Huntzinger, D.N.; Reichstein, M.; Berthier, G.; Ciais, P.; Cook, R.B.; El-Masri, B. et al. Impact of large-scale climate extremes on biospheric carbon fluxes: An intercomparison based on MsTMIP data. Glob. Biogeochem. Cycles; 2014; 28, pp. 585-600. [DOI: https://dx.doi.org/10.1002/2014GB004826]

52. Zhang, X.; Liu, K.; Wang, S.; Wu, T.; Li, X.; Wang, J.; Wang, D.; Zhu, H.; Tan, C.; Ji, Y. Spatiotemporal evolution of ecological vulnerability in the Yellow River Basin under ecological restoration initiatives. Ecol. Indic.; 2022; 135, 108586. [DOI: https://dx.doi.org/10.1016/j.ecolind.2022.108586]

53. Zhang, X.; Liu, K.; Li, X.; Wang, S.; Wang, J. Vulnerability assessment and its driving forces in terms of NDVI and GPP over the Loess Plateau, China. Phys. Chem. Earth Parts A/B/C; 2022; 125, 103106. [DOI: https://dx.doi.org/10.1016/j.pce.2022.103106]

54. Wang, X.; Lu, C.; Fang, J.; Shen, Y. Implications for development of grain-for-green policy based on cropland suitability evaluation in desertification-affected north China. Land Use Policy; 2007; 24, pp. 417-424. [DOI: https://dx.doi.org/10.1016/j.landusepol.2006.05.005]

55. Liu, J.; Xu, X.; Shao, Q. Grassland degradation in the “Three-River Headwaters” region, Qinghai Province. J. Geogr. Sci.; 2008; 18, pp. 259-273. [DOI: https://dx.doi.org/10.1007/s11442-008-0259-2]

56. Cui, Y. Preliminary estimation of the realistic optimum temperature for vegetation growth in China. Environ. Manag.; 2013; 52, pp. 151-162. [DOI: https://dx.doi.org/10.1007/s00267-013-0065-1]

57. Guha, A.; Han, J.; Cummings, C.; McLennan, D.A.; Warren, J.M. Differential ecophysiological responses and resilience to heat wave events in four co-occurring temperate tree species. Environ. Res. Lett.; 2018; 13, 065008. [DOI: https://dx.doi.org/10.1088/1748-9326/aabcd8]

58. Shi, H.; Chen, J. Characteristics of climate change and its relationship with land use/cover change in Yunnan Province, China. Int. J. Climatol.; 2018; 38, pp. 2520-2537. [DOI: https://dx.doi.org/10.1002/joc.5404]

59. Li, Z.; Chen, Y.; Wang, Y.; Fang, G. Dynamic changes in terrestrial net primary production and their effects on evapotranspiration. Hydrol. Earth Syst. Sci.; 2016; 20, pp. 2169-2178. [DOI: https://dx.doi.org/10.5194/hess-20-2169-2016]

60. Liu, Y.; Xiao, J.; Ju, W.; Zhou, Y.; Wang, S.; Wu, X. Water use efficiency of China’s terrestrial ecosystems and responses to drought. Sci. Rep.; 2015; 5, 13799. [DOI: https://dx.doi.org/10.1038/srep13799]

61. Peng, C.; Zhou, X.; Zhao, S.; Wang, X.; Zhu, B.; Piao, S.; Fang, J. Quantifying the response of forest carbon balance to future climate change in Northeastern China: Model validation and prediction. Glob. Planet. Chang.; 2009; 66, pp. 179-194. [DOI: https://dx.doi.org/10.1016/j.gloplacha.2008.12.001]

62. Sirvydas, A.; Kučinskas, V.; Kerpauskas, P.; Nadzeikienė, J.; Kusta, A. Solar Radiation Energy Pulsations in a Plant Leaf. J. Environ. Eng. Landsc. Manag.; 2010; 18, pp. 188-195. [DOI: https://dx.doi.org/10.3846/jeelm.2010.22]

63. Durand, M.; Murchie, E.H.; Lindfors, A.V.; Urban, O.; Aphalo, P.J.; Robson, T.M. Diffuse solar radiation and canopy photosynthesis in a changing environment. Agric. For. Meteorol.; 2021; 311, 108684. [DOI: https://dx.doi.org/10.1016/j.agrformet.2021.108684]

64. Feng, X.; Fu, B.; Piao, S.; Wang, S.; Ciais, P.; Zeng, Z.; Lü, Y.; Zeng, Y.; Li, Y.; Jiang, X. et al. Revegetation in China’s Loess Plateau is approaching sustainable water resource limits. Nat. Clim. Chang.; 2016; 6, pp. 1019-1022. [DOI: https://dx.doi.org/10.1038/nclimate3092]

65. Wang, S.; Li, G.; Fang, C. Urbanization, economic growth, energy consumption, and CO₂ emissions: Empirical evidence from countries with different income levels. Renew. Sustain. Energy Rev.; 2018; 81, pp. 2144-2159. [DOI: https://dx.doi.org/10.1016/j.rser.2017.06.025]

66. Devaraju, N.; Bala, G.; Caldeira, K.; Nemani, R. A model based investigation of the relative importance of CO₂-fertilization, climate warming, nitrogen deposition and land use change on the global terrestrial carbon uptake in the historical period. Clim. Dyn.; 2015; 47, pp. 173-190. [DOI: https://dx.doi.org/10.1007/s00382-015-2830-8]

67. He, Q.; Zhang, M.; Huang, B. Spatio-temporal variation and impact factors analysis of satellite-based aerosol optical depth over China from 2002 to 2015. Atmos. Environ.; 2016; 129, pp. 79-90. [DOI: https://dx.doi.org/10.1016/j.atmosenv.2016.01.002]

68. Levy, R.C.; Mattoo, S.; Munchak, L.A.; Remer, L.A.; Sayer, A.M.; Patadia, F.; Hsu, N.C. The Collection 6 MODIS aerosol products over land and ocean. Atmos. Meas. Tech.; 2013; 6, pp. 2989-3034. [DOI: https://dx.doi.org/10.5194/amt-6-2989-2013]

69. Li, X.; Zhang, C.; Li, W.; Liu, K. Evaluating the Use of DMSP/OLS Nighttime Light Imagery in Predicting PM2.5 Concentrations in the Northeastern United States. Remote Sens.; 2017; 9, 620. [DOI: https://dx.doi.org/10.3390/rs9060620]

70. Li, X.; Liang, H.; Cheng, W. Spatio-Temporal Variation in AOD and Correlation Analysis with PAR and NPP in China from 2001 to 2017. Remote Sens.; 2020; 12, 976. [DOI: https://dx.doi.org/10.3390/rs12060976]

71. Li, X.; Seth, A.; Zhang, C.; Feng, R.; Long, X.; Li, W.; Liu, K. Evaluation of WRF-CMAQ simulated climatological mean and extremes of fine particulate matter of the United States and its correlation with climate extremes. Atmos. Environ.; 2020; 222, 117181. [DOI: https://dx.doi.org/10.1016/j.atmosenv.2019.117181]

72. Li, X.; Zhang, C.; Li, W.; Anyah, R.O.; Tian, J. Exploring the trend, prediction and driving forces of aerosols using satellite and ground data, and implications for climate change mitigation. J. Clean. Prod.; 2019; 223, pp. 238-251. [DOI: https://dx.doi.org/10.1016/j.jclepro.2019.03.121]

73. Soni, K.; Parmar, K.S.; Kapoor, S.; Kumar, N. Statistical variability comparison in MODIS and AERONET derived aerosol optical depth over Indo-Gangetic Plains using time series modeling. Sci. Total Environ.; 2016; 553, pp. 258-265. [DOI: https://dx.doi.org/10.1016/j.scitotenv.2016.02.075] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/26925737]

74. Gampe, D.; Zscheischler, J.; Reichstein, M.; O’Sullivan, M.; Smith, W.K.; Sitch, S.; Buermann, W. Increasing impact of warm droughts on northern ecosystem productivity over recent decades. Nat. Clim. Chang.; 2021; 11, pp. 772-779. [DOI: https://dx.doi.org/10.1038/s41558-021-01112-8]