Changes in the trends of vegetation net primary productivity in China between 1982 and 2015

We estimated the NPP using the Carnegie-Ames-Stanford Approach (CASA) model. The CASA model is a powerful tool for large-scale NPP estimation based on plant physiology and carbon sequestration process, whose input includes meteorological, soil and vegetation data (Potter et al 1993, Field et al 1995). In this model, NPP (gram carbon, gC) is expressed as the product of the amount of photosynthetically active radiation (APAR, megajoule, MJ) and light utilization efficiency (ε, gC MJ⁻¹) (equation (1)):

$$\begin{eqnarray}&&{\rm{NPP}}\left(x,t\right)={\rm{A}}{\rm{P}}{\rm{A}}{\rm{R}}(x,t)\times \varepsilon (x,t),\end{eqnarray} \tag{ 1 }$$

where NPP(x, t) is net primary production for a specific location (x) and period (t); APAR(x, t) and ε(x, t) can be calculated by equations (2) and (3) respectively.

$$\begin{eqnarray}&&{\rm{A}}{\rm{P}}{\rm{A}}{\rm{R}}\left(x,t\right)={\rm{F}}{\rm{P}}{\rm{A}}{\rm{R}}(x,t)\times {\rm{S}}{\rm{O}}{\rm{L}}(x,t)\times 0.5\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&\varepsilon \left(x,t\right)={T}_{1}(x,t)\times {T}_{2}(x,t)\times W(x,t)\times {\varepsilon }_{{\max }},\end{eqnarray} \tag{ 3 }$$

where FPAR(x, t) is the fraction of photosynthetically active radiation, which can be derived from vegetation type and normalized differential vegetation index (NDVI); SOL(x, t) is the total solar radiation; T₁(x, t) and T₂(x, t) are effects of temperature stress, which are related to extreme temperatures and deviation from optimum temperature respectively; W(x, t) is calculated by the soil moisture submodel and accounts for the effect of water stress; ε_max is the maximum value of light utilization efficiency and is assigned as 0.405 gC MJ⁻¹ in this study (Field et al 1995).

We selected China except for the non-vegetated regions marked by 1:1000 000 vegetation map (ECVMC 2007) as the study area, which covers 91.3% of China's landmass. To meet the requirements of the CASA model and clearly present the spatial differentiation of NPP variations, the study area was classified into ten vegetation types and four climatic regions (He et al 2019), namely temperate continental, temperate monsoonal, subtropical–tropical monsoonal and the Tibetan Plateau (figure S1 is available online at stacks.iop.org/ERL/14/124009/mmedia). Moreover, soil texture data for calculating the water stress coefficient were provided by the Cold and Arid Region Sciences Data Center (http://westdc.westgis.ac.cn/).

We used meteorological data from 2480 conventional and 130 radiation monitoring stations (http://data.cma.cn) to create gridded monthly temperature, precipitation, and solar radiation products. Spatial interpolation was implemented by Anusplin 4.4, which is based on improved thin plate spline method and has been proven to have better interpolation results than inverse distance weighted and kriging methods in China (Qian et al 2010, Hutchinson and Xu 2013). We tried to change the order of the spline and introduce different covariables to obtain the optimal interpolation models. Parameters and deviations derived from general cross validation of all tested models are shown in table S1. Accordingly, we selected a quadratic model using elevation as a covariate for temperature interpolation and a quadratic but non-covariate model for solar radiation. Although designating the distance to the coastline as a covariate could help to obtain a lower error, it also led to an unusual pattern in Northwest China (figure S2). Therefore, we selected quadratic but non-covariate model for precipitation interpolation in practice. Elevation data used for spatial interpolation were provided by CGIAR-CSI (http://srtm.csi.cgiar.org/).

GIMMS NDVI data with high temporal (15 d) and spatial (8 km or 0.083°) resolution have been widely used for monitoring terrestrial vegetation activity (Beck and Goetz 2011, Li et al 2014). In this study, we applied the maximum value composition method to denoise the dataset and obtained monthly data between 1982 and 2015 (Holben 1986).

Using the above data, mean annual NPP values for each pixel derived from CASA model were highly consistent with those of the MOD17A3 dataset (R² = 0.70, p < 0.01, figure S3(a)). In addition, annual averages of the two datasets were quite synchronous (r = 0.684, p < 0.01) (figure S3(b)). Therefore, the estimation in this study was relatively reliable and able to accurately depict the spatiotemporal changes of NPP in China.

2.3.1. Piecewise regression model

Under the influence of climate change and human activity, ordinary linear regression model (equation (4)) is not adequate to describe the change in long time series vegetation activity (Peng et al 2011, Piao et al 2011). Therefore, we introduced the piecewise regression model with a piecewise point to identify the phasic change features of NPP. Unlike the application scenario for which the model was first proposed (Toms and Lesperance 2003), changes in NPP may not continue (Liang et al 2015), so we deformed the model as shown in equation (5)

$$\begin{eqnarray}&&y={\beta }_{0}+{\beta }_{1}\times t+\varepsilon ,\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&y=\left\{\begin{array}{l}{\beta }_{2}+{\beta }_{3}\times t+\varepsilon \,t\leqslant \alpha \\ {\beta }_{4}+{\beta }_{5}\times t+\varepsilon \,t\gt \alpha \end{array}\right.,\end{eqnarray} \tag{ 5 }$$

where y is annual or seasonal NPP; t represents for the year; ε is an error term; α is the piecewise point or the SHY of the NPP trend in our case; β₁ in equation (4) is the average annual increment of NPP during the entire study period; β₃ and β₅ in equation (5) are NPP trends before and after the SHY, respectively; β₀, β₂ and β₄ are constant terms.

The Akaike information criterion (AIC), which comprehensively reflects fitting error and the number of parameters, was employed to compare the performances of the above two models (equation (6))

$$\begin{eqnarray}&&\begin{array}{l}{\rm{AIC}}=n\times \,\mathrm{log}({\rm{R}}{\rm{S}}{\rm{S}}/n)+2\times k+2\times k\\ \,\,\,\,\times \,(k+1)/(n-k-1),\end{array}\end{eqnarray} \tag{ 6 }$$

where RSS is the residual sum of squares; n is the sequence length; k is the number of parameters, which equal to 2 and 5 for equations (4) and (5) respectively.

Burnham and Anderson (2002) noted that the difference between AICs ($\delta {\rm{AIC}}={{\rm{AIC}}}_{1}-{{\rm{AIC}}}_{2}$) and the threshold of −2 can be used to choose the better of two models. Specifically, in this study, we used piecewise regression for model 1 and ordinary linear regression for model 2. Accordingly, $\delta {\rm{AIC}}\leqslant -2$ implied that the piecewise regression model was better and there was a significant change in NPP trend. While $\delta {\rm{AIC}}\gt -2$ suggested that there was no significant difference in the NPP trend before and after the SHY.

2.3.2. Multiple regression model

Taking temperature, precipitation and solar radiation as inputs, multiple linear regression model can commendably simulate changes in NPP in China (Wang et al 2017). Thus, equations (7) and (8) were applied to evaluate the independent impacts of these variables on NPP trends (He et al 2019)

$$\begin{eqnarray}&&{\rm{NPP}}={\beta }_{C}+{\beta }_{T}\times {X}_{T}+{\beta }_{P}\times {X}_{P}+{\beta }_{S}\times {X}_{S}+\varepsilon ,\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&{C}_{i}=d({\beta }_{i}\times {X}_{i})/dt\,i=T,P,S,\end{eqnarray} \tag{ 8 }$$

where X_T, X_P and X_S are sequences of temperature, precipitation and solar radiation respectively; β_T, β_P and β_S are regression coefficients corresponding to three climatic factors; β_C is a constant term; ε is an error term; t represents for the year; C_i is the derivative of the product term with respect to time and denotes the trend of NPP which is directly affected by a certain variable. It is noteworthy that C_T (C_P, C_S) is usually expressed as the contribution of temperature (precipitation, solar radiation) to the NPP trend in this study.

The CASA model and data analyses were all implemented with R (R core team 2019).

The mean annual NPP from 1982 to 2015 was 2530 TgC in China. During the entire study period, the annual NPP significantly increased (R² = 0.68, p < 0.01), with a mean annual increment of 12.63 TgC, or an increase rate of 0.5% (figure 1(a)). Because non-evergreen vegetation, whose NPP was close to zero in winter (December–February), occupies most of China's landmass, we concentrated only on the dynamics of vegetation NPP during the growing season (figures 1(b)–(d)). The mean annual NPP in spring, summer and autumn were 550 TgC (21.7% of annual value), 1279 TgC (50.4%) and 543 TgC (21.4%) respectively. Over the 34 years, NPP in all three seasons increased significantly (p < 0.01). Spring NPP had the largest growth rate (0.83%), followed by that of autumn (0.49%) and summer (0.30%).

Zoom In Zoom Out Reset image size

The results of piecewise regression indicated that there was a significant change in the trend of annual NPP ($\delta {\rm{AIC}}\leqslant -2$). The mean annual increment decreased from 24.74 to 15.06 TgC yr⁻¹ in 1998 (figure 1(a)). At the seasonal scale, significant SHYs were also found in summer and autumn. Before and after 1997, summer NPP increased rapidly by 11.98 TgC yr⁻¹ and then increased slowly by 5.43 TgC yr⁻¹. A sudden change occurred in 1992, when the previous decreasing trend (−1.68 TgC yr⁻¹) of autumn NPP reversed to an increasing trend (1.60 TgC yr⁻¹). Moreover, the shift in spring NPP trend appeared in 1998, but there was a significant difference in trend values before and after the SHY.

From 1982 to 2015, significant SHYs of annual NPP were identified in 63.3% of the country. SHYs in the Tibetan Plateau appeared early, distributed between 1992 and 1995, while those in the temperate continental region occurred relatively late, mostly after 2000 (figure 2(a)). Furthermore, the changes in the magnitude of NPP trends also varied from different regions. NPP growth was accelerated in over 60% of the total area in both the temperate continental region and the subtropical–tropical region, whereas the proportion was 48.8% in the temperature monsoonal region and 51.0% in the Tibetan Plateau (figure 2(b)).

Zoom In Zoom Out Reset image size

As the largest component of annual NPP, change features of summer NPP were similar to those of the annual value. More precisely, the acceleration of NPP growth was mainly concentrated in the boundary between the temperate continental region and the temperate monsoonal region as well as some parts of southern China. The regions with reduced growth rates after SHYs were centered in northeast and southwest China (figures 2(a), (b), 3(c), (d)). Although SHYs of summer and autumn NPP were different from those of the annual NPP, the spatial distribution of the regions where NPP growth was accelerated was highly consistent with the annual situation. Moreover, it was also noteworthy that the trends of autumn NPP were increased after SHYs in some parts of northeast and southwest China (figures 3(a), (b), (e), (f)). More detailed statistical results of phasic change features in four regions can be found in table S2.

Zoom In Zoom Out Reset image size

Our results confirmed the negative effect of the warming hiatus on NPP growth (He et al 2019). SHYs of annual, spring and summer NPPs all appeared in 1997 or 1998, which is similar to the time when global average temperature reached the peak (Easterling et al 2000, Karl et al 2000). Furthermore, the growth rate of annual temperature in China declined from 0.056 °C yr⁻¹ before 1998 to 0.009 °C yr⁻¹, and the NPP trend directly attributed to temperature changed from 9.429 to 1.004 TgC yr⁻¹ in the meantime, which accounted for the largest part of the trend change (tables 1, 2).

In contrast to temperature, China's annual, spring and autumn precipitation increased to varying degrees after SHYs, which created conditions for the acceleration of NPP growth. In autumn, especially, the precipitation trend increased from −3.131 mm yr⁻¹ before 1992 to 0.797 mm yr⁻¹, and its contribution increased synchronously from −1.928 to 0.635 TgC yr⁻¹, which was up to 78.2% of NPP trend change. In contrast, precipitation, whose contribution was the largest among three climatic factors, was proven to be the main cause of decelerated summer NPP growth (tables 1, 2).

Some studies have shown that the increase in solar radiation can aggravate the water stress on vegetation growth in arid areas (Han et al 2012). However, once the influence of temperature was removed statistically, vegetation NPP typically had positive correlation with solar radiation (Wang et al 2017). In our case, annual radiation at the national scale increased by 27.252 MJ yr⁻¹ before 1998 and by 6.322 MJ yr⁻¹ afterwards, which made it another important factor, in addition to temperature, slowing down the annual NPP growth. Moreover, radiation also played a prominent role in shifting the trend of autumn NPP, with its contribution accounting for 28.9% of the value, second only to precipitation (tables 1, 2).

Summarizing pixel-level results of multiple regression by region, we found marked spatiotemporal differences in the response of NPP trends to climatic factors (figures 4, 5). The warming hiatus was the main cause of the stalled annual NPP trend in all regions except for the Tibetan Plateau, where the proportion and amplitude of the reduction in temperature contribution were both at high levels (figure 4(a)). The trends of annual NPP were more sensitive to precipitation variation in the temperate continental region and the temperate monsoonal region, with the proportions of positive changes in both precipitation trend and its contribution reaching to 34.2% and 67.5%, respectively. In contrast, such response was relatively weak in the subtropical–tropical monsoonal region (figure 4(b)). In addition, our results indicated that solar radiation was the determinant of annual NPP change in the subtropical–tropical monsoonal region and the Tibetan Plateau because the proportions of changes with the same direction in trend and contribution were the highest among all three factors (figure 4(c)).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Due to the differences in climatic conditions, NPP trend also had some seasonal characteristics in response to the changes in various factors (figure 5). Warming can advance the onset and delay the end of the growing season in cold areas, increasing NPP (Myneni et al 1997, Richardson et al 2010). For this reason, in the Tibetan Plateau, the region where vegetation growth is chronically restricted by temperature, the proportion of the decreased NPP trend caused by the warming hiatus in spring and summer was highest among all regions (figures 5(a), (d)). The high values of this proportion in the temperate monsoonal region during spring and summer as well as in the subtropical–tropical monsoonal region during summer also implied that the warming hiatus indirectly influenced the change in the NPP trend by regulating the length of the growing season (figures 5(a), (d), (g)). However, warming hiatus was found to facilitate summer NPP growth in two monsoonal regions, for the regions with negative changes in temperature trend and positive changes in contribution accounted for 35.8% and 37.5% of the total areas respectively (figure 5(d)). This suggested that the warming hiatus could partially alleviate the inhibition of vegetation growth caused by summer drought (Zhao and Running 2010). It is noteworthy that in the temperate continental region where vegetation growth is theoretically susceptible to water, the proportion of negative changes in both summer temperature trend and contribution remained considerable (56.9%) (figure 5(d)). The reasons for this are two-fold. First, the elevation of this area is high, and vegetation growth is still mainly limited by temperature. Second, the rising temperature may increase the amount of water available, as glacier melting is a key source of water other than precipitation in this region (Barnett et al 2005, Zhang et al 2016).

The regional scale analysis results showed that the increase in autumn precipitation in recent years was mainly caused by the strengthening of the East Asian monsoon (figure 5(h)), which was consistent with numerous predictions based on climate models (Bao 2012, IPCC 2013). Increased monsoonal precipitation had become the most important factor triggering the acceleration of autumn NPP growth in Eastern China. The proportion of positive changes in both trend and contribution was up to 68.8% in the subtropical–tropical monsoonal region and was 62.5% in the temperate monsoonal region (figure 5(h)). Precipitation is the key factor limiting photosynthesis in the growing season (Poulter et al 2014, Ahlström et al 2015). Moderate increases in precipitation can promote vegetation growth. While excess rainfall, on the one hand, may cause water-logging, resulting in a severe oxygen deficiency in the root zone and the decline in NPP, and on the other hand, it may indirectly affect carbon sequestration of plants by weakening solar radiation (Perata et al 2011, Liu et al 2018). Therefore, in the subtropical–tropical monsoonal region with abundant rainfall, NPP trends in spring and summer were not sensitive to the change in precipitation, and the proportions of opposite changes in trend and contribution were also larger than that in autumn (figures 5(b), (e), (h)). However, although the relationships between precipitation trend and contribution in spring and summer were consistent in this region, their causes were different (figures 5(b), (e)). Under the influence of Indian Ocean monsoon, spring precipitation accounts for a small part of the annual total in southwest China (Yu et al 2018). As a result, vegetation growth there is limited by water in spring and is likely to be inhibited in summer due to excess precipitation. By comparison, the distribution of seasonal precipitation is more balanced in southeast China, where it is controlled by the Pacific monsoon (Yang and Lau 2004). Thus, vegetation growth in this region may experience a water deficit in summer due to extensive evaporation. Multiple regression results based on standardized data fully confirmed the above viewpoint. The responses of NPP to increasing precipitation were exactly opposite in spring and summer (figure S4).

The impacts of solar radiation on vegetation growth were two-sided. A strong positive correlation between the changes in solar radiation and contribution was found in the subtropical–tropical region. Especially in spring and summer, the contributions of solar radiation to the NPP trend were much higher than those of other climatic factors (figures 5(c), (f)). In contrast, the increase in solar radiation inhibited on NPP growth in the western part of China, and the proportions of opposite changes in trend and contribution were higher than those in the eastern part (figures 5(c), (f), (i)). Such spatial differences were probably related to the level of radiation (figure S5). Under low exposure, the increase in solar radiation has more of a promoting effect on vegetation growth than inhibitory effect. However, when the level is high, the opposite is true.

In recent years, human activity, mainly in the form of afforestation and grass planting, is profoundly affecting the carbon sequestration functioning and patterns of terrestrial ecosystems in China (Fang et al 2014, Lu et al 2018). Accordingly, the spatial distribution of several major ecological projects was collected and used to explore the connection between the implementation of these projects and the phasic change in annual NPP (figures S6, S7). The results showed that the start time of projects was partially consistent with SHYs of annual NPP. The regions where the absolute values of time difference were less than three accounted for 38.2% of the project area (figure S7(a)). In addition, the positive proportion and the median value of NPP trend changes were both higher than those in other regions (figure S7(b)). Further calculation showed that the regions with increased NPP trends after the implementation of projects accounted for 63.8% of the project area, but only 8.9% of these relationships were significant (figures S7(c), (d)), which indicated that these ecological projects played limited roles in promoting NPP growth. Moreover, significant stalled trends of NPP in parts of northwest China also confirmed that unsuitable restoration projects might harm local vegetation (figures S7(c), (d)) (Feng et al 2015).

In order to identify climatic or anthropogenic drivers controlling the changes in the NPP trends, it is necessary to exclude the effect of the potential inconsistency of GIMMS NDVI dataset (Tian et al 2015). The T-test was carried out to compare the 5-year average NDVI at the national scale around the year 2000, because the sensor was changed from AVHRR-2 to AVHRR-3 (figure S8) at around the year. The result indicated that there was no significant difference between the two values (t = 1.030, p > 0.05) and the dataset was continuous in 2000. We also applied piecewise regression to the sequence of NDVI, and found that the SHY appeared in 1990, which was consistent with the previous study by Peng et al (2011) and was different from that of NPP (1998) (figures 1, S8). Therefore, the potential inconsistency of GIMMS NDVI is likely playing a limited role in triggering the shifting of NPP, while climate changes and human activities could be the key factors. Additionally, we blended different NPP products to verify the stability of SHY. Using the combined sequence which consists of CASA NPP (1982–1999) and corrected MODIS NPP (2000–2014), we found that the SHY of annual NPP still appeared in 1998 (figure S9). This result also suggested that spectral data itself didn't significantly influence the SHY of NPP.