Divergent historical GPP trends among state-of-the-art multi-model simulations and satellite-based products

Understanding historical changes in gross primary productivity (GPP) is essential for better predicting the future global carbon cycle. However, the historical trends of terrestrial GPP, due to the CO₂ fertilization effect, climate, and land-use change, remain largely uncertain. Using long-term satellite-based near-infrared radiance of vegetation (NIRv), a proxy for GPP, and multiple GPP datasets derived from satellite-based products, dynamic global vegetation model (DGVM) simulations, and an upscaled product from eddy covariance (EC) measurements, here we comprehensively investigated their trends and analyzed the causes for any discrepancies during 1982–2015. Although spatial patterns of climatological annual GPP from all products and NIRv are highly correlated (r>0.84), the spatial correlation coefficients of trends between DGVM GPP and NIRv significantly decreased (with the ensemble mean of r=0.49) and even the spatial correlation coefficients of trends between other GPP products and NIRv became negative. By separating the global land into the tropics plus extratropical Southern Hemisphere (Trop+SH) and extratropical Northern Hemisphere (NH), we found that, during 1982–2015, simulated GPP from most of the models showed a stronger increasing trend over Trop+SH than NH. In contrast, the satellite-based GPP products indicated a substantial increase over NH. Mechanistically, model sensitivity experiments indicated that the increase of annual global total GPP was dominated by the CO₂ fertilization effect (83.9 % contribution), however, with the largest uncertainty in magnitude in individual simulations among the three drivers of CO₂ fertilization, climate, and land-use change. Interestingly, the spatial distribution of inter-model spreads of GPP trends resulted mainly from climate and land-use change rather than CO₂ fertilization effect. After 2000, trends from satellite-based GPP products were different from the full time series, suggesting weakened rising trends over NH and even significantly decreasing trends over Trop+SH, while the trends from DGVMs and NIRv kept increasing. The inconsistencies of GPP trends are very likely caused by the contrasting performance between satellite-derived and DGVM simulated vegetation structure parameter (leaf area index, LAI). Therefore, the uncertainty in satellite-based GPP products induced by highly uncertain LAI data in the tropics undermines their roles in assessing the performance of DGVM simulations and understanding the changes of global carbon sinks. The higher consistency between DGVM GPP and NIRv suggests that the trends from a DGVM ensemble might even have better performance than satellite-based GPP products.

1 Introduction

The gross primary productivity (GPP) is the largest carbon flux in the terrestrial carbon cycle. Quantifying terrestrial GPP and understanding its variations are vital in the global and regional carbon cycle (Ryu et al., 2019). To date, there are multiple global GPP products, mainly including the upscaled products from the eddy covariance (EC) flux data by machine learning techniques (Beer et al., 2010; Jung et al., 2020), satellite-based estimates by light-use efficiency (LUE) model (Running et al., 2004; Yuan et al., 2010; Joiner et al., 2018; Zheng et al., 2020), and simulations by the state-of-the-art dynamic global vegetation models (DGVMs) (Huntzinger et al., 2013; Sitch et al., 2015).

The machine learning FLUXCOM GPP products based on the global FLUXNET network, remote sensing, and meteorological input (Jung et al., 2020; Pastorello et al., 2020) are widely used in terrestrial carbon cycle studies. Taking FLUXCOM GPP as a benchmark, researchers have explored the interannual variation, seasonal cycle, and climatology pattern of global and regional GPP (Chen et al., 2017; Jia et al., 2020; Zhang and Ye, 2021). However, due to the lack of the CO₂ fertilization effect, the performance of this product on the long-term GPP trend is not realistic (Jung et al., 2020). Based on the LUE principle and derived from the Advanced Very High-Resolution Radiometer (AVHRR) and the Moderate Resolution Imaging Spectroradiometer (MODIS) datasets, the satellite-based GPP estimates include MOD17, GLASS, GIMMS, FluxSat, WECANN, and the revised EC_LUE GPP product (Running et al., 2004; Yuan et al., 2007; Smith et al., 2016; Alemohammad et al., 2017; Joiner et al., 2018; Zheng et al., 2020). These GPP products capture the seasonal variation, spatial distribution, and interannual variation to a large extent (Wang et al., 2014) but do not always account for the CO₂ fertilization effect (O'Sullivan et al., 2020). For DGVM simulations, different forcing datasets, parameterizations, and processes considered can make the surprising differences in model representation of responses of photosynthesis to CO₂ concentration, soil moisture, temperature, and water vapor deficit (Rogers, 2014; Rogers et al., 2017; Ito et al., 2017; M. Wang et al., 2021). These differences caused large inter-model spreads in GPP simulations (Ito et al., 2017). Hence, many efforts have been made to constrain the global GPP magnitude based on the satellite observations like solar-induced chlorophyll fluorescence (SIF) (Macbean et al., 2018; Bacour et al., 2019; Norton et al., 2019; J. Wang et al., 2021a).

The application of satellite-derived GPP proxy datasets provides a breakthrough for estimating global GPP (Running et al., 2004; Badgley et al., 2019; Piao et al., 2020). Many GPP proxy indices, such as normalized difference vegetation index (NDVI), enhanced vegetation index (EVI), and SIF, have been widely used to estimate the global GPP (Frankenberg et al., 2011; Guanter et al., 2014). However, each of them has its shortcomings. For example, NDVI can be saturated in tropical regions, demonstrating its nonlinear relationship with GPP (Badgley et al., 2017, 2019; Camps-Valls et al., 2021). The EVI index improves the NDVI algorithm, but this index has not entirely solved the saturation problem (Huete et al., 2002). Without dealing with the problem of distinguishing whether the signal comes from the plant or other interference factors, satellite retrieval of SIF measures the light emitted by chlorophyll in leaves and can be used as a robust proxy of GPP (Frankenberg et al., 2011; Mohammed et al., 2019). However, the time range of global SIF products is short, with direct observations only available from 2007. Representing the proportion of reflected near-infrared radiation attributable to vegetation, long-term satellite-derived near-infrared radiance of vegetation (NIRv) is a relatively recent GPP proxy (Badgley et al., 2017). Compared to NDVI and EVI, the saturation problem of NIRv and GPP in the tropical region is weakened because the mixed effects of background brightness, leaf area, and the distribution of canopy photosynthetic capacity with depth are largely eliminated. Since NIRv can be directly obtained from observational datasets of the AVHRR sensors, it can be derived from 1982 to the present. Moreover, previous studies have shown that NIRv and SIF are closely related and indicated that NIRv could well represent changes in GPP (Badgley et al., 2017, 2019; Camps-Valls et al., 2021; Wang et al., 2020).

Although there have been a lot of studies focusing on extreme anomalies, the seasonal cycle, interannual variation, and the climatological pattern of global and regional GPP based on the multiple GPP products and proxy indices (Chen et al., 2017; Madani et al., 2020; J. Wang et al., 2021b), few efforts have been devoted to evaluate the long-term GPP trends across different GPP sources and to analyze the causes of uncertainties. This study will comprehensively investigate historical GPP trends during 1982–2015, based on the satellite-derived GPP proxy (NIRv), simulations from process-based models, machine-learning products, satellite-based estimates, and site-level observations. Hereinafter, Sect. 2 describes the datasets and statistical methods used. The comparison of GPP trends among GPP proxy, DGVM simulations, and satellite-based GPP products is in Sect. 3.1. The mechanisms of the trend attributions are explored in Sect. 3.2. The uncertainties in GPP trends are discussed in Sect. 3.3. At last, the main conclusions of the results are summarized in Sect. 4.

2 Datasets and methods

2.1 TRENDYv6 multi-model-simulated GPP

We used the model simulation results conducted under the auspices of the “Trends and drivers of the regional scale sources and sinks of carbon dioxide” (TRENDY) Project (Sitch et al., 2015). We used 10 DGVMs in the TRENDYv6 project for the period of 1982–2015, including CABLE (Haverd et al., 2018), CLASS-CTEM (Melton and Arora, 2016), CLM4.5 (Oleson et al., 2013), DLEM (Tian et al., 2015), ISAM (Jain et al., 2013), OCN (Zaehle et al., 2010), ORCHIDEE-MICT (Guimberteau et al., 2018), ORCHIDEE (Krinner et al., 2005), VEGAS (Zeng et al., 2005), and VISIT (Kato et al., 2013). There is a suite of experimental protocols in the TRENDY project, and we here explored GPP trends and their mechanisms using the GPP outputs from three simulations. In detail, DGVMs were run under the varying CO₂ concentration, and constant climate conditions and land-use change in experiment S1 in the Supplement; the varying CO₂ concentration and climate conditions, with constant land-use change in experiment S2; the varying CO₂ concentration, climate conditions, and land-use change in experiment S3. Hence, experiment S1 represents the impact of the CO₂ fertilization effect. The contributions of climate change and land-use change (hereafter “LUC”) are calculated through the differences between experiments S2 and S1 and experiments S3 and S2, respectively. These modeling details are listed in Table 1.

In Sect. 3.2 and 3.3.1, we calculated the ensemble mean of the seven model simulations, which included all scenarios as the DGVM ensemble GPP and calculated their standard deviation to represent inter-model spread across these models. In other sections which only need results from experiment S3, we used the ensemble mean simulations from 10 models.

Table 1 Information of TRENDYv6 models used in this study. S1 represents the impact of the CO₂ fertilization effect, S2 represents the impact of the CO₂ fertilization effect and climate change, and S3 represents the impact of the CO₂ fertilization effect, climate change, and LUC.

^a Simulation datasets in the corresponding experiments (S1, S2, and S3) as available for models indicated with the notation of “✓”.

2.2 FLUXCOM GPP

The FLUXCOM datasets comprised of 120 global carbon flux products generated by nine machine learning algorithms based on site-level observed GPP measured by EC associated with remote sensing information and meteorology data but did not take the CO₂ fertilization effect into account (Jung et al., 2020). This research used the ensemble mean of GPP datasets forced by the Climatic Research Unit – Japanese Reanalysis (CRU JRA) climate data and generated from three machine learning techniques (random forest, artificial neural network, and multivariate adaptive regression splines) from 1982 to 2015. The original spatial resolution of this dataset is $\mathrm{0.5}{}^{\circ }\times \mathrm{0.5}{}^{\circ }$.

2.3 Satellite-based GPP products

In this study, the GLASS GPP and revised EC-LUE GPP estimates were used as representatives of long-term satellite-based GPP products from 1982 to 2015. GLASS GPP originated from the eddy covariance–light-use efficiency (EC-LUE) model (Yuan et al., 2007), which considered various impact factors (NDVI, photosynthetically active radiation, temperature, CO₂ concentrations, the Bowen ratio of sensible to latent heat flux, water vapor pressure deficit (VPD), direct radiation fluxes, and diffuse radiation fluxes) and nine ecosystem types to accurately estimate the long-term change of GPP (Yuan et al., 2019). The original spatial resolution of this dataset is $\mathrm{0.05}{}^{\circ }\times \mathrm{0.05}{}^{\circ }$.

The revised EC-LUE GPP is a long-term GPP dataset based on the LUE equation. Zheng et al. (2020) generated the revised EC-LUE GPP using the following formula:

$$\begin{array}{}\text{(1)}& \begin{array}{rl}\mathrm{GPP}& =\left({\mathit{\epsilon }}_{\mathrm{msu}}\times \mathrm{APARsu}+{\mathit{\epsilon }}_{\mathrm{msh}}\times {\mathrm{APAR}}_{\mathrm{sh}}\right)\\ & \times \mathrm{Cs}\times \min \left(Ts,\mathrm{Ws}\right),\end{array}\end{array}$$

where ε_msu is the maximum LUE of sunlit leaves and ε_msh is the maximum LUE of shaded leaves; APARsu is the PAR absorbed by sunlit leaves, and APAR_sh is the PAR absorbed by shaded leaves. C_s, T_s, and W_s represent the downward regulation scalars of atmospheric CO₂ concentration, temperature, and VPD on LUE with the range from 0 to 1. Specifically, the direct effect of CO₂ fertilization on GPP is determined by the following equations:

$$\begin{array}{}\text{(2)}& {\displaystyle }{C}_{\mathrm{s}}={\displaystyle \frac{C_{\mathrm{i}}-\mathit{\phi }}{C_{\mathrm{i}}+\mathrm{2}\mathit{\phi }}}\text{(3)}& {\displaystyle }{C}_{\mathrm{i}}=C_{\mathrm{a}}\times \mathit{\chi },\end{array}$$

where C_i represents the CO₂ concentration inside the leaf, C_a is the atmospheric CO₂ concentration, φ means the CO₂ compensation point in the absence of dark respiration (ppm), and χ means the ratio of CO₂ concentration inside the leaf to that in the atmosphere (Farquhar et al., 1980). After adding the effect of CO₂ fertilization, Zheng et al. (2020) suggested that the generated GPP was closer to the site observation data of more than 5 years (R²=0.44) than that of other LUE models (R² ranged from 0.06 to 0.30). The original spatial resolution of this dataset is $\mathrm{0.05}{}^{\circ }\times \mathrm{0.05}{}^{\circ }$.

The spatial pattern and temporal changes of the GLASS GPP and revised EC-LUE GPP are highly consistent (Figs. S1, S2, and S3). Therefore, for simplicity, we averaged them to represent satellite-based GPP products.

2.4 NIRv dataset

The long-term (1982–2015) satellite growing season NIRv dataset used in this study as the GPP proxy (Badgley et al., 2019) originates from the AVHRR sensors. Following the previous study, the NIRv was first calculated as a function of monthly NDVI and near-infrared reflection of the total pixel (NIRT) via the equation of $\mathrm{NIRv}=(\mathrm{NDVI}-\mathrm{0.08})\times \mathrm{NIRT}$ (Badgley et al., 2017). By defining the growing season as monthly average temperature higher than 0 ^∘C, the NIRv for each grid cell during the growing season was then aggregated into annual value. Wang et al. (2020) have shown that AVHRR NIRv could explain about 60 % of the monthly variances in ground observational GPP from the FLUXNET2015 dataset, and NIRv trends could be used as the proxy of GPP trends under the different definitions of the growing season.

2.5 Site-level GPP observations

We also adopted 20 EC sites from the FLUXNET2015 dataset (Pastorello et al., 2020) with an observation period longer than 15 years to evaluate the performance of different global GPP products. These sites, all located over Northern Hemisphere, included five vegetation types: evergreen broadleaf forest (EBF), evergreen needleleaf forest (ENF), deciduous broadleaf forest (DBF), grassland (GRA), and mixed forest (MF) (Table 2). The GPP variable used in this study is GPP_NT_VUT_REF. When evaluating the global gridded GPP datasets with the site observations, the bilinear interpolation method was used to interpolate the gridded data to the specific site locations.

Figure 1 Geographical distributions of linear trends of NIRv and GPP during 1982–2015. (a) AVHRR NIRv, (b) ensemble mean of TRENDY multi-model-simulated GPP, (c) FLUXCOM GPP, (d) mean of satellite-based products from revised EC-LUE and GLASS GPP. Striped areas indicate that the trends are significant with p<0.05 following the non-parametric Mann–Kendall trend test. The trends of NIRv and GPP are unitless and in kgC m⁻² yr⁻², respectively. Additionally, due to a lack of the CO₂ fertilization effect in FLUXCOM GPP (c), we used a smaller scale than in panels (b) and (d).

Table 2 FLUXNET sites used in this study. The vegetation types include evergreen broadleaf forest (EBF), evergreen needleleaf forest (ENF), deciduous broadleaf forest (DBF), grassland (GRA), and mixed forest (MF).

2.6 Leaf area index

GLASS LAI version 05 was used to compare the TRENDY model ensemble LAI (S3) because it is an input parameter for GLASS GPP and revised EC-LUE GPP. This dataset originated from version 4 of the Long-Term Data Record (LTDR) AVHRR surface reflectance product before 2001 with a spatial resolution of $\mathrm{0.05}{}^{\circ }\times \mathrm{0.05}{}^{\circ }$ and MODIS surface reflectance product (MOD09) after 2001 with a spatial resolution of 1 km × 1 km. The spatial-average method was used to aggregate the dataset to $\mathrm{0.05}{}^{\circ }\times \mathrm{0.05}{}^{\circ }$. Biome-specific general regression neural networks were used to fuse these two datasets to generate a long-term LAI product (1982–2018), which improved performance than the original datasets. Its spatial and temporal resolutions are $\mathrm{0.05}{}^{\circ }\times \mathrm{0.05}{}^{\circ }$ and 8 d, respectively (Xiao et al., 2016). The previous study has shown that this product performed well than other long-time LAI estimations based on the evaluation of high-resolution reference maps at VAlidation of Land European Remote sensing Instruments (VALERI) sites (Xiao et al., 2017).

Figure 2 Latitudinal profiles of trends of annual zonal total NIRv and GPP (0.5^∘ latitudinal bins). Results for the DGVM ensemble mean (blue), FLUXCOM GPP (green), satellite-based products (orange), and NIRv (red) during 1982–2015 (a) and 2001–2015 (b), respectively. The GPP values are given on the bottom axis and the NIRv values on the top axis. The shaded areas represent the standard deviation of the individual TRENDY model-simulated GPP trend. The units of the NIRv and GPP trends are unitless and GtC yr⁻², respectively. Panel (c) represents the latitudinal total vegetated land areas (m²).

2.7 Statistical methods used

Due to the difference among spatial resolution of each product, we resampled all GPP datasets into $\mathrm{0.5}{}^{\circ }\times \mathrm{0.5}{}^{\circ }$ through the first-order conservative remapping method:

$$\begin{array}{}\text{(4)}& {\overline{F}}_k={\displaystyle \frac{\mathrm{1}}{A_k}}\int f\mathrm{dA},\end{array}$$

where ${\overline{F}}_k$ is the area-averaged destination quantity, A_k indicates the area of grid k, and f is the quantity in an old grid with an overlapping area with the destination grid. We resampled NIRv and LAI into $\mathrm{0.5}{}^{\circ }\times \mathrm{0.5}{}^{\circ }$ using bilinear interpolation and generated the annual datasets according to the weights of days. To detect the historical changes in GPP in each dataset, we calculated the global and regional total GPP and their linear trends. We also calculated the linear trends of each dataset at the pixel level to generate the spatial patterns of GPP trends. The linear trend was calculated as

$$\begin{array}{}\text{(5)}& y=kx+b+\mathit{\epsilon },\end{array}$$

where k represents the linear trend of the time series, b is the intercept, and ε is the error term.

Finally, non-parametric Mann–Kendall trend tests were used to evaluate the level of significance for each GPP time series because it does not require the data to follow the normal distribution (Hamed, 1998).

Figure 3 Linear trends of global and regional total GPP. Changes of annual total GPP relative to 1982, based on DGVM ensemble mean (blue), FLUXCOM (green), and satellite-based products (orange), compared to AVHRR-NIRv (red), for global (a), Trop+SH (b), and NH (c) measurements. The shaded areas denote the TRENDY inter-model 1σ spread. (d–e) Global and regional GPP trends in individual models and products for the period of 1982–2015 and 2001–2015, respectively. Asterisks indicate that the trend is significant with p<0.05 following the non-parametric Mann–Kendall trend test.

3 Results and discussion

3.1 Different GPP trends in DGVMs and satellite-based products

During 1982–2015, the spatial patterns of climatological annual GPP from different products are highly correlated with satellite-derived NIRv with their spatial correlation coefficients ranging from 0.84 to 0.95 (Table 3). However, the spatial distributions of the various GPP and NIRv trends are quite different. The spatial patterns of trends of NIRv, the DGVM ensemble, FLUXCOM, and satellite-based GPP during 1982–2015 are presented in Fig. 1. NIRv clearly shows increasing trends in most land regions, especially in the northwest parts of Eurasia, and shows decreasing trends over Alaska and Kazakhstan (Fig. 1a). The global distribution of the DGVM ensemble GPP trends is generally consistent with satellite-derived NIRv with their spatial correlation coefficient (r) of 0.49. However, the increasing trends of the DGVM ensemble GPP in the tropical regions (Amazon and equatorial Africa) are higher than those in the boreal zone. Further, for DGVM ensemble GPP, there are about 59.8 % of the global land regions showing significant positive GPP trends and 3.8 % showing significant negative GPP trends. For NIRv, 88.2 % of the global land had positive GPP trends (Table 3).

Although the DGVM ensemble GPP trends are close to those of NIRv, inconsistencies exist in spatial distribution and magnitude of GPP trends among individual model simulations. Firstly, the spatial correlations among individual models and NIRv range from 0.15 to 0.48. Secondly, the GPP simulated by the DLEM shows increasing trends in about three-fourths (77.3 %) of the global land areas, while the GPP simulated by the CLASS-CTEM has increasing trends in only about one-third (33.9 %) of the land area (Table 3). Finally, in magnitude, the trends of VEGAS GPP appear generally weaker than those of other models (Figs. S2, S4a).

The spatial distributions of the trends between NIRv and the remaining non-DGVM products are quite different, ranging from uncorrelated to negatively correlated (Table 3). For FLUXCOM, due to the lack of CO₂ fertilization effect, the GPP trend pattern generally shows no significant trends over 74.1 % of the land areas. And the most striking differences between NIRv and satellite-based GPP products are located in the low latitudes, especially over the Amazon and Indonesia, with the latter indicating significant decreasing trends over these regions (Fig. 1d).

Table 3 Spatial information for NIRv and different GPP products.

Figure 4 Attributions of global total GPP trends for TRENDYv6 simulations from 1982 to 2015: CO₂ fertilization effect (S1), climate (S2 – S1), and land-use change (S3 – S2). “All” gives the values of the reference simulation that includes the effect of all three drivers (S3). Asterisks indicate that the trend is significant with p<0.05 following the non-parametric Mann–Kendall trend test.

Figure 5 The spatial distribution patterns and the inter-model spreads of GPP trends from the DGVM ensemble. (a–b) GPP trends and spreads due to the CO₂ fertilization effect; (c–d) GPP trends and spreads due to climate change; (e–f) GPP trends and spreads owing to LUC; (g–h) GPP trends and spreads from the combined effects of all drivers (S3). Striped areas indicate that the trend is significant with p<0.05 following the non-parametric Mann–Kendall trend test.

As mentioned above, the trends of DGVM ensemble GPP and NIRv show relatively consistent patterns for 1982–2015. However, the DGVM ensemble mean GPP shows slightly stronger trends over the tropics than over the Northern Hemisphere, whereas NIRv has relatively stronger increasing trends over the Northern Hemisphere (Fig. 2a). The tropical region shows the most extensive inter-model spread for the DGVM ensemble, with the strongest trends in CABLE and the weakest trends in VEGAS (Figs. S4a, S5b). The latitudinal distribution of satellite GPP products is quite different from DGVM ensemble and NIRv. Satellite-based GPP shows a significant decrease over the tropical region, rebounding to the most substantial increase located between 15 and 25^∘ S. The GPP increases in the middle and high latitudes of the Northern Hemisphere, but its magnitude is weaker than that of the DGVM ensemble in the most northern regions.

The GPP trends are quite different for the long-term (1982–2015) and recent short-term (2001–2015) periods (Figs. 2b, S4, S6). Compared to the comparable trend magnitudes of DGVM ensemble mean GPP and NIRv over NH during these two periods, the DGVM ensemble mean GPP trends show much stronger increase but NIRv appears a little weakened increase over tropical regions during 2001–2015 (Fig. 2a and b). Additionally, satellite-based GPP products indicate a much stronger GPP decrease over the tropical regions and no noticeable trends in midlatitudes of the Northern Hemisphere after 2000 (Fig. 2b). Also, FLUXCOM GPP trends are near zero in most latitudinal bands during these two periods (Fig. 2a and b) due to the lack of CO₂ fertilization effect (Jung et al., 2020).

The relative changes and linear trends of annual total GPP among the DGVM simulations, FLUXCOM, and satellite-based GPP products vary substantially both globally and regionally (Figs. 3, S3). Based on the analysis over the 34 years, the trend of global GPP was about 0.37±0.08 (DGVM ensemble mean ± 95 % confidence intervals), 0.0 (FLUXCOM), and 0.18 GtC yr⁻² (satellite-based GPP products), respectively (Fig. 3d). Before 2001, the GPP trend of satellite-based products was generally stronger than that of DGVMs (Fig. 3a). Separating the global land into the tropics plus extratropical Southern Hemisphere (Trop+SH: 90^∘ S–23^∘ N) and extratropical Northern Hemisphere (NH: 23–90^∘ N), results show that more than half of the increase of GPP in the DGVM ensemble is from the Trop+SH (57 %). In comparison, the increase of GPP in satellite-based GPP products is mainly attributed to the NH (60 %) rather than Trop+SH (40 %) (Fig. 3d). In individual models, with the exception of ORCHIDEE-MICT and VEGAS, the others indicate that Trop+SH largely contributes to the global GPP trend (ranging from 54.3 % to 65.3 %) (Fig. 3d).

After 2000, there were obvious differences in the trend between DGVM ensemble and satellite-based GPP products, with the satellite-based GPP products showing an obvious turning point at the year of 2000 (Fig. 3a–c), confirming the previous studies (Yuan et al., 2019; Madani et al., 2020). Both GLASS and revised EC-LUE GPP changed from significant increasing trends to significant decreasing trends, resulting mainly from Trop+SH (Fig. 3a–c). Studies based on satellite-based GPP products suggested that this transition was mainly due to the increased atmospheric vapor pressure deficit in the tropical zones (Yuan et al., 2019; Madani et al., 2020). Meanwhile, the increasing GPP trend in the NH was greatly weakened (from 0.10 to 0.01 GtC yr⁻² for GLASS and from 0.11 to 0.05 GtC yr⁻² for revised EC-LUE). In contrast, DGVM ensemble mean GPP and NIRv kept increasing. However, in detail, 4 out of 10 DGVM models (CLASS-CTEM, OCN, ORCHIDEE-MICT, and VEGAS) simulated weakened GPP increasing signals primarily from the Trop+SH (Fig. 3d and e).

In general, we found that, for the two study periods, the trends from DGVM GPP have higher consistencies with the trends of NIRv at the global and regional scales compared with satellite-based GPP. It may suggest that long-term trends of GPP from DGVM ensemble results have better performance than satellite-derived GPP products.

Figure 6 The spatial distributions of LAI trends (m² m⁻² yr⁻²) from the (a) DGVM ensemble mean and (b) GLASS from 2001 to 2015. Striped areas indicate that the trend is significant with p<0.05 following the non-parametric Mann–Kendall trend test.

3.2 Trend attributions in DGVMs

We analyzed the contributions of three drivers (CO₂ fertilization effect, climate change, and LUC) to the GPP trends during 1982–2015 by using the results of the TRENDY sensitivity experiments (Fig. 4). Globally, DGVM ensemble results suggest that the CO₂ fertilization effect is the dominant driver to the increasing GPP (0.29 GtC yr⁻² accounting for 83.9 %), followed by climate change (0.09 GtC yr⁻² accounting for 26.5 %). Additionally, LUC has little effect on the trend of GPP (−0.04 GtC yr⁻² accounting for −10.4 %). For individual model simulations, the contributions from CO₂ fertilization effect, climate change, and LUC range from 65.7 % to 116.3 %, 2 % to 50.4 %, and −18.4 % to 5.5 %, respectively. The model simulations for the period of 2001–2015 show similar results (Fig. S7).

The spatial distributions of GPP trends indicate that the CO₂ fertilization effect consistently increases global GPP, especially in tropical rainforest areas (Fig. 5a). Climate change has inhomogeneous effects on GPP due to different regional changes of climate elements (Figs. 5c, S8a and b) associated with different vegetation sensitivities to each climate element (Wang et al., 2016; Jung et al., 2017). For instance, in the high northern latitudes, global warming dominates the increase of GPP (Figs. 5c, S8a, and S9c); the increase in temperature and decrease in precipitation over the Amazon lead to the GPP decrease; the increase of GPP over equatorial Africa appears to be consistent with the increase of the precipitation based on visual comparison against the trends for the CRU observational dataset (Figs. 5c and S8b). Concentrated in the Trop+SH, LUC basically weakens GPP (Figs. 5e and S9a), mainly due to deforestation (Friedlingstein et al., 2019).

3.3 Uncertainties in GPP trends

3.3.1 DGVM simulations

By analyzing the simulation results driven by each factor, it can be found that though the CO₂ fertilization effect has the most considerable contribution to the global GPP trend, it has the largest inter-model uncertainty (σ=0.11 GtC yr⁻²) among three drivers (Fig. 4). A previous study showed that some models might overestimate the CO₂ fertilization effect on stomatal closure (Anav et al., 2015). The spatial pattern of the standard deviation among each model indicates that the inter-model spreads are mainly located in the Trop+SH (Fig. 5h). The inter-model spread attributed to the CO₂ fertilization effect shows a consistently positive effect on GPP at the global scale (Fig. 5a and b). By contrast, the inter-model spreads driven by climate change and LUC over Trop+SH outweigh the CO₂ fertilization effect (Fig. 5b, d, and f). However, because the inhomogeneous impacts on GPP from climatic elements can offset each other to a large extent (Fig. 5c), it makes the largest inter-model uncertainty of CO₂ fertilization to the global GPP increase rather than the climate effect. Meanwhile, the largest uncertainties in the impact of LUC on GPP trends concentrate over 20^∘–40^∘ S in South America.

We further calculated the spatial correlation coefficients among the GPP trends of each model experiment S3 simulation to quantify their spatial consistencies. The correlation coefficients in pairs among individual models varied from 0.16 to 0.61 (Fig. S10), implying that large uncertainties existed in the distribution of GPP trends among models, which were caused by differences in model structures and parameterizations (Rogers, 2014; Rogers et al., 2017). Furthermore, studies have shown that the global GPP increase can be largely overestimated without nitrogen (N) constraints, especially in the tropical region, where the nitrogen limitation will reduce the photosynthetic capacity of vegetation and weaken its response to the increasing atmospheric CO₂ concentration (He et al., 2017; Terrer et al., 2019). Phosphorus (P) availability also limits the extent to which plants respond to the CO₂ fertilization effect, which is especially relevant in the Amazon rain forest (Fleischer et al., 2019). Therefore, DGVM ensemble GPP may overestimate the increasing GPP trend in the tropical regions since not all of the models used in this study take the effect of N limitation into consideration, and no model includes P limitation in the CO₂ fertilization effect.

Figure 7 Comparisons of annual GPP over different FLUXNET2015 sites (black), FLUXCOM (green), satellite-based product (orange), DGVM ensemble (blue), and NIRv (red). The global GPP datasets were interpolated into the locations of these 20 sites according to the bilinear interpolation method. Observation sites with significant trends are marked with values. Single (^∗) and double (${}^{\ast \ast }$) asterisks indicate that the trend is significant with p<0.1 and p<0.05 following the non-parametric Mann–Kendall trend test. The units of GPP and GPP trend are kgC m⁻² yr⁻¹ and kgC m⁻² yr⁻², respectively.

3.3.2 Satellite-based GPP products

In this study, GLASS GPP and revised EC-LUE GPP were used as a representative of long-term satellite-based GPP products. As an essential input in the LUE model (Eq. 1), absorbed photosynthetically active radiation (APAR) is a function of LAI, suggesting that LAI is a key parameter in satellite-derived GPP. The GLASS LAI product was used in calculating both the GLASS GPP and the revised EC-LUE GPP products. The spatial correlation coefficient of trends between GLASS LAI and satellite-based GPP (i.e., the mean of GLASS GPP and revised EC-LUE GPP) is 0.42. A previous study over China has shown that satellite-based LAI datasets play a more important role in GPP estimation than meteorological data for all land cover types (Liu et al., 2014). Also, the spatial distribution of trends between LAI and GPP simulated from the DGVM ensemble is even more consistent (r=0.77), confirming the previous studies showing that the trends of GPP and LAI are highly correlated in biome models (Ito et al., 2017; Liu et al., 2019), and the changes of GPP and LAI are consistent in Earth system models from CMIP5 (Hashimoto et al., 2019).

Figure 6 compares the LAI trends of GLASS and the DGVM ensemble during 2001–2015. The spatial distribution of DGVM LAI indicates significant increasing trends over the boreal forest regions, Indonesia, equatorial Africa, and India, and significant decreasing trends over Kazakhstan, Southeast Asia, and western Australia. The trends of GLASS LAI were obviously weaker than those of the DGVM ensemble, especially in the equatorial African region, northern Amazon region, Indonesia, and northern high latitudes. The large inconsistencies between the LAI from the DGVMs and those observed in the satellite products could lead to substantial uncertainties in generating/simulating the global GPP (Figs. 1 and 6). Xiao et al. (2017) suggested that the trends of four satellite-derived LAI products showed large discrepancies in equatorial Africa from 1982 to 2011 and differed across each vegetation type. Jiang et al. (2017) revealed that NOAA satellite orbit changes and MODIS sensor degradation might cause long-term satellite-derived LAI products inconsistent with each other. Xie et al. (2019) also suggested that satellite-derived LAI datasets can cause uncertainties in GPP estimations through model structure and the complexity of the terrain. Hence, the long-term trends of satellite GPP products based on satellite-derived LAI remain highly uncertain (Smith et al., 2016; Jiang et al., 2017; Liu et al., 2018).

3.4 Evaluations at site level

We further adopted 20 sites with observations longer than 15 years from FLUXNET2015 datasets to evaluate long-term GPP trends from global products and simulations. Compared to the site observations, the magnitudes of GPP at most of site locations were underestimated by satellite-based GPP products, FLUXCOM, and the DGVM ensemble mean. Also, the interannual variation of GPP at sites are more obvious than those of global GPP products and NIRv (Fig. 7). This is possibly due to the fact that climate variables at 0.5^∘ grid cells are smoothed out compared to those recorded at individual sites, which leads to a moderate GPP variation.

More than half of the sites indicate that FLUXNET GPP has increased on a long-time scale (Figs. 7, S11), which was mainly caused by rising LUE due to the CO₂ fertilization effect and increased green vegetation cover (Cai and Prentice, 2020). Although FLUXCOM was upscaled from FLUXNET datasets, it did not capture the trends of GPP observed at sites, which was also mentioned by Anav et al. (2015). Sites with significant increasing GPP trends were all captured by DGVM ensemble mean but some were missed by satellite-based GPP (Fig. 7b, d, k, l, m, s, and t). Furthermore, none of sites with decreasing GPP trends were reflected in the global GPP products and NIRv (Fig. 7f, h, i, j, n, q, r), which may be in part due to the different spatial representativeness between a tower fetch and a model or satellite grid point. Therefore, the uncertainty remained when using the site-level observed GPP to evaluate the GPP trends of the DGVM simulations and satellite-based GPP products. We also selected sites with observation longer than 12 years and found similar results (Fig. S12). It is worth mentioning that sites with more than 15 years of observations were all located at NH (from 39.32 to 50.96^∘ N) (Table 2), and only three sites with more than 12 years of observations were located at Trop+SH (Table S2). Therefore, it is hardly possible for us to evaluate the GPP trends of global products over the Trop+SH by using the GPP observations from sites.

4 Conclusions

Based on five kinds of GPP or GPP-related datasets, including satellite-based products, the upscaled product from EC measurements, DGVM simulations, satellite-observed proxy (NIRv), and site-level observations, we comprehensively assessed the global and regional GPP trends during 1982–2015. The simulated spatial pattern of GPP trends from the DGVM ensemble is highly consistent with NIRv but shows considerable inconsistency with satellite-based GPP products, especially in the tropical regions. After 2000, the GPP generated by the satellite-based GPP products decreased significantly in Trop+SH, and the increasing trend in NH also weakened. However, the results of DGVMs showed that global GPP kept increasing after 2000, even in the tropical regions, which was closer to the performance of NIRv. By analyzing the impacts of different drivers in DGVM simulations, the results indicate that the CO₂ fertilization effect has the dominant contribution to the global GPP. Spatially, the CO₂ fertilization effect makes the global GPP increase consistently, while climate has an inhomogeneous impact on GPP trends over different regions.

We further explored the uncertainties in GPP trends among these different datasets. For the DGVM ensemble, globally, the CO₂ fertilization effect causes the largest inter-model spread. At the grid cell level, the uncertainties in simulated GPP trends concentrate over the Trop+SH, which results mainly from climate and LUC. Furthermore, the large discrepancy in the GPP trends between DGVM ensemble and satellite-based GPP products is, to a large extent, induced by the difference of vegetation canopy structure parameter (LAI). Therefore, the highly uncertain satellite-derived LAI data in the tropical regions increase the uncertainty of satellite GPP products and weaken their reliability in explaining changes in the global GPP.

Finally, GPP trends from satellite-based products and DGVM simulations were evaluated by using the FLUXNET2015 dataset. Results show that all of sites with significant increasing GPP trends can be captured by DGVM ensemble mean, but some of them were missed by satellite-based GPP. Also none of sites with decreasing GPP trends were reflected in the global GPP products. Therefore, uncertainty remained when using the FLUXNET observed GPP to evaluate the GPP trends of the global GPP products.

Generally, the differences among models, observations, and products suggest the importance of the research on the GPP trend and make us cautious to interpret the mechanisms of the global carbon cycle by using the long-term GPP products.