Ecosystem carbon use efficiency at global scale from upscaling eddy-covariance data with machine learning and MODIS products

The carbon use efficiency (CUE) is an important ecological indicator that quantifies the capacity of terrestrial ecosystems to act as sinks for carbon transferred from the atmosphere. It is considered as one of the three key axes that capture most of the variability within ecosystem functions¹; the other two axes are the maximum ecosystem productivity and the ecosystem water-use efficiency. The knowledge of the three axes at global scale is crucial to understand and quantify the response of ecosystems to climatic and other environmental changes and to map the overall ecosystem functioning.

“Carbon use efficiency is a conceptually simple parameter. The reality of measuring CUE is, however, methodologically challenging.” This was stated by Bradford & Crowther² in 2013, who point out that CUE requires measurement of both carbon uptake and associated growth. The gross primary productivity (GPP) –which is not measured in situ directly but obtained from eddy covariance (EC) measurements– is frequently used in the literature to estimate the amount of carbon uptake: it quantifies the carbon fixed by plants through photosynthesis. CUE is usually defined for a plant community as the ratio of net primary productivity (NPP) to GPP, being NPP = GPP – R_a, where R_a refers to the autotrophic (plant) respiration. GPP represents the capacity of plants to transform carbon into new biomass. A higher value of CUE translates to greater growth per unit of carbon mass acquired. However, since carbon can be also stored in soils and thus released to atmosphere by both autotrophs and heterotrophs (microbial or heterotrophic respiration, R_h), the “ecosystem carbon use efficiency” is introduced^3,4,5 and defined as the ratio of the net ecosystem productivity (NEP = GPP – (R_a + R_h)) to GPP, CUE = NEP/GPP. Under non-steady-state conditions, NEP is retained in the ecosystem as phytomass and/or soil organic carbon⁵ and, thus, the ecosystem carbon use efficiency represents the potential carbon sink capacity of the ecosystem. According to this, ecosystem CUE values are lower than those of plant communities. From now on, we use CUE to refer to that obtained using the total ecosystem respiration (R_a + R_h) and CUE’ to refer that including only R_a (Fig. 1). The plant-centric view of CUE’ overlooks the critical role of R_h, which can represent a significant pathway for carbon return to the atmosphere. Although R_h is usually smaller than the other carbon fluxes, it can represent a significant percentage of total ecosystem respiration: 23%–32% for forests, 47%–57% for prairies and tundra, and 48%–54% for agricultural systems⁶; therefore CUE(:<)CUE’, and this must be considered when comparing carbon use efficiency data from different sources. An exclusive focus on CUE’ may therefore provide an incomplete picture of the ecosystem’s overall capacity to retain carbon. To address this conceptual gap, this study aims to provide a global quantification of ecosystem-level CUE. By defining CUE based on NEP, our approach explicitly incorporates both autotrophic and heterotrophic respiration. This holistic metric moves beyond the producer level to assess the carbon use efficiency of the entire ecosystem offering a more comprehensive assessment of how effectively terrestrial ecosystems assimilate and sequester atmospheric carbon.

Ecosystem-atmosphere carbon fluxes can show both spatial and temporal variability^7,8 due to changes in ecosystem type and composition, weather patterns, and phenology, among others⁹. Therefore, both spatial and temporal changes in CUE might be expected¹⁰. CUE shows less variability than the carbon fluxes used to compute it. While net and gross productivities and ecosystem respiration vary by two orders of magnitude across biomes¹¹, CUE values remain within a narrower range (if we exclude ecosystems with negative NEP). However, this does not imply a single, universal CUE value for any given ecosystem type. Jin et al.¹², using EC data, propose a global average CUE of 0.50 ± 0.13, with the greatest values corresponding with croplands and the lowest with mixed forests. Chen & Yu¹³ proposed an averaged value of CUE’ = 0.537 ± 0.114 across the whole China region, which showed large spatial variations associated to climate factors (mean annual temperature and) precipitation, and to ecosystem types. Particularly, for forest ecosystems –which are the main carbon sinks–, substantial variation in CUE’ has been reported in the literature¹⁴ because tree respiration is not a constant fraction of photosynthesis but it depends on forest age and type. The highest CUE’ values are found in young forests and in temperate deciduous forests.

Both the processes of photosynthesis and respiration are affected by environmental conditions such as elevated temperatures, water shortage and increased drought stress typical of global warming¹⁵. Therefore, a CUE temporal variability is expected¹². CUE is usually related to temperature⁵ –which is the main controlling factor for CUE variations¹⁶– and to precipitation and water availability¹³. For example, drought events that affected southeast Europe during the 2000–2014 period reduced the CUE by 10 to 20% and, as a result, the region shifted from a carbon sink to a carbon source¹⁷.

Despite all the variations in spatial and temporal scales, terrestrial carbon cycle models have frequently assumed fixed values of CUE’ or CUE, that is, fixed values for the ratio R_a/GPP or (R_a + R_h)/GPP, respectively. However, the literature has sufficiently shown the weakness of this assumption, highlighting that the knowledge of CUE variability would contribute to improve the carbon fluxes estimates and their prediction to understand ecosystem functioning^14,16. This highlights the need to map annual CUE at global scale (establishing differences between ecosystem types) and to analyze temporal and spatial variations to detect areas vulnerable to degradation under external disturbances (negative CUE trend), or even areas that can shift from being a carbon sink to a carbon source.

In situ data provide insights into CUE patterns; however, their scarcity, driven by vegetation heterogeneity and limited sampling, constrains accurate CUE mapping at the global scale. For example, Desai et al.¹⁸ evaluate the variation in NEP in the upper Great Lakes region of the United States of America (USA) and show the challenges of upscaling the fluxes across space and the limits of using data from very tall EC towers as proxies for regional fluxes. These in situ carbon flux measurements are only representative of an area determined by the tower’s footprint area (local scale): it can range from a few hundred meters to a few kilometers depending on tower height, canopy characteristics (heterogeneity), and wind velocity¹⁹. The upscaling of these punctual measurements to global scale can be achieved by exploiting machine learning and deep learning approaches such as random forests²⁰, neural networks²¹ and Gaussian processes (GPs)^7,22, among others.

In this framework, the main goal of this manuscript is to provide a tool to obtain ecosystem CUE at the global scale through the upscaling of in situ EC CUE exploiting Gaussian Processes Regression (GPR) and remote sensing (RS) observations from the MODerate resolution Imaging Spectroradiometer (MODIS). To our knowledge, no 1-km global-scale and long-term multitemporal validated carbon use efficiency maps adopting the ecosystem carbon use efficiency definition are currently available. The application of this definition in the existing literature has been predominantly confined to site-specific or regional analyses, thereby lacking the spatial coverage and temporal consistency required for a global-scale assessment. Our approach computes global, long-term CUE estimates, accounting for the contribution of total ecosystem respiration, thereby filling an important gap in ecosystem carbon research. The results are evaluated as a function of climatic zones and biome types, considering –in addition to literature data– the CUE’ that is obtained from the MOD17A3HGF product. To highlight the usefulness of the proposed CUE, changes in annual CUE over the 2001–2023 period are analyzed to determine CUE potential for detecting particularly sensitive or vulnerable areas at the global scale.

This manuscript is organized into 4 sections. Both in situ and RS data are described in the section “Materials and methods”, which also includes the description of the upscaling approach to obtain global CUE. The section “Results and discussion” contains the results and discussion concerning CUE at global scale including the spatial variability in terms of climate and ecosystem types, and the annual CUE map obtained throughout time series of 22 years to detect trends. Finally, the section “Conclusions” summarizes the main conclusions. Further details of the results are included in Supplementary Information online.

In situ data

The CUE upscaling process implies the training and execution of a machine learning model, which is fed by MODIS observations that are used as model inputs whereas CUE is the model output. In this study, in situ CUE was computed using concomitant GPP and NEP from the FLUXNET network (https://fluxnet.org). The in situ GPP and NEP measurements are based on the use of the EC technique, which is the standard method to measure directly (in situ) trace gas fluxes between ecosystems and atmosphere¹⁹. In particular, the most recent FLUXNET data product is used, the FLUXNET2015 Dataset (https://fluxnet.org/data/fluxnet2015-dataset/). Detailed descriptions of this dataset, and the method to partition NEP into GPP and ecosystem respiration are found in Pastorello et al.²³ and Reichstein et al.²⁴, respectively. The FLUXNET network allowed us to gather daily data from 211 EC flux towers distributed around the world. The towers are located in almost all the representative ecosystems (mainly forests and grasslands) from 77º N to 57º S¹⁹. Some of these sites have been collecting data for several decades, allowing the study of ecosystems over time. Supplementary Fig. S1 available online shows the location, dominant ecosystem, and climate type over every flux tower. Finally, daily data were temporally aggregated to match the temporal resolution (8-day) of the MODIS products. It is worth mentioning that the temporal coverage of EC flux towers in the Fluxnet2015 dataset is not uniform: while some sites only provide between one and four years of observations, others cover more than a decade. A quality control was applied to retain only the highest quality samples, selecting those that simultaneously met the highest quality EC tower data and MODIS observations. This filtering process yielded a total of 2912 samples available for model training and evaluation. The temporal distribution of the samples spans different years and months, with a lower proportion during the December-March period (see Supplementary Fig. S2 online). Although this number may seem limited for a global application, it is comparable to datasets used in other studies that have successfully implemented GPR models for global mapping²⁵.

Remote sensing observations

The inputs selected in this study are RS observations of the MODIS sensor, which are related with environmental and physiological variables, and biosphere-atmosphere interactions. In particular, eight MODIS-based predictors from five different products are selected (see Table 1): the MCD18C2 product²⁶ that provides Photosynthetically Active Radiation (PAR); the MOD11A2 product²⁷, which provides day- and night-time land surface temperature (LST_D, LST_N); the MOD16A2 product²⁸, which derives Evapotranspiration (ET) and Potential Evapotranspiration (PET), and allows to compute a water stress factor (C_ws) as C_ws = ET/PET; the MCD15A2H product²⁹ that provides the Leaf Area Index (LAI); and the kernel version of the normalized difference vegetation index (kNDVI)³⁰ computed from the MCD43A4 product³¹. All these products were downloaded at 1-km spatial resolution using Google Earth Engine (GEE) over the locations of the EC towers from 2001 to 2023. As previously mentioned, every product’s quality flag was used to filter only best quality pixels. Eventually, an annual CUE’, computed as CUE’ = NPP/GPP from the MOD17A3HGF product³² has been also downloaded from GEE to be compared with the CUE retrieved by the GPR.

Gaussian processes regression

GPs are Bayesian tools for discriminative machine learning³³. When applied to regression, they are known as Gaussian process regression, which provides a probabilistic approach to nonparametric kernel-based regression methods. GPR assumes a Gaussian process prior governs the set of possible unobserved (latent) functions. Both the observations and the likelihood yield posterior probabilistic estimates. The standard regression model expresses the output as the sum of an unknown latent function f(x) of the inputs and independent Gaussian noise of the form (:mathcal{N}(0,{sigma:}_{n}^{2})). GPR adopts a Bayesian, non-parametric framework by placing a zero-mean Gaussian process prior (:GPleft(0,{k}_{{uptheta:}}right)) on the latent function. Here, (:{k}_{{uptheta:}}) is a covariance function parameterized by hyperparameters (:varvec{theta:}). Given this prior, samples of f(x) evaluated at input locations (:mathbf{X}=left{{mathbf{x}}^{left(1right)},{mathbf{x}}^{left(2right)},dots:,{mathbf{x}}^{left(Nright)}right}), where N is the number of samples, each one having P features, follow a joint multivariate Gaussian distribution with zero mean and covariance matrix K with (:{left[mathbf{K}right]}_{ij}=:{k}_{{uptheta:}}:({mathbf{x}}^{left(iright)},{mathbf{x}}^{left(jright)}).).

For a test input (:{mathbf{x}}_{mathbf{*}}) associated to a scalar output (:{text{y}}_{text{*}}), the GPR induces a prior distribution between the observations y and (:{mathbf{x}}_{mathbf{*}}) with the corresponding output (:{text{y}}_{mathbf{*}}). Given a dataset (:mathcal{D}equiv:{mathbf{X},y}), the posterior distribution over a, unknown output (:{text{y}}_{text{*}}) can be derived analytically as:

where (:{mathbf{k}}_{*} = [k({text{x}}_{*} ,{text{x}}^{{left( 1 right)}} ), ldots :,k({text{x}}_{*} ,{text{x}}^{{left( N right)}} )]^{{{ top }:}}) is an N-dimensional vector, (:{k}_{**}=kleft({text{x}}_{*}{,text{x}}_{*}right)), and (:{mathbf{I}}_{text{n}}) is the N×N identity matrix. Note that the GPR not only offers pointwise estimations, (:{mu:}_{GP*})but also confidence estimates (:{sigma:}_{GPR*}^{2}). The relation between the input and the output is established as:

where α_i is the weight assigned to each training predictor, α₀ is the bias in the regression function, and k_θ is the aforementioned covariance (kernel) function that evaluates the similarity between test and training data. In this paper, we selected the automatic relevance determination kernel:

where ν is a scaling factor, (:{sigma:}_{n}) accounts for the noise standard deviation, P is the number of predictors, and (:{sigma:}_{p}) is a dedicated parameter controlling the spread of the relations for each particular predictor p. The inverse of (:{sigma:}_{p}) represents the relevance of each predictor p. Model hyperparameters are collectively grouped in (:{varvec{theta:}=[nu:,sigma:}_{n,:}{sigma:}_{1},:dots:,:{sigma:}_{p}]) and model weights can be automatically optimized by maximizing the marginal likelihood in the training set³³:

Algorithm performance and CUE estimates

The GPR model is evaluated with 2912 data pairs (8 predictors–CUE) obtained matching the in situ CUE data with the RS products. As the predictors exhibit a substantial range of variation among them, the inputs are scaled to 0–1 range to mitigate the risk of suboptimal performance. The calibration of the model is carried out by means of varying the number of training–validation data. Eleven cases of training–validation percentages are considered: 1% – 99%, 5% – 95%, 10% – 90%, 20% – 80%, 30% – 70%, 40% – 60%, 50% – 50%, 60% – 40%, 70% – 30%, 80% – 20%, and 90% – 10%. For every case, a hundred random initializations are conducted, and both the root mean square error (RMSE) and the coefficient of determination (R²) are computed to assess the GPR performance. Validation statistics were averaged across iterations, ensuring that all sites and years were represented in both training and testing datasets, thereby accounting for spatial and temporal variability. Subsequently, in view of the obtained results (see Results section), model evaluation was conducted using the 70%–30% split, enabling both biome-specific validation and assessment of predictor relevance according to the (:{sigma:}_{p}) values provided by the GPR. Finally, the GPR model (trained with all samples) is executed to obtain 1–km multitemporal CUE estimates at global scale from 2001 to 2023. This allows to report average CUE values and trends over different climate classes and biomes.

Model performance over tower data

Figure 2 shows the assessment of the GPR model. The effect of model training/testing size reveals that for small training datasets, the GPR model does not perform optimally (Fig. 2(a)). Increasing the ratio of training/testing samples results in higher R² values and lower RMSE. The values are quite constant when datasets containing ≥ 70% of training samples are used for model training. Therefore, the 70% – 30% dataset is used to assess the accuracy of the CUE predictions provided by the GPR model. Figure 2(b) shows the assessment of the CUE predictions over the validation set composed of in situ EC data never used during model training. The accuracy metrics demonstrate strong agreement between the estimates and in situ measurements. Notably, the R² reaches 0.84 and the RMSE remains low (0.10), indicating strong predictive performance. The GPR CUE predictions reveal low mean error (ME = 0.01) and mean absolute error (MAE = 0.06). The GPR model provides the relevance of every predictor as the inverse of (:{sigma:}_{p}). Figure 2(c) shows this relevance normalized with the value of the maximum relevance predictor. The LAI is identified as the most relevant predictor for CUE, followed by the ET and LST_D. This is in accordance with other studies such as that by Liu et al.³. Through a hierarchical partitioning analysis between the explanatory variables and CUE, these authors found that LAI exerted greater impact than other factors related to climate (as temperature and precipitation). The accuracy metrics are disaggregated per biome type as shown in Table 2. The best results in all metrics are obtained over wetlands and evergreen needleleaf forests but must be biased by the low number of samples compared with the rest of biomes. The poorest agreement between estimates and in situ data is obtained over savannas.

Global CUE

The execution of the GPR model provided estimates at a 1-km spatial resolution and 8-day temporal frequency, which were subsequently aggregated into annual means for the 2001–2023 period. In this section, the mean global CUE over this period is presented.

The global mean CUE computed by the GPR model (0.43 ± 0.08) is, as expected, lower than the CUE’ obtained from MODIS (0.49 ± 0.06) during the same period. Similar underestimation behavior is found when compared with other studies that compute CUE’. For instance, Zhang et al.³⁵ report a CUE’ value of 0.52 from MODIS data in the 2000–2003 period. He et al.³⁶ employed five process-based models to estimate a global average CUE’ of 0.45 ± 0.05, slightly higher than the global mean obtained in the present study using GPR. In addition, the same study also reported a CUE’ value of 0.48 ± 0.05 from MODIS data in the 2000–2012 period, which is very similar to the one obtained for MODIS in the present study for the 2001–2023 period. Tang et al.³⁷ and Jin et al.¹² found a CUE’ value of 0.488 ± 0.136 and 0.50 ± 0.13, respectively.

Figure 3 exhibits the global spatial distribution of the mean CUE, associated GPR pixel-wise uncertainty, and the spatial difference concerning CUE’ computed from the MOD17A3HGF product. A considerable global CUE variation, roughly ranging from 0.3 to 0.9 is found (Fig. 3, top). The sigma value provided by the GPR (Fig. 3, middle) ranges from 0.01 to 0.05. The spatial variation of CUE reported in Fig. 3 (top) highlights a clear latitudinal path, suggesting a dependence on climate variability. CUE increases with latitude: from minima around tropical zones (including tropical rainforests and tropical savannas), as also found by Jin et al.¹² and Gang et al.³⁸, to reaching maximum values in high northern latitudes (e. g., subarctic climate areas in Canada and Siberia). Street et al.³⁹ analyzed the European subarctic region and showed the influence of mosses –that have CUE values greater than those observed for vascular plants– to increase the ecosystem CUE. The highest northern latitudes show CUE lower than the subarctic areas. Consistent with other global upscaling approaches, the sparse spatial distribution of flux tower data is a known source of uncertainty in global upscaling. The GPR framework directly addresses this limitation by providing explicit uncertainty, thereby identifying regions where model estimates are likely to be less reliable. The spatial pattern of the uncertainty is quite constant along the latitudes as shown in Fig. 3 (middle), with higher uncertainties over very high latitudes and some zones over the tropics.

Figure 3 (bottom) shows the spatial difference ((:{varDelta:}_{text{C}text{U}text{E}})) between mean carbon use efficiencies obtained by the GPR model and MODIS. For the sake of clarity, from now on, CUE_GPR and CUE’_MODIS will refer to CUE and CUE’ values obtained by the GPR model and MODIS product (NPP and GPP from the MOD17A3HGF product), respectively. In general, a systematic CUE_GPR < CUE’_MODIS behavior is found, as expected by the different CUE definitions. However, there exist northern zones where very similar values are found. This can be partly due to very low soil heterotrophic respiration reported by other studies^40,41 over northern zones, which reduces the differences between CUEs. In addition, our findings suggest CUE_GPR ≥ CUE’_MODIS and high uncertainty over these northern zones (see Fig. 3 (middle)). Other studies⁴² also reported largest uncertainty over high latitudes.

The upscaling process is built upon the FLUXNET network, which has denser coverage in North America and Europe. To address this well-known limitation, the GPR model operates in the environmental feature space rather than geographic space. This approach leverages globally predictor variables to apply learned eco-physiological relationships by the GPR model. Despite the robustness of this methodology, the uncertainties provided by the GPR model could be larger in regions and biomes that are sparsely sampled. Expanding the in-situ monitoring network in these areas remains a critical priority for the global carbon cycle community. In addition, it is crucial to note that the CUE_GPR framework is not intended as a numerical ‘correction’ for previous CUE’_MODIS or similar estimates, but rather as a complementary metric that answers a different and broader ecological question. While CUE’ accurately reflects efficiency at the producer level, CUE evaluates the net outcome of carbon assimilation and total ecosystem respiration. The global patterns of CUE presented here therefore offer a new benchmark for the scientific community, particularly for evaluating and constraining the performance of Earth System Models, which must accurately simulate the complete ecosystem carbon balance. Although the GPR framework does not explicitly propagate uncertainties from flux tower observations and satellite inputs to the global scale —an interesting aspect that remains for future research— the dataset presented here contributes to improving our understanding of the terrestrial carbon sink. Both CUE_GPR and CUE’_MODIS are further analyzed in terms of the climatic areas according to the Köppen–Geiger climate classification⁴³, and as a function of the biome type.

Mean CUE for the different climatic classes

Given the latitudinal gradient observed in CUE values, the first analysis is in terms of the Köppen-Geiger climate classification, which was updated by Beck et al.⁴³ using high-resolution, observation-based climatologies. Supplementary Table S1 available online summarizes the findings (mean annual values and standard deviation) over the five major climatic classes: tropical, arid (excluding main deserts, masked out in Fig. 3), temperate, cold and polar. As shown CUE_GPR < CUE’_MODIS for all the climatic classes, except for cold areas, where it shows a slightly higher value. In addition, Fig. 4 shows the distributions (represented using boxplots) of annual CUE_GPR and CUE’_MODIS. The edges of each box represent the first quartile (25th percentile) and the third quartile (75th percentile), while the whiskers extend to the 5th and 95th percentiles. The black line within the boxes stands for the median CUE value. Note that the median CUE_GPR is always lower than the median CUE’_MODIS, even in cold areas. Complementarily, results for the sub-classes according to the Köppen-Geiger climate classification are provided on Supplementary Figs. S3–S7 online. These kinds of comparisons between different carbon use efficiency definitions are carried out to illustrate the ecological significance of incorporating R_h. As dictated by their theoretical relationship (CUE = CUE’ – R_h/GPP), CUE must be less than or equal to CUE’. The reported results, therefore, do not validate the GPR model’s predictive power, but rather provide a quantitative assessment of the R_h/GPP ratio’s impact on carbon retention efficiency across different climate zones. The differences reveal how ecosystems belonging to different climate types vary in their capacity to retain assimilated carbon after accounting for all respiratory losses.

As mentioned above, CUE values increase from tropical to cold climates and decrease again towards the polar zones. Concerning climatic differences during the 2001–2023 period, Supplementary Table S2 available online reports the mean CUE_GPR and CUE’_MODIS for every of the 30 sub-climatic classes. The values obtained by the GPR model are systematically lower than those reported by MODIS. However, over the temperate – dry winter and hot summer (Cwa), and cold – dry winter and hot summer (Dwa) climatic zones the mean values of both approaches are coincident (0.50 and 0.49, respectively). The highest mean CUE_GPR values are found on temperate zones dominated by dry winter and cold summer (Cwc) (0.57 ± 0.06), and subarctic zones such as cold – dry and cold summer (Dsc) (0.57 ± 0.11), cold – dry summer and very cold winter (Dsd) (0.57 ± 0.08), and cold – dry and very cold winter (Dwd) (0.57 ± 0.04). In contrast, the lowest CUE_GPR values are found over the tropical – rainforest climate class (Af) (0.35 ± 0.05). This similar behavior is shown by other studies¹². Over the tundra polar climatic class, CUE_GPR = 0.49 ± 0.03, which is very similar to 0.5 (computed as NEP/GPP) found by Reichle⁴⁴.

CUE mean values for the different biomes

For a given climatic zone, differences in CUE demonstrate that vegetation type has an obvious impact on CUE¹³. In fact, maintenance metabolism (R_a) costs are related to the biome type and to its adaptation to the environmental conditions. The biome type also determines the percentage of total ecosystem respiration that the R_h represents (lower for forests than for prairies and crops)⁶. Carbon fluxes, and therefore CUE values, are related to vegetation type. In this subsection, CUE values are analyzed as a function of the biome type (Fig. 5; Table 3). To account for the impact of land cover dynamics on CUE, we employed a year-specific methodology using the MODIS Land Cover Type (MCD12Q1) Version 6.1 product. For each year within our study period, we utilized the corresponding annual MCD12Q1 land cover map for that same year. For a given year, we calculated the CUE for every pixel. Then, using the land cover map for that year, we spatially aggregated these pixel-level CUE values based on the International Geosphere-Biosphere Programme (IGBP) classification scheme. This produced an annual time series of mean CUE for each IGBP biome. The overall aggregated results presented in this paper were then computed by combining these annual means over the entire 2001–2023 period. The analysis of temporal trends was conducted using this year-by-year time series of CUE for each biome. This dynamic approach ensures that our findings incorporate the influence of inter-annual land cover changes, such as those caused by deforestation, fires, and land use conversion, on ecosystem carbon dynamics. As mentioned in the introductory section, low CUE values imply that a little amount of carbon is converted to biomass and biological products, i.e., less carbon is retained in the organism and more is released⁴. It should be noted that, while climate drives the vegetation activity, a given vegetation type or biome can be found in areas showing different climatic conditions. Similar to the climatic zones results, median CUE_GPR < median CUE’_MODIS in all biomes (see Fig. 5). In OSH (open shrublands) and WET (wetlands), the whiskers of the boxplots suggest that, although the median CUE_GPR is lower than CUE’_MODIS, the mean value of the GPR could be equal or higher than that from MODIS, which is confirmed by values in Table 3.

In relation to the CUE dependence on biome type, Jin et al.¹² and Tang et al.³⁷ found mean CUE’ values for CRO equal to 0.58 ± 0.12 and 0.566 ± 0.145, respectively. Luo et al.⁴² reported CUE’ values of 0.50 ± 0.09, 0.46 ± 0.10 and 0.39 ± 0.10 for CRO, DBF and ENF, respectively, and identical values for EBF and SAV (0.32 ± 0.12). Our results suggest higher values for ENF, EBF and SAV (CUE_GPR = 0.47 ± 0.09, CUE_GPR = 0.37 ± 0.04 and CUE_GPR = 0.41 ± 0.08, respectively), but lower for CRO and DBF (CUE_GPR = 0.42 ± 0.08 and CUE_GPR = 0.43 ± 0.09). Tang et al.³⁷ reported a CUE’ value of 0.464 ± 0.127 for forests, and Jin et al.¹² 0.44 ± 0.13 for MF, which are lower than the CUE reported in this study for MF (CUE_GPR = 0.52 ± 0.09). Regarding GRA, CUE_GPR = 0.44 ± 0.07, which is slightly lower than the CUE’ obtained by Tang et al.³⁷ (0.457 ± 0.109).

Our results show that, for CSH (closed shrublands), the mean CUE_GPR = 0.42 ± 0.06 is clearly lower than that for OSH (open shrublands) (0.53 ± 0.12), indicating that the density or “openness” of the shrub canopy influences its CUE. For OSH, CUE_GPR = 0.53 ± 0.12 is slightly higher than the mean CUE’ (:(sim)0.51) reported by Jin et al.¹² and the one obtained in the present study with MODIS (CUE’_MODIS = 0.48 ± 0.14). Table 3 shows that CUE_GPR < CUE’_MODIS, except for OSH, which is consistent with the spatial pattern found in CUE (previous section). In the case of savanna biomes, CUE_GPR = 0.41 ± 0.08 and CUE_GPR = 0.45 ± 0.11 for SAV (savannas) and for WSAV (woody savannas) are obtained, respectively. This suggests that savanna CUE may be lower due to the prevalence of C4 photosynthesis in grasslands, which implies higher energy and carbon costs compared to C3 photosynthesis⁴⁵. Woody components in savannas tend to use the C3 pathway, which could contribute to higher CUE. Regarding WET (permanent wetlands), CUE_GPR = 0.50 ± 0.10 is reported, which is lower than the mean CUE’ value found by Tang et al.³⁷ (0.607 ± 0.133). However, our result is consistent with the CUE value of 0.504 over aquatic ecosystems reported by Reichle⁴⁴.

Trend

The computation of annual multitemporal estimates from 2001 to 2023 allows the provision of CUE trends at a global scale. All the trends (S) in this section have been quantified using the Theil-Sen slope estimator^46,47,48. Note that all slopes are significant according to Hamed and Rao modified Mann-Kendall test^49,50. S_GPR and S_MODIS will refer to the trend slopes obtained from CUE_GPR and CUE’_MODIS, respectively.

Figure 6 shows the pixel-wise trend (slope) of the annual CUE_GPR in the 2001–2023 period at global scale. The mean slope is (:{S}_{text{GPR}}=left(-1.2pm:0.3right)times:{10}^{-2}:{text{d}text{e}text{c}text{a}text{d}text{e}}^{-1}). This reveals a decreasing global carbon use efficiency during the last two decades, which coincides with the decline rate values expected at the end of the 21 st century as estimated by Chen et al.¹⁶, mainly attributed to the temperature increase. Zhang et al.⁵¹ also revealed a negative trend for global CUE from 2000 to 2009. The spatial pattern of the trend obtained from 2001 to 2023 is in accordance with Gang et al.³⁸, but disagrees from Zhang et al.⁵¹ in some areas (mainly South Africa). The discrepancy can be caused partly by the different and shorter temporal period used in that study (2000 to 2009), which highlights the importance of continuously assessing CUE to detect changes in its trend. The spatial pattern shows that the areas exhibiting significant positive changes ((:{S}_{text{GPR}}>0,:)in red color) are found in humid zones (as those in South America and in central-South Africa) as well as in arid zones (Australia and western USA) according to Zomer et al.⁵². In any case, an accelerated aridification has been identified by Sardans et al.⁵³ in these zones, who enhanced central Africa as a new hotspot. Areas with negative trends ((:{S}_{text{GPR}}<0)) are also found in both humid and arid zones: the highest decreasing CUE is reported over the Amazon (humid), India, and northwestern Australia (more arid).

Some of the above spatial patterns of negative trends (in South America, central Africa, and Southeast Asia) coincide with the forest carbon losses (including aboveground and belowground biomass carbon loss and soil organic carbon loss) identified by Feng et al.⁵⁴ across the tropics. Most of them are attributed to agricultural practices and expansion (replacing the rainforests), confirming a dominant role of agriculture in long-term pan-tropical carbon reductions on formerly forested landscapes. Wang et al.⁵⁵ show a net increase of the above-ground biomass carbon stock in Africa, but the rainforests present a loss. These authors highlight the human-induced deforestation and water stress (especially the vapor pressure deficit) as the most important variables explaining the spatial and temporal above-ground biomass variations. The negative CUE trend obtained in most part of Africa (not only in tropical forest but also in the transition zone to savanna woodlands and even in savannas) suggests the necessity of considering not only GPP but also respiration fluxes to further characterize the ecosystem state.

In South America, a region with negative CUE trend is shown (Fig. 6), which mainly comprises tropical and subtropical forests (Amazonia), grasslands and savannas (for example, the north of the Brazilian Cerrado). The Brazilian Cerrado, which is considered the most biodiverse savanna on the planet, is affected by a degradation process due to both human activity (agricultural expansion and forest fires) and climate variability⁵⁶. Simulated data by Delgado et al.⁵⁷ show that tropical and subtropical forest, and savannas in South America will continue to show a decrease in vegetation activity due to the expected increase in air temperature as well as increasing fires during the dry season. Globally, both positive and negatives trends span different climatic zones and very different ecosystems, which is analyzed in the next sections.

CUE trend for the different climatic classes

An analysis per climate type (Fig. 7; Table 4) revealed decreasing trends over all major climatic classes except in arid zones where low positive trend is found. The negative trends are reported for both the GPR model and MODIS ((:{S}_{text{GPR}}<0,:{:S}_{text{MODIS}}<0)), however, over temperate zones (:left|{S}_{text{GPR}}right|>left|{S}_{text{MODIS}}right|).

Complementarily, online Supplementary Figs. S8 to S12 show the CUE trend reported for all 30 climatic classes. All three tropical sub-types follow a negative trend (see Supplementary Fig. S8 online), whereas in the case of the arid sub-types the arid – steppe cold (BSk) shows a negative trend contrary to the rest of arid zones (see Supplementary Fig. S9 online). In the case of temperate sub-types there is a mixture of behaviors in trends (see Supplementary Fig. S10 online) that makes the final trend of temperate zones to be slightly negative. The cold climatic zone is largely dominated by the sub-type cold – no dry season and cold summer (Dfc) (see Supplementary Fig. S1 online), which makes that the Dfc trend (see Supplementary Fig. S11 online) has a major impact on the overall cold trend. Finally, the polar CUE trend is mainly driven by the sub-type polar – tundra (ET) as shown in Supplementary Fig. S12 online.

CUE trend for the different biomes

Figure 8 shows the temporal evolution of CUE_GPR and CUE’_MODIS per biome type. In general, the CUE_GPR temporal evolution resulted in slightly noisier temporal estimates than CUE’_MODIS. This can be also observed in the slope’s error, being generally higher in the case of the GPR model (see Table 5). These interannual variabilities are usually induced by machine learning models in long-term trends of carbon uptake⁵⁸. Both CUE_GPR and CUE’_MODIS show positive trends ((:{S}_{text{GPR}}>0,:{:S}_{text{MODIS}}>0)) for SAV, WSAV, and CRO, and negative trends ((:{S}_{text{GPR}}<0,:{:S}_{text{MODIS}}<0)) for DBF, EBF, ENF, DNF, GRASS, OSH, and CSH. In the case of MF (mixed forests), (:{S}_{text{GPR}}<0), while MODIS shows a very slight positive trend ((:{:S}_{text{MODIS}}>0)) but two order of magnitude lower in absolute value ((:sim)10^–6 yr^–1 vs. 10^–4 yr^–1). For WET, the trends reported by the two models are opposite, but also low: (:{S}_{text{GPR}}sim)10^–5 yr^–1, and (:{:S}_{text{MODIS}}sim) − 10^–6 yr^–1). Both (:{S}_{text{GPR}})and (::{S}_{text{MODIS}}) show the highest negative values over evergreen forests (ENF and EBF). The highest positive trend in CUE_GPR is found over croplands (CRO), whereas the highest positive trend in CUE_MODIS is found over wooded savannas. There is a lack of studies reporting CUE trends over different biomes at global scale, and only few studies address this subject. For instance, Yang et al.⁵⁹ find positive CUE trends over SAV and WSAV, and a negative trend for non-woody grasslands at global scale during the period 2000–2013. Du et al.⁶⁰ report a negative CUE trend of (:sim) − 5 × 10⁻⁴ yr⁻¹ for GRASS in the 2001–2017 period over the Ningxia province (northwest China), which is similar to the trend reported in this study at global scale (see Fig. 8). Lei et al.⁶¹ evaluate the 2001–2023 period and show a negative CUE’ trend for forests, grasslands, and croplands, and a positive trend over shrublands in the Nanling Mountains (700 km × 400 km region in China).

Having efficient tools to characterize and quantify the functioning of ecosystems at a global scale is critical for climate change adaptation. The main indicator used by the plant community is based on autotrophic CUE’, a metric of plant-level efficiency that does not account for subsequent carbon losses due to heterotrophic respiration. This omission can lead to an overestimation of the efficiency with which entire ecosystems retain carbon. The present study bridges this gap by calculating and presenting the first global dataset of an ecosystem CUE, providing a robust measure of global terrestrial carbon sequestration efficiency. This study develops a potent, robust tool to map ecosystem carbon use efficiency at global scale, which allows to quantify CUE changes over time. The approach relies on an upscaling data-driven methodology combining in situ data with remote sensing observations into a GPR regression algorithm to retrieve global CUE. Results provide strong evidence of the algorithm’s high performance assessed over in situ data never used during model training.

The GPR model fed with 1-km multitemporal annual data from 2001 to 2023 allows to provide a global mean annual CUE map, whose overall spatial pattern is in accordance with the studies related to atmosphere-biosphere carbon exchanges reported by the scientific community. Results show that the mean CUE values increase from tropical to cold climates and decrease again towards the polar zones. The findings revealed that the mean spatial-temporal CUE estimates obtained by the GPR model are generally lower than the mean CUE’ of MODIS reported by the present paper. This can be mainly explained by different definitions of carbon use efficiency: CUE from GPR model is defined as NEP/GPP, whereas MODIS CUE’ accounts for NPP/GPP. This distinction is crucial when comparing values from different approaches that use varying definitions for carbon use efficiency. A preliminary analysis of estimated mean CUE per biome type supports the hypothesis that the proportion and type of woody vegetation in savannas significantly influences its CUE, with woody savannas potentially being more efficient due to different photosynthetic pathways or carbon allocation strategies compared to non-woody savannas. This highlights the importance of accurately classifying and distinguishing between subtypes of savanna, as well as different shrubs, in remote sensing studies to understand biome-specific carbon dynamics.

The GPR model reported a negative trend in the global mean CUE value, with a slope (decline rate) equal to (:left(-1.2pm:0.3right)times:{10}^{-3}:{text{y}text{r}}^{-1}), which is in accordance with future scenarios of global carbon use efficiency. The reported negative global CUE trend suggests a decreasing ecosystems’ capacity to sequester atmospheric CO₂, which highlights the importance of assessing the terrestrial carbon uptake in the future.