A physiology‐based Earth observation model indicates stagnation in the global gross primary production during recent decades

Abstract Earth observation‐based estimates of global gross primary production (GPP) are essential for understanding the response of the terrestrial biosphere to climatic change and other anthropogenic forcing. In this study, we attempt an ecosystem‐level physiological approach of estimating GPP using an asymptotic light response function (LRF) between GPP and incoming photosynthetically active radiation (PAR) that better represents the response observed at high spatiotemporal resolutions than the conventional light use efficiency approach. Modelled GPP is thereafter constrained with meteorological and hydrological variables. The variability in field‐observed GPP, net primary productivity and solar‐induced fluorescence was better or equally well captured by our LRF‐based GPP when compared with six state‐of‐the‐art Earth observation‐based GPP products. Over the period 1982–2015, the LRF‐based average annual global terrestrial GPP budget was 121.8 ± 3.5 Pg C, with a detrended inter‐annual variability of 0.74 ± 0.13 Pg C. The strongest inter‐annual variability was observed in semi‐arid regions, but croplands in China and India also showed strong inter‐annual variations. The trend in global terrestrial GPP during 1982–2015 was 0.27 ± 0.02 Pg C year−1, and was generally larger in the northern than the southern hemisphere. Most positive GPP trends were seen in areas with croplands whereas negative trends were observed for large non‐cropped parts of the tropics. Trends were strong during the eighties and nineties but levelled off around year 2000. Other GPP products either showed no trends or continuous increase throughout the study period. This study benchmarks a first global Earth observation‐based model using an asymptotic light response function, improving simulations of GPP, and reveals a stagnation in the global GPP after the year 2000.


S2 Order of the environmental constraints on GPP
To determine in which order the constraining environmental variables (Tair, CO2, SWC, and VPD) were to be introduced in the model for the different biomes, we used a regression tree analysis. It is a robust statistical tool to analyse complex, nonlinear relationships and interactions between a single response variable and several explanatory variables (De'ath & Fabricius, 2000). The data set was repeatedly split into more and more homogeneous subgroups, each categorized by values of both the dependent and the independent variables. Splitting continues until a tree is created, which is then pruned back to a proper size according a cross validation procedure. In each analysis, we separated the data into ten subgroups of approximately equal size. Each subgroup was left out once, and ten trees were created with the remaining nine subgroups, which were evaluated against the left-out subgroup. The error was summed for all ten trees and for each of the different tree sizes. The smallest tree with the minimum error was selected. The cross-validation was repeated 100 times and the most common tree size was used in the final analysis. We repeated the analysis allowing various number of days in each subgroup (bin size) for further splitting the tree, and then used the bin size generating the lowest RMSE in the final analysis. The analysis picks out the environmental variables which makes it easiest to place GPP within a certain easily defined range. It was therefore used for deciding the order in which the environmental constraining variables were applied in the GPP model for the different biomes (Table S2). Table S2. Statistics for the regression tree analysis studying relationships between seasonal dynamics in the gross primary production and the explanatory variables for the various biomes. The R 2 is the coefficient of determination and RMSE is the root-mean-square-error for the regression trees. Tair is air temperature at 2 m height (°C), NDVI is the normalized difference vegetation index, PAR is daily mean photosynthetically active radiation (W m −2 ), SWC is soil water content at 0.07 m soil depth (% volumetric water content), and VPD is vapour pressure deficit (kPa). The values in brackets are relative predictor importance. The bin size is minimum data size for each subgroup, and pruning level is the number of splits of the regression tree. S3 Impact of meteorological and hydrological stress on GPP S3.1 Evergreen needleleaf forest Figure S1 Relationships between GPP variability and environmental variables for evergreen needleleaf forest. a-d) Relationships between direct measures of GPP and the environmental variables, and e-h) relationship between the ratio between modelled and measured GPP (the GPP scalar) and the environmental variables. The environmental variables are: air temperature at 2 m height (Tair; °C), atmospheric CO2 concentrations, soil water content at 0.0-0.07 m soil depth (SWC %volumetric water content), and vapour pressure deficit (VPD; kPa). The variables follow the order of introduction into the model based on the regression tree analysis (Table S2). Included are the impact on the upper threshold (98 th percentile; red) and the median (blue). Figure S2 Relationships between GPP variability and environmental variables for evergreen broadleaf forest. a-d) Relationships between direct measures of GPP and the environmental variables, and e-h) relationship between the ratio between modelled and measured GPP (the GPP scalar) and the environmental variables. The environmental variables are: air temperature at 2 m height (Tair; °C), atmospheric CO2 concentrations, soil water content at 0.0-0.07 m soil depth (SWC %volumetric water content), and vapour pressure deficit (VPD; kPa). The variables follow the order of introduction into the model based on the regression tree analysis (Table S2). Included are the impact on the upper threshold (98 th percentile; red) and the median (blue).

S3.3 Deciduous needleleaf forest
Deciduous needleleaf forest is mainly larch forest located at the northern hemisphere in the Siberian tundra ( Fig S1). It is the biome with least data in FLUXNET 2015 data base generating a highly uncertain analysis for this biome.

Figure S3
Relationships between GPP variability and environmental variables for deciduous needleleaf forest. a-d) Relationships between direct measures of GPP and the environmental variables, and e-h) relationship between the ratio between modelled and measured GPP (the GPP scalar) and the environmental variables. The environmental variables are: air temperature at 2 m height (Tair; °C), atmospheric CO2 concentrations, soil water content at 0.0-0.07 m soil depth (SWC %volumetric water content), and vapour pressure deficit (VPD; kPa). The variables follow the order of introduction into the model based on the regression tree analysis (Table S2). Included are the impact on the upper threshold (98 th percentile; red) and the median (blue). Figure S4 Relationships between GPP variability and environmental variables for deciduous needleleaf forest. a-d) Relationships between direct measures of GPP and the environmental variables, and e-h) relationship between the ratio between modelled and measured GPP (the GPP scalar) and the environmental variables. The environmental variables are: air temperature at 2 m height (Tair; °C), atmospheric CO2 concentrations, soil water content at 0.0-0.07 m soil depth (SWC %volumetric water content), and vapour pressure deficit (VPD; kPa). The variables follow t the order of introduction into the model based on the regression tree analysis (Table S2). Included are the impact on the upper threshold (98 th percentile; red) and the median (blue). Figure S5 Relationships between GPP variability and environmental variables for mixed forest. ad) Relationships between direct measures of GPP and the environmental variables, and e-h) relationship between the ratio between modelled and measured GPP (the GPP scalar) and the environmental variables. The environmental variables are: air temperature at 2 m height (Tair; °C), atmospheric CO2 concentrations, soil water content at 0.0-0.07 m soil depth (SWC %volumetric water content), and vapour pressure deficit (VPD; kPa). The variables follow the order of introduction into the model based on the regression tree analysis (Table S2). Included are the impact on the upper threshold (98 th percentile; red) and the median (blue). Figure S6 Relationships between GPP variability and environmental variables for savanna/shrublands. a-d) Relationships between direct measures of GPP and the environmental variables, and e-h) relationship between the ratio between modelled and measured GPP (the GPP scalar) and the environmental variables. The environmental variables are: air temperature at 2 m height (Tair; °C), atmospheric CO2 concentrations, soil water content at 0.0-0.07 m soil depth (SWC %volumetric water content), and vapour pressure deficit (VPD; kPa). The variables follow the order of introduction into the model based on the regression tree analysis (Table S2). Included are the impact on the upper threshold (98 th percentile; red) and the median (blue). Figure S7 Relationships between GPP variability and environmental variables for grasslands a-d) Relationships between direct measures of GPP and the environmental variables, and e-h) relationship between the ratio between modelled and measured GPP (the GPP scalar) and the environmental variables. The environmental variables are: air temperature at 2 m height (Tair; °C), atmospheric CO2 concentrations, soil water content at 0.0-0.07 m soil depth (SWC %volumetric water content), and vapour pressure deficit (VPD; kPa). The variables follow the order of introduction into the model based on the regression tree analysis (Table S2). Included are the impact on the upper threshold (98 th percentile; red) and the median (blue). Figure S8 Relationships between GPP variability and environmental variables for grasslands. a-d) Relationships between direct measures of GPP and the environmental variables, and e-h) relationship between the ratio between modelled and measured GPP (the GPP scalar) and the environmental variables. The environmental variables are: air temperature at 2 m height (Tair; °C), atmospheric CO2 concentrations, soil water content at 0.0-0.07 m soil depth (SWC %volumetric water content), and vapour pressure deficit (VPD; kPa). The variables follow the order of introduction into the model based on the regression tree analysis (Table S2). Included are the impact on the upper threshold (98 th percentile; red) and the median (blue).

S4 Constraining GPP based on environmental stress
For evergreen needleleaf forest, deciduous broadleaf forest, and mixed forest, GPP was modelled as: Where GPP is the final modelled GPP product, GPPopt is the modelled GPP based on the light response function with NDVI-modelled Fopt and α (Eq. 1-3 in main text). STair is a scalar constraining GPP at low air temperature modelled using Eq. 4 in main text. However, air temperature can only constrain GPP, never enhance it and as such STair larger than one was set to one. The TSWC is the threshold of the upper GPP boundary for a certain soil water condition (SWC) condition, and it was modelled based on Eq. 5 in main text. GPP at a certain SWC can never reach higher values than as estimated by the upper boundary, and modelled GPP that was larger was therefore set to TSWC.
For evergreen broadleaf forest and croplands, there were only strong relationships between GPP and air temperature and GPP for these biomes were therefore modelled as: Strong relationships between upper threshold of GPP and air temperature and SWC were seen for deciduous needleleaf forest, and GPP for this biome was modelled as: Where TTair is the threshold of the upper GPP for a certain air temperature, and it was modelled based on Eq. 4 in main text.
For savanna/shrublands and grasslands, we additionally saw strong relationships between the upper boundary of GPP and vapour pressure deficit, and modelled GPP was therefore limited to never reach above this threshold: Where TVPD is the upper GPP threshold for a certain vapour pressure deficit condition, and it was modelled based on Eq. 4 in main text.

S5 Model sensitivity analysis
S5.1 Sensitivity of LRF-GPP to Fopt and α variability Figure S9 Sensitivity of LRF-GPP to variability in optimized carbon uptake at light saturation (Fopt) and quantum efficiency (α) for evergreen needle leaf forest. Impact on a) GPP; b) bias (modelled GPP subtracted from field measured); c) root-mean-square-error (RMSE); d) coefficient of determination (R 2 ); and e) slope and f) intercept from an ordinary least square linear regression fitted between modeled and field measured GPP. Note that scales differ for the different biomes ( Fig. S9-S16).

Figure S10
Sensitivity of the LRF-GPP to variability in optimized carbon uptake at light saturation (Fopt) and quantum efficiency (α) for evergreen broadleaf forest. Impact on a) GPP; b) bias (modelled GPP subtracted from field measured); c) root-mean-square-error (RMSE); d) coefficient of determination (R 2 ); and e) slope and f) intercept from an ordinary least square linear regression fitted between modeled and field measured GPP. Note that scales differ for the different biomes ( Fig. S9-S16).
Figure S11 Sensitivity of LRF-GPP to variability in optimized carbon uptake at light saturation (Fopt) and quantum efficiency (α) for deciduous needleleaf forest. Impact on a) GPP; b) bias (modelled GPP subtracted from field measured); c) root-mean-square-error (RMSE); d) coefficient of determination (R 2 ); and e) slope and f) intercept from an ordinary least square linear regression fitted between modeled and field measured GPP. Note that scales differ for the different biomes ( Fig. S9-S16).
Figure S12 Sensitivity of LRF-GPP to variability in optimized carbon uptake at light saturation (Fopt) and quantum efficiency (α) for deciduous broadleaf forest. Impact on a) GPP; b) bias (modelled GPP subtracted from field measured); c) root-mean-square-error (RMSE); d) coefficient of determination (R 2 ); and e) slope and f) intercept from an ordinary least square linear regression fitted between modeled and field measured GPP. Note that scales differ for the different biomes ( Fig. S9-S16).

Figure S13
Sensitivity of the LRF-GPP to variability in optimized carbon uptake at light saturation (Fopt) and quantum efficiency (α) for mixed forest. Impact on a) GPP; b) bias (modelled GPP subtracted from field measured); c) root-mean-square-error (RMSE); d) coefficient of determination (R 2 ); and e) slope and f) intercept from an ordinary least square linear regression fitted between modeled and field measured GPP. Note that scales differ for the different biomes ( Fig. S9-S16).
Figure S14 Sensitivity of LRF-GPP to variability in optimized carbon uptake at light saturation (Fopt) and quantum efficiency (α) for savanna/shrublands. Impact on a) GPP; b) bias (modelled GPP subtracted from field measured); c) root-mean-square-error (RMSE); d) coefficient of determination (R 2 ); and e) slope and f) intercept from an ordinary least square linear regression fitted between modeled and field measured GPP. Note that scales differ for the different biomes ( Fig. S9-S16).
Figure S15 Sensitivity of LRF-GPP to variability in optimized carbon uptake at light saturation (Fopt) and quantum efficiency (α) for grasslands. Impact on a) GPP; b) bias (modelled GPP subtracted from field measured); c) root-mean-square-error (RMSE); d) coefficient of determination (R 2 ); and e) slope and f) intercept from an ordinary least square linear regression fitted between modeled and field measured GPP. Note that scales differ for the different biomes ( Fig. S9-S16).
Figure S16 Sensitivity of LRF-GPP to variability in optimized carbon uptake at light saturation (Fopt) and quantum efficiency (α) for croplands. Impact on a) GPP; b) bias (modelled GPP subtracted from field measured); c) root-mean-square-error (RMSE); d) coefficient of determination (R 2 ); and e) slope and f) intercept from an ordinary least square linear regression fitted between modeled and field measured GPP. Note that scales differ for the different biomes ( Fig. S9-S16).
S5.2 Sensitivity of the LRF-GPP to STair and SVPD variability Figure S17 Sensitivity of LRF-GPP to variability in environmental scalars (air temperature (Tair); soil water content (SWC)) for evergreen needleleaf forest. The percentiles of each scalar were extracted and the constraining models were fitted. Shown is impact of this percentile fitting on a) GPP; b) bias (modelled GPP subtracted from measured); c) root-mean-square-error (RMSE); d) coefficient of determination (R 2 ); and e) slope and f) intercept from an ordinary least square linear regression fitted between modeled and measured GPP. Note that scales differ for the different biomes ( Fig. S17-S26).
Figure S18 Sensitivity of LRF-GPP to variability in the air temperature scalar for evergreen broadleaf forest. The percentiles of each scalar were extracted and the constraining models were fitted. Shown is impact of this percentile fitting on a) GPP; b) bias (modelled GPP subtracted from measured); c) root-mean-square-error (RMSE); d) coefficient of determination (R 2 ); and e) slope and f) intercept from an ordinary least square linear regression fitted between modeled and field measured GPP. Note that scales differ for the different biomes ( Fig. S17-S26).
Figure S19 Sensitivity of LRF-GPP to variability in the air temperature scalar for deciduous broadleaf forest. The percentiles of each scalar were extracted and the constraining models were fitted. Shown is impact of this percentile fitting on a) GPP; b) bias (modelled GPP subtracted from field measured); c) root-mean-square-error (RMSE); d) coefficient of determination (R 2 ); and e) slope and f) intercept from an ordinary least square linear regression fitted between modeled and field measured GPP. Note that scales differ for the different biomes ( Fig. S17-S26).

Figure S20
Sensitivity of LRF-GPP to variability in air temperature scalar for mixed forest. The percentiles of each scalar were extracted and the constraining models were fitted. Shown is impact of this percentile fitting on a) GPP; b) bias (modelled GPP subtracted from field measured); c) rootmean-square-error (RMSE); d) coefficient of determination (R 2 ); and e) slope and f) intercept from an ordinary least square linear regression fitted between modeled and field measured GPP. Note that scales differ for the different biomes ( Fig. S17-S26).

Figure S21
Sensitivity of LRF-GPP to variability in air temperature scalar for savanna/shrublands. The percentiles of each scalar were extracted and the constraining models were fitted. Shown is impact of this percentile fitting on a) GPP; b) bias (modelled GPP subtracted from field measured); c) root-mean-square-error (RMSE); d) coefficient of determination (R 2 ); and e) slope and f) intercept from an ordinary least square linear regression fitted between modeled and field measured GPP. Note that scales differ for the different biomes ( Fig. S17-S26).

Figure S22
Sensitivity of LRF-GPP to variability in the air temperature (Tair) and vapor pressure deficit (VPD) scalars for grasslands. The percentiles of each scalar were extracted and the constraining models were fitted. Shown is impact of this percentile fitting on a) GPP; b) bias (modelled GPP subtracted from field measured); c) root-mean-square-error (RMSE); d) coefficient of determination (R 2 ); and e) slope and f) intercept from an ordinary least square linear regression fitted between modeled and field measured GPP. Note that scales differ for the different biomes ( Fig. S17-S26).

Figure S23
Sensitivity of LRF-GPP to variability in the air temperature scalar for croplands. The percentiles of each scalar were extracted and the constraining models were fitted. Shown is impact of this percentile fitting on a) GPP; b) bias (modelled GPP subtracted from field measured); c) rootmean-square-error (RMSE); d) coefficient of determination (R 2 ); and e) slope and f) intercept from an ordinary least square linear regression fitted between modeled and field measured GPP. Note that scales differ for the different biomes ( Fig. S17-S26). Table S3. Table of (4 S7 Evaluation of global-scale Earth observation GPP products S7.1 Evaluation against FLUXNET2015 GPP  Figure 25 Evaluation of seasonal dynamics in modelled GPP (black) against GPP from the FLUXNET2015 sites (red) for the different biomes; a) evergreen needle leaf forest; b) evergreen broadleaf forest; c) deciduous needleaf forest; d) deciduous broadleaf forest; e) mixed forest; f) savanna/shrublands; g) grasslands; and h) croplands. The uncertainty is ±one standard deviation.

Figure S26
Evaluation of seasonal dynamics in FLUXCOM GPP (black) against GPP from the FLUXNET2015 sites (red) for the different biomes; a) evergreen needle leaf forest; b) evergreen broadleaf forest; c) deciduous needleaf forest; d) deciduous broadleaf forest; e) mixed forest; f) savanna/shrublands; g) grasslands; and h) croplands. The uncertainty is ±one standard deviation.

Figure S27
Evaluation of seasonal dynamics in Kolby Smith GPP (black) against GPP from the FLUXNET2015 sites (red) for the different biomes; a) evergreen needle leaf forest; b) evergreen broadleaf forest; c) deciduous needleaf forest; d) deciduous broadleaf forest; e) mixed forest; f) savanna/shrublands; g) grasslands; and h) croplands. The uncertainty is ±one standard deviation.

Figure S28
Evaluation of seasonal dynamics in p-model GPP (black) against GPP from the FLUXNET2015 sites (red) for the different biomes; a) evergreen needle leaf forest; b) evergreen broadleaf forest; c) deciduous needleaf forest; d) deciduous broadleaf forest; e) mixed forest; f) savanna/shrublands; g) grasslands; and h) croplands. The uncertainty is ±one standard deviation.

Figure S29
Evaluation of seasonal dynamics in MOD17 GPP (black) against GPP from the FLUXNET2015 sites (red) for the different biomes; a) evergreen needle leaf forest; b) evergreen broadleaf forest; c) deciduous needleaf forest; d) deciduous broadleaf forest; e) mixed forest; f) savanna/shrublands; g) grasslands; and h) croplands. The uncertainty is ±one standard deviation.

Figure S30
Evaluation of seasonal dynamics in SMAP GPP (black) against GPP from the FLUXNET2015 sites (red) for the different biomes; a) evergreen needle leaf forest; b) evergreen broadleaf forest; c) deciduous needleaf forest; d) deciduous broadleaf forest; e) mixed forest; f) savanna/shrublands; g) grasslands; and h) croplands. The uncertainty is ±one standard deviation.

Figure S31
Evaluation of seasonal dynamics in VPM GPP (black) against GPP from the FLUXNET2015 sites (red) for the different biomes; a) evergreen needle leaf forest; b) evergreen broadleaf forest; c) deciduous needleaf forest; d) deciduous broadleaf forest; e) mixed forest; f) savanna/shrublands; g) grasslands; and h) croplands. The uncertainty is ±one standard deviation.
S7.2 Evaluation against SIF data extracted from the FLUXNET2015 sites Figure 32 Evaluation of seasonal dynamics in modelled GPP (black) against Solar Induced Fluorescence (SIF; red) for the different biomes; a) evergreen needle leaf forest; b) evergreen broadleaf forest; c) deciduous needleaf forest; d) deciduous broadleaf forest; e) mixed forest; f) savanna/shrublands; g) grasslands; and h) croplands. The uncertainty is ±one standard deviation. Figure S33Evaluation of seasonal dynamics in FLUXCOM GPP (black) against Solar Induced Fluorescence (SIF; red) for the different biomes; a) evergreen needle leaf forest; b) evergreen broadleaf forest; c) deciduous needleaf forest; d) deciduous broadleaf forest; e) mixed forest; f) savanna/shrublands; g) grasslands; and h) croplands. The uncertainty is ±one standard deviation.

Figure S34
Evaluation of seasonal dynamics in Kolby Smith GPP (black) against Solar Induced Fluorescence (SIF; red) for the different biomes; a) evergreen needle leaf forest; b) evergreen broadleaf forest; c) deciduous needleaf forest; d) deciduous broadleaf forest; e) mixed forest; f) savanna/shrublands; g) grasslands; and h) croplands. The uncertainty is ±one standard deviation.

Figure S35
Evaluation of seasonal dynamics in the p-model GPP (black) against Solar Induced Fluorescence (SIF; red) for the different biomes; a) evergreen needle leaf forest; b) evergreen broadleaf forest; c) deciduous needleaf forest; d) deciduous broadleaf forest; e) mixed forest; f) savanna/shrublands; g) grasslands; and h) croplands. The uncertainty is ±one standard deviation.

Figure S36
Evaluation of seasonal dynamics in MOD17 GPP (black) against Solar Induced Fluorescence (SIF; red) for the different biomes; a) evergreen needle leaf forest; b) evergreen broadleaf forest; c) deciduous needleaf forest; d) deciduous broadleaf forest; e) mixed forest; f) savanna/shrublands; g) grasslands; and h) croplands. The uncertainty is ±one standard deviation.

Figure S37
Evaluation of seasonal dynamics in SMAP GPP (black) against Solar Induced Fluorescence (SIF; red) for the different biomes; a) evergreen needle leaf forest; b) evergreen broadleaf forest; c) deciduous needleaf forest; d) deciduous broadleaf forest; e) mixed forest; f) savanna/shrublands; g) grasslands; and h) croplands. The uncertainty is ±one standard deviation.

Figure S38
Evaluation of seasonal dynamics in VPM GPP (black) against Solar Induced Fluorescence (SIF; red) for the different biomes; a) evergreen needle leaf forest; b) evergreen broadleaf forest; c) deciduous needleaf forest; d) deciduous broadleaf forest; e) mixed forest; f) savanna/shrublands; g) grasslands; and h) croplands. The uncertainty is ±one standard deviation.

S8 Impact of land cover change
An important factor affecting the spatiotemporal dynamics in the LRF modelled GPP budgets are land cover changes (Tagesson et al., 2020). In the current model set-up, land cover was set as static 2001 land cover from the MCD12C1 product. The main reason was that, as far as the authors are aware of, there is no dynamic land cover product covering the full study period. Other GPP products based on the GIMMS data also used static land cover in their GPP simulations (Jung et al., 2011;Kolby Smith et al., 2015;Stocker et al., 2019).
In an attempt to study the impact of land cover change, we made an analysis of change in land cover 2001-2015, i.e. the period over which we have a dynamic land cover product (MCD12C1). The biome with the largest negative land cover change between 2001 and 2015 was evergreen broadleaf forest, whereas the most positive changes were seen for cropland, grassland, and savanna/shrublands (Fig S31b). To show the regions most influenced by land cover change, we subtracted the biome coverage percentage 2015 from that of 2001 ( Fig S31c). The LRF-GPP simulations of the pixels with larger land cover changes (dark in Fig S39c) may be more uncertain than those with those with less land cover change (orange in Fig S39c).
We did a rough estimate of impact of land cover change on the global LRF-GPP budget by taking biome-wise average GPP based on the static land cover data set, and then applied these to the biomes of the dynamic land cover data (Fig S39d). The impact of land cover change on the average global GPP budget 2001-2015 was minor (annual GPP with static land cover: 123.95±1.33 Pg C (average±one standard deviation based on inter-annual variability); and with dynamic land cover: 123.55±1.29 Pg C). However, the impact of land cover change was larger on the estimated trends 2001-2015, (trend in GPP based on static land cover: 0.082±0.080 Pg C y -1 ; and based on dynamic land cover: 0.042±0.079 Pg C y -1 ). Still, both trends 2001-2015 were insignificantly different from 0, and there was no significant difference between the two estimates, indicating that the uncertainty by a lack of land cover change is within the inter-annual variability of the simulated GPP budgets.
On the other hand, it can clearly be seen that the GPP including a dynamic land cover is getting increasingly lower than static land cover GPP 2001-2015 (Fig S39d). Taking the difference between the two generates an estimate of impact of land cover change on global GPP 2001-2015. It is then clear that land cover change has had a significant impact on GPP 2001-2015 (-0.038±0.006 Pg C y -1 ; R 2 =0.74; p-value<0.001). For future studies on the impact of factors affecting the spatiotemporal GPP patterns, it is necessary to include a dynamic land cover data set in the LRF-GPP simulations.