Background & Summary

Gross primary production (GPP) is a vital component of the terrestrial carbon budget and plays a prominent role in the global carbon cycle1,2,3,4. Accurate estimation of terrestrial GPP is critical for understanding the interaction between the terrestrial biosphere and the atmosphere in the context of global climate change5,6, projecting future change7, and informing climate policy decisions8. Therefore, characterizing the spatiotemporal variation of GPP9 is a key issue in the climate change study.

GPP is closely related to vegetation types10,11,12, meteorological factors13,14,15,16,17, soil moisture18,19, and other factors. In particular, GPP is affected by vegetation canopy structures12,20,21, e.g., sunlit and shaded leaves. Sunlit leaves can absorb direct and diffuse radiation simultaneously, and light saturation is easy to occur when the radiation is high, while shaded leaves can only absorb diffuse radiation and the absorbed radiation intensity is generally between the light compensation point and the light saturation point22,23,24. The two components of GPP derived from sunlit (GPPsun) and shaded leaves (GPPshade) have drawn increasing attentions recently due to two reasons. First, commonly used sun-induced chlorophyll fluorescence (SIF), which is strongly correlated with GPP at various scales25,26,27,28,29,30, is mainly emitted from sunlit leaves31,32,33,34 and GPPsun can also be used to retrieve key photosynthetic parameters such as maximum carboxylation velocity (Vcmax)35. Besides, the contribution of GPPshade to the total GPP increased with the increases of leaf area index (LAI) and diffuse radiation ratio36,37. Therefore, it is of great importance to distinguish GPPshade and GPPsun respectively for building an improved SIF-GPP relationship and for obtaining high-precision photosynthetic parameters to feed carbon cycle models. Second, some process-based terrestrial biosphere models simulate GPPsun and GPPshade individually, but there is a lack of credible, large-scale, and long-time series GPPsun and GPPshade products for validating model outputs. Thus, such products can support exploring the similarities and differences of sunlit and shaded leaves contributing to GPP or SIF, to further excavate the interior ecological mechanism of different carbon cycle processes and advance carbon cycle modelling.

Light use efficiency (LUE) models have the advantages of few required parameters, and easy to digest remote sensing data38,39. They have been established as a popular method to estimate regional and global carbon fluxes40,41,42,43,44. By incorporating satellite-derived land surface variables into LUE models, several global data products (e.g. MOD17A2) have been produced45,46,47. Considering the difference in solar radiation absorption and LUE of leaves within a canopy, the two-leaf LUE (TL-LUE) model divides the vegetation canopy into shaded and sunlit leaves and calculates GPP separately for these two portions48. Therefore, the TL-LUE model significantly reduces the sensitivity to sky conditions, effectively alleviates the systematic underestimation of GPP under low solar radiation by the MOD17 model, and improves the simulation accuracy of GPP48,49.

In the previous version of the TL-LUE model, the CO2 fertilization effect (CFE), the enhancement of vegetation productivity by the increase of CO2 concentration, is not included. It’s well-known that global atmospheric CO2 concentration has continued to rise, increasing about 17% during 1992–2020, and the increase of atmospheric CO2 substantially enhance global GPP16. CFE on the global terrestrial carbon exchange has attracted widespread attentions4,16,50, but has rarely been considered in LUE models. Thus, it is imperative to include the change of atmospheric CO2 concentration in GPP estimation with LUE models. Based on the previous version of TL-LUE model, this study additionally includes atmospheric CO2 concentration regulation scalar (CS) and replaces temperature regulation scalar (TS) with that used in the Terrestrial Ecosystem Model (TEM) model51 to capture negative effects of low and high temperatures on GPP. Eddy covariance carbon flux measurements of 68 sites (480 site years) and 25 sites (170 site years) from the FLUXNET2015 dataset were used to calibrate and validate the revised TL-LUE model, respectively.

This study employs various remote sensing data as model inputs, including European Space Agency Climate Change Initiative Land Cover (ESA-CCI land cover) from 1992 to 2020 and GLOBMAP leaf area index (GLOBMAP-LAI), in conjunction with meteorological data provided by the Climatic Research Unit and Japanese reanalysis (CRUJRA v2.2) dataset. The GPPs are calculated by the revised TL-LUE model (Fig. 1). The temporal resolutions of the dataset are 8-day, monthly and annual, and the spatial resolution is 0.05°×0.05°. During 1992–2020, the estimated mean annual totals of global GPP, GPPsun, and GPPshade are 125.0 Pg C a−1, 50.5 Pg C a−1, and 74.5 Pg C a−1, respectively.

Fig. 1
figure 1

Workflow of the study. CCI represents European Space Agency Climate Change Initiative Land Cover. CRUJRA represents the Climatic Research Unit and Japanese reanalysis. VPD is vapor pressure deficit. Ta is air temperature. Ts and Cs are regulation scalars for temperature and CO2 concentration. The dswrf is downward solar radiation flux. TL-LUE is two-leaf light use efficiency model. GPPsun and GPPshade are GPP derived by sunlit and shaded leaves.

Full size image

Methods

Model description

The model used in this study is a revised version of the two-leaf light use efficiency (TL-LUE) model48. The revised TL-LUE model adds the atmospheric CO2 concentration regulation scalar and modifies air temperature regulation scalar. GPP is divided into GPPshade and GPPsun52. It is described as Eqs. (1–3):

$${\rm{GPP}}={{\rm{GPP}}}_{{\rm{shade}}}+{{\rm{GPP}}}_{{\rm{sun}}}$$
(1)
$${{\rm{GPP}}}_{{\rm{shade}}}=({{\rm{\varepsilon }}}_{{\rm{msh}}}\times {{\rm{APAR}}}_{{\rm{sh}}})\times {{\rm{T}}}_{{\rm{s}}}\times {{\rm{W}}}_{{\rm{s}}}\times {{\rm{C}}}_{{\rm{s}}}$$
(2)
$${{\rm{GPP}}}_{{\rm{sun}}}=({{\rm{\varepsilon }}}_{{\rm{msu}}}\times {{\rm{APAR}}}_{{\rm{su}}})\times {{\rm{T}}}_{{\rm{s}}}\times {{\rm{W}}}_{{\rm{s}}}\times {{\rm{C}}}_{{\rm{s}}}$$
(3)

where εmsu and εmsh are the maximum light use efficiency of sunlit and shaded leaves (detailed in Table 1), respectively. APARsu and APARsh are the absorbed PAR of sunlit and shaded leaves, respectively. They were expressed as:

$${{\rm{APAR}}}_{{\rm{sh}}}=(1-{\rm{\alpha }})\times \left(\frac{{{\rm{PAR}}}_{{\rm{dif}}}-{{\rm{PAR}}}_{{\rm{dif}},{\rm{u}}}}{{\rm{LAI}}}+{\rm{C}}\right)\times {{\rm{LAI}}}_{{\rm{sh}}}$$
(4)
$${{\rm{APAR}}}_{{\rm{su}}}=(1-{\rm{\alpha }})\times \left(\frac{{{\rm{PAR}}}_{{\rm{dir}}}\times \cos \,{\rm{\beta }}}{\cos \,\theta }+\frac{{{\rm{PAR}}}_{{\rm{dif}}}-{{\rm{PAR}}}_{{\rm{dif}},{\rm{u}}}}{{\rm{LAI}}}+{\rm{C}}\right)\times {{\rm{LAI}}}_{{\rm{su}}}$$
(5)
$${{\rm{LAI}}}_{{\rm{su}}}=2\,\cos \,\theta \times \left(1-{{\rm{e}}}^{-\frac{{\rm{LAI}}\times \Omega }{2\cos \theta }}\right)$$
(6)
$${{\rm{LAI}}}_{{\rm{sh}}}={\rm{LAI}}-{{\rm{LAI}}}_{{\rm{su}}}$$
(7)
$${\rm{C}}=0.07\times \Omega \times {{\rm{PAR}}}_{{\rm{dir}}}\times (1.1-0.1{\rm{LAI}}){{\rm{e}}}^{-\cos {\rm{\theta }}}$$
(8)
$${{\rm{PAR}}}_{{\rm{dif}}}={\rm{PAR}}\times (0.7527+3.8453{\rm{R}}-16.316{{\rm{R}}}^{2}+18.962{{\rm{R}}}^{3}-7.0802{{\rm{R}}}^{4})$$
(9)
$${{\rm{PAR}}}_{{\rm{dir}}}={\rm{PAR}}-{{\rm{PAR}}}_{{\rm{dif}}}$$
(10)
$${{\rm{PAR}}}_{{\rm{dif}},{\rm{u}}}={{\rm{PAR}}}_{{\rm{dif}}}\times {{\rm{e}}}^{-\frac{0.5\Omega {\rm{LAI}}}{\cos \bar{{\rm{\theta }}}}}$$
(11)
$$\cos \,\bar{{\rm{\theta }}}=0.537+0.025\times {\rm{LAI}}$$
(12)

where α is albedo; θ is the solar zenith angle; β is the leaf angle, which is set as 60°; Ω is clumping index (detailed in Table 1); C is multiple scattered radiation (unit: W m−2); PARdir and PARdif (unit: W m−2) are the incoming direct and diffuse photosynthetically active radiation above the canopy; PARdif,u (unit: W m−2) denotes diffuse PAR under the canopy48,53; LAIsu and LAIsh are the LAI of sunlit and shaded leaves; R represents the sky clearness index and equals to S/(S0 cosθ), where S is solar radiation in W m−2, and S0 is the solar constant (1367 W m−2); \(\bar{{\rm{\theta }}}\) is the representative zenith angle for diffuse radiation transmission53.

Table 1 Parameters used in the revised TL-LUE model.
Full size table

The regulation scalars of temperature (Ts)52, water stress (Ws)48, and atmospheric CO2 concentration (Cs)54 are calculated as follows:

$${{\rm{T}}}_{{\rm{s}}}=\frac{({\rm{T}}-{{\rm{T}}}_{{\rm{\max }}})\times ({\rm{T}}-{{\rm{T}}}_{{\rm{\min }}})}{({\rm{T}}-{{\rm{T}}}_{{\rm{\max }}})\times ({\rm{T}}-{{\rm{T}}}_{{\rm{\min }}})-{({\rm{T}}-{{\rm{T}}}_{{\rm{opt}}})}^{2}}$$
(13)
$${{\rm{W}}}_{{\rm{s}}}=\frac{{{\rm{VPD}}}_{{\rm{\max }}}-{\rm{VPD}}}{{{\rm{VPD}}}_{{\rm{\max }}}-{{\rm{VPD}}}_{{\rm{\min }}}}$$
(14)
$${{\rm{C}}}_{{\rm{s}}}=\frac{{{\rm{C}}}_{{\rm{i}}}-{\Gamma }^{* }}{{{\rm{C}}}_{{\rm{i}}}+2{\Gamma }^{* }}$$
(15)

where the maximum (Tmax) and minimum temperatures for vegetation photosynthesis (Tmin) were set to 313.15 K and 273.15 K, respectively49. The optimum temperature for vegetation photosynthesis (Topt) is set as the average of different types in Huang et al.55 (details in Table 1). VPDmax and VPDmin are the VPD when GPP reaches the maximum and minimum, respectively38. If the value of VPD is greater than or equal to VPDmax, Ws is equal to 0 and if the value of VPD is less than or equal to VPDmin, Ws is set to 148. Γ* is the CO2 compensation point in the absence of dark respiration (calculated by Eq. 16)54; Ci is the intercellular concentration of CO2 (ppm).

$${\Gamma }^{* }=4.22{\times {\rm{e}}}^{\frac{37830({\rm{T}}-298.15)}{298.15{\rm{RT}}}}$$
(16)
$${{\rm{C}}}_{{\rm{i}}}={{\rm{C}}}_{{\rm{a}}}\times {\rm{\chi }}$$
(17)
$${\rm{\chi }}=\frac{{\rm{\xi }}}{{\rm{\xi }}+\sqrt{{\rm{VPD}}}}$$
(18)
$${\rm{\xi }}=\sqrt{\frac{356.51{\rm{K}}}{1.6{{\rm{\eta }}}^{\ast }}}$$
(19)
$${\rm{K}}={{\rm{K}}}_{{\rm{c}}}\times \left(1+\frac{{{\rm{P}}}_{{\rm{o}}}}{{{\rm{K}}}_{{\rm{o}}}}\right)$$
(20)
$${{\rm{K}}}_{{\rm{c}}}=39.97\times {{\rm{e}}}^{\frac{79.43\times ({\rm{T}}-298.15)}{298.15{\rm{RT}}}}$$
(21)
$${{\rm{K}}}_{{\rm{o}}}=27480\times {{\rm{e}}}^{\frac{36.38\times ({\rm{T}}-298.15)}{298.15{\rm{RT}}}}$$
(22)

where Ca is the atmospheric CO2 concentration (using NOAA global monthly mean CO2 concentration at the unit of ppm), and χ is the ratio of leaf-internal to ambient CO256. K is the Michaelis-Menten coefficient of Rubisco and η* is the viscosity of water relative to its value at 25 °C (0.8903)57. Kc and Ko are the Michaelis-Menten coefficients of Rubisco for CO2 and O2, respectively, and Po is the partial pressure of O2 (21 kPa)56. R is the molar gas constant (8.314 J mol−1 K−1).

Input data and processing

Eddy covariance flux measurements from the FULLSET daily product in the FLUXNET2015 dataset (https://fluxnet.org) were used for model calibration and validation. We selected site years data according to the following requirements: the missing observations of air temperature (TA_F), vapor pressure deficit (VPD_F), CO2 mole fraction (CO2_F_MDS), incoming photosynthetic photon flux density (PPFD_IN), or shortwave radiation (SW_IN_F), gross primary production (GPP_DT_CUT_MEAN) in one year are less than 2 months. A linear interpolation was applied to fill the individual missing values, which accounted for about 2% of the total measurements. About 75% of the sites were randomly selected to calibrate the revised TL-LUE model parameters for each vegetation type, and the remaining sites were applied for parameters validation. The sites and years selected for calibration and validation are detailed in Table S1. The shortwave-to-PAR conversion parameter in global was estimated to vary between 0.39 to 0.5358,59,60. With the comparison between measurements of PAR and shortwave radiation in sites, 0.43 is most suitable for this study. The spatial distribution of calibration and validation sites is shown in Figure S1.

The land cover dataset we used is European Space Agency Climate Change Initiative Land Cover (ESA-CCI land cover) at a 300 m spatial resolution for every year from 1992 to 2020 (https://cds.climate.copernicus.eu/cdsapp#!/dataset/satellite-land-cover). We resampled the raw data to 0.05°×0.05° using nearest neighbour resampling. ESA-CCI land cover dataset uses the United Nations Land Cover Classification System (LCCS), thus, we converted them to match the International Geosphere-Biosphere Program land cover scheme (IGBP)61. In particular, ESA-CCI land cover provides land cover data for the years before 2000, which makes it possible to study changes in GPP caused by changes in land cover types over long-term series.

GLOBMAP leaf area index (GLOBMAP-LAI)62 as a model input is available at 0.0727° spatial resolution for every 8 days (2001–2020) and half-month (1992–2000) from 1992 to 2020. GLOBMAP LAI (Version 3) provides a consistent long-time LAI product (1981–2020, continuously updated) by quantitative fusion of Moderate Resolution Imaging Spectroradiometer (MODIS) and historical Advanced Very High Resolution Radiometer (AVHRR) data. The long-term LAI series was made up by combination of AVHRR LAI (1981–2000) and MODIS LAI (2001–2020). MODIS LAI series was generated from MODIS land surface reflectance data (MOD09A1)63 based on the GLOBCARBON LAI algorithm64. The relationships between GIMMS NDVI and MODIS LAI were established pixel by pixel using the two data series during overlapped period (2000–2006). And the AVHRR LAI65 back to 1981 was estimated from historical AVHRR observations based on these pixel-level relationships. The clumping effects was considered at the pixel level by employing global clumping index map at 500 m resolution66. The cloud mask for the MOD09A1 data were created by a new cloud detection algorithm based on time series surface reflectance observations67. GLOBMAP-LAI has been smoothed by locally adjusted cubic spline capping approach68. We resampled them to 0.05°×0.05° for the model. Additionally, we extracted the LAI of each site from GLOBMAP-LAI (500 m, 8-day) for model calibration and validation.

The Climatic Research Unit and Japanese reanalysis (CRUJRA) version 2.2 data69 provided 6-hourly at 0.5° resolution meteorological variables, such as downward solar radiation flux (dswrf, unit: J m−2 6 h−1), specific humidity(spfh, unit: kg kg−1), the temperature at 2 m (tmp, unit: K), pressure (pres, unit: Pa), from 1992 to 2020. Daily dswrf used in this study (unit: J m−2 d−1) was converted from the CRUJRA dswrf dataset by summing four 6-hourly data per day. Temperature (unit: K) and vapor pressure deficit (VPD, unit: kPa) by taking the mean of the 6-hourly data of each day in CRUJRA. The daily dswrf was corrected with site shortwave radiation data. The global monthly CO2 concentration (unit: ppm) is available on www.esrl.noaa.gov/gmd/ccgg/trends/. VPD was calculated using specific humidity, pressure, and temperature as Eqs. (23, 24):

$${\rm{VPD}}={{\rm{VP}}}_{{\rm{sat}}}-{\rm{spfh}}\times \frac{{\rm{pres}}/1000}{0.622+{\rm{spfh}}\times 0.378}$$
(23)
$${{\rm{VP}}}_{{\rm{sat}}}=0.61121\times {{\rm{e}}}^{\left(18.678-\frac{{{\rm{t}}}_{{\rm{air}}}}{234.5}\right)\times \frac{{{\rm{t}}}_{{\rm{air}}}}{257.14+{{\rm{t}}}_{{\rm{air}}}}}$$
(24)

where tair is the air temperature at the unit of °C. VPsat means the saturated vapor pressure (kPa); spfh represents specific humidity(kg kg−1); press is pressure(Pa).

Calibration of model parameters

The maximum LUEs of the sunlit (εmsu) and shaded (εmsh) leaves spatially differ due to the changes in vegetation canopy structures and vegetation species types70, which leads to the distinct εmsh and εmsu of different vegetation types. To reduce the GPP simulation bias caused by εmsh and εmsu of sunlit and shaded leaves, the daily data of 68 sites (480 site years) in the FLUXNET2015 dataset (details in Table S1) were used for parameter optimization with the shuffled complex Evolution-University of Arizona method49,71. The agreement index (d) was used as the optimization criterion. This index was widely used in parameter optimization49,72. It ranges from 0 (complete disagreement) to 1 (complete agreement). Parameter values identified when d maximizes are optimization results. The calculation of d is as Eq. (25):

$${\rm{d}}=1-{\sum }_{{\rm{n}}=1}^{{\rm{m}}}{\left({{\rm{E}}}_{{\rm{n}}}-{{\rm{M}}}_{{\rm{n}}}\right)}^{2}/{\sum }_{{\rm{n}}=1}^{{\rm{m}}}{\left(| {{\rm{E}}}_{{\rm{n}}}-\bar{{\rm{M}}}| +| {{\rm{M}}}_{{\rm{n}}}-\bar{{\rm{M}}}| \right)}^{2}$$
(25)

where m is the number of all measurements; En and Mn are the nth estimation and measured GPP, respectively. \(\bar{{\rm{M}}}\) represents the mean of all measured GPP values.

The average εmsh and εmsu values with standard deviation of all vegetation types are shown in Table 1. The parameters of deciduous needleleaf forests (DNF) were consistent with DBF settings in modelling. The R2 of observed GPP (GPPEC) and estimated GPP (GPP) in WSA and SAV were 0.39 and 0.41, respectively, and the R2 values of other vegetation types were all above 0.5, in which the R2 of DBF was the highest (0.89). The calibration results of different vegetation types are shown in Figure S2.

Data Records

The dataset provides global gridded 0.05° GPP, GPPshade and GPPsun at three temporal resolutions (8-day, monthly, annual) from 1992 to 2020. The units are g C m−2 (8day)−1, g C m−2 month−1 and g C m−2 a−1, respectively. We divided the dataset into 29 folders by year, where each folder contains the data for the year at three temporal scales. The files were named as “GGGG_v21_TTTT.tif” and stored in the GeoTiff format, where GGGG represents GPP, GPPshade and GPPsun. For the 8-day scale, TTTT represents year and DOY (e.g. Shade_GPP_v21_1999_249.tif). For the one-month scale, TTTT represents year and month (e.g. Sun_GPP_v21_1999_01.tif). For the annual scale, TTTT represents only year (e.g. GPP_v21_1999.tif). The scale factor of the monthly data is 0.1, that of the 8-day data is 0.01. The data type of the monthly and 8-day data is 16-bit integer, and that of the annual data is double. All datasets are publicly available from the DRYAD repository https://doi.org/10.5061/dryad.dfn2z352k73.

The high annual average GPP values were mainly distributed in the Amazon, Indonesia and Congo Basin in the low latitude regions, which was about 3500 g C m−2 a−1. The relatively high GPP (~2000 g C m−2 a−1) occurred in Southeast Asia, the southeast of North and South America, and south central Africa, and the GPP that ranged from 0 to 1500 g C m−2 a−1 accounted for 74.7% of global vegetation cover. The shaded and sunlit GPP (GPPshade and GPPsun) had similar spatial distribution with GPP, but the value of GPPsun was lower than GPPshade near the equator (Fig. 2b,c). GPPshade that ranged from 0 to 1000 g C m−2 a−1 accounted for 78.7% of global vegetation cover, and GPPsun that ranged from 0 to 750 g C m−2 a−1 accounted for 80.8% of global vegetation cover.

Fig. 2
figure 2

Spatial distribution of the annual average (a) GPP, (b) GPPshade, and (c) GPPsun from 1992 to 2020.

Full size image

From 1992 to 2020, the GPP of Malaysia, Southeast Asia, Indian Peninsula, central Africa, north and southeast South America, and central Europe all showed an obvious increasing trend, and the rate was close to 20 g C m−2 a−2. The GPP scattered in central South America, east Africa, Central Asia, and Southeast Asia showed a significant decreasing trend, and the change rate was near to −20 g C m−2 a−2 (Fig. 3a). For GPPshade, eastern and southern Asia, central Europe, central and western Africa, northwest and southeast South America showed a significant increasing trend, and the change rate was above 10 g C m−2 a−2, which was consistent with the trend of GPP. Besides, areas with reduced GPPshade (trend < −10 g C m−2 a−2) were sporadically distributed in Central Asia, Central South America, Eastern Africa, and Southeast Asia (Fig. 3b). For GPPsun, central and southeastern South America, southeastern Asia, and Europe showed the most obvious increase, with a rate of change of around 10 g C m−2 a−2, which was similar to the spatial distribution of GPP and GPPshade trends. However, the value of GPPsun was lower than GPPshade (Fig. 3b,c). The vast majority of global vegetation cover with significant change exhibits an increasing trend (GPP trend >0: 92.0%, GPPshade trend >0: 91.2%, and GPPsun trend >0: 88.7%). In general, most global GPP, GPPshade, and GPPsun revealed increasing trends (Fig. 3).

Fig. 3
figure 3

Spatial distribution of the trend of (a) GPP, (b) GPPshade, and (c) GPPsun during 1992 to 2020. The results have removed the value which is not significant (p > 0.05).

Full size image

The mean annual totals of global GPP, GPPsun, and GPPshade from 1992 to 2020 are 125.0 Pg C a−1, 50.5 Pg C a−1 and 74.5 Pg C a−1, respectively. Overall, the GPP proportions of individual vegetation types to total were similar for global GPP, GPPsun, and GPPshade. Among the 11 vegetation types, EBF contributed most GPP, followed by CRO. These two types accounted for more than half of the total GPP, GPPsun, and GPPshade (Fig. 4). The GPP of SAV, MF, WET, and WSA were 1.0 Pg C a−1, 2.0 Pg C a−1, 2.8 Pg C a−1, and 2.9 Pg C a−1, respectively, which were relatively low (Fig. 4a). The GPPsun were close to GPPshade for CRO, GRA, and WET, while the GPPsun were higher than GPPshade for SAV and WSA. The GPPshade of ENF, EBF, DNF, DBF, MF, and OSH with relatively complicated canopy structures were much higher than their GPPsun (Fig. 4b,c). Overall, GPPshade played a key role in GPP for forest types, while GPPsun was greater than GPPshade for non-forest types. In total, GPPshade contributes more to GPP than GPPsun. In addition, the GPP of all vegetation types showed an increasing trend. With one exception that the GPPsun of SAV showed a decreasing trend, the GPPsun and GPPshade of other vegetation types showed an increasing trend, among which GPP (also GPPsun, and GPPshade) of CRO showed the distinctively greatest increasing trend. Except WSA, the increasing trend of GPPshade of other vegetation types is greater than that of GPPsun (Fig. 4d). It’s worth noting that for DBF, the increasing trend of GPP is relatively large, and GPPshade overwhelmingly outweighs GPPsun. MF shows a similar phenomenon, but the overall trend is smaller.

Fig. 4
figure 4

The mean annual totals of (a) GPP, (b) GPPsun, (c) GPPshade and (d) the trend of GPPshade and GPPsun from 1992 to 2020 for different vegetation types over the globe.

Full size image

Technical Validation

Validation of model parameters

Carbon flux data of 25 sites (170 site years) from the FLUXNET2015 dataset were selected (detailed in Table S1) for model validation. The comparison between GPP that estimated by the revised TL-LUE model with optimized εmsu and εmsh, and GPP measurements (GPPEC) at each flux site is shown in Fig. 5. The revised TL-LUE model performed well in estimation of GPP for all vegetation types. All sites have R2 above 0.5, except for AU-Gin (WSA) (R2 = 0.49) and BR-Sa3 (EBF) (R2 = 0.36).

Fig. 5
figure 5

Validation of daily GPP estimated by the revised TL-LUE (GPP) model with tower measurements (GPPEC) at 25 FLUXNET sites. The revised TL-LUE model was driven by average optimized parameters for different vegetation types, tower-based meteorological data, and smoothed LAI.

Full size image

Comparisons with other global GPP products

Previous studies have shown that the differences in GPP estimation are usually caused by different model structures, parameter settings, and input data74,75,76. Here, we compare our global annual GPP with several global GPP products derived from remote sensing models, including data-driven models and LUE models (Fig. 6).

Fig. 6
figure 6

Comparison of global annual GPP totals estimated from a set of global remote sensing GPP products.

Full size image

For data-driven GPP products, Li et al.77 used SIF-GPP relationship and SIF observed by the Orbiting Carbon Observatory-2 (OCO-2) to obtain global GPP, which named GOSIF GPP (135.5 ± 8.8 Pg C a−1, 2001–2017). And the WECANN product was produced using an artificial neural network (ANN) with SIF and other data sources as inputs78. The WECANN-GPP showed an obvious decreasing trend between 2007 and 2015, with a range of 110.1 to 118.2 Pg C a−1. FluxSat GPP79, which was generated using satellite data-driven models based on the LUE framework, ranged from 148.1 to 156.7 Pg C a−1 during 2001 to 2019. In addition, a few recent studies have identified that there is a strong spatio-temporal correlation between near-infrared reflectance (NIRv) and GPP, suggesting that NIRv can be employed to estimate the GPP of vegetation13,79. The global total of GPP estimated from NIRv by Wang et al.80 was 130.5 ± 3.2 Pg C a−1 during 1992–2018, close to our estimate. For the LUE models, the range of GPP obtained by the improved EC-LUE model81 was 106.4~118.3 Pg C a−1 (1992–2018), and that by VPM model was 131.2~140.0 Pg C a−1 (2000–2016). The range of MOD17A2H.00682 was 94.8~106.2 Pg C a−1 (2001–2020). The global annual GPP of FluxSat was the highest, and that of MOD17A2 was the lowest. Our estimated global GPP ranged from 120.02 to 132.65 Pg C a−1, placing at the middle of the various GPP products. Anav et al.76 suggested that previous observation-constrained estimates of global GPP45,83 (e.g. based on either δ18O measurements of atmospheric CO2 or eddy covariance flux upscaling) was at 120 Pg C a−1 for the period before 2010. Thus, our estimate agreed reasonably well with such estimates.

In addition, our estimated global GPP showed an overall increasing trend, which is consistent with most other GPP products. Only WECANN-GPP and EC-LUE GPP (after 2000) showed a significant declining trend. The declining trend of WECANN obviously associated with the degradation of GOME-2 SIF sensor, and the GOME-2 SIF data used for training was not corrected84.

Uncertainties

Previous studies have indicated that different LAI products lead to clear differences in estimated GPP85, and the uncertainties of various LAI products are higher in low-latitude areas86,87. Since the land cover used by GLOBMAP-LAI is different from that used in this study, the corresponding LAI values of a small amount pixels (<0.01%) for SAV and WSA in low latitude areas are relatively high, which lead to abnormally high estimated GPP. These anomalies are not processed because of remaining the authenticity of data.

The model parameters varied in different areas, due to the plant species included in the same vegetation type. In particular, the LUE of C3 and C4 crops was greatly discrepant as many previous studies proved88,89,90,91. According to the restriction of ESA-CCI land cover data, C3 and C4 plants could not be distinguished, so the GPP of C3 and C4 crops cannot be gained separately. The LAI of each site extracted from GLOBMAP-LAI with a resolution of 500 m cannot completely match the flux tower data scale, which cause uncertainty in the parameter optimization. Simultaneously, the eddy covariance measurements also have some uncertainties, which inevitably affected the parameter optimization.

Potential benefits and usages of this dataset

Facilitate researches on the relationship between SIF and GPP

As is known, SIF signals come mainly from sunlit leaves31,92. Most of the current studies on the relationship between SIF and the photosynthesis of sunlit leaves are at canopy and leaf scales, while similar studies at large scales are currently rare due to the lack of publicly available global or regional GPP datasets for sunlit leaves. In addition, it is known that there is a link between SIF and GPP across biome types34, but the relationship is not well quantified. The GPPsun and GPP provided in this study may help to explore the relationship between SIF and GPPsun or GPP in different ecosystem types at large scales.

Facilitate researches on the interactions between solar radiation and terrestrial carbon cycling

Compared with direct radiation, the increase of diffuse radiation can effectively promote carbon fixation93. Shaded leaves make more effective use of diffuse radiation48, and GPPshade plays a major role in the vegetation areas with more sheltered leaves or cloudy conditions94. The GPPshade and GPPsun cannot be directly measured over large regions. The GPPshade and GPPsun estimated by the revised TL-LUE model can capture the contribution of sunlit and shaded leaves to GPP in long-term and at large scales, and make it possible to quantify the carbon fixation increase (or decrease) influenced by the change in radiation fraction over a long period and at large scales, which facilitates further investigations on the interactions between radiation and terrestrial carbon cycling.

Facilitate researches on the dynamics, processes and drivers of GPP at large scales

This study provides 8-day and monthly GPP from 1992 to 2020, which allows for studying the changes in seasonal cycles (e.g. amplitude and phase changes of growing season) of GPP and its processes (GPPsun and GPPshade) over many years. In addition, we employed a long-term ESA-CCI land cover data since 1992 (while most other dataset uses MODIS land cover since 2000) with the consideration of year to year land cover changes, enabling it to characterize the impact of land cover change on GPP. The dataset would help to dig the underlying mechanisms of climate and human impacts on global terrestrial GPP.

Usage Notes

In the dataset, in order to ensure the authenticity, we did not delete or modify a small number of abnormally high values. Therefore, when using this dataset, you can set thresholds to remove the anomalies.