X-BASE: the first terrestrial carbon and water flux products from an extended data-driven scaling framework, FLUXCOM-X

Mapping in situ eddy covariance measurements of terrestrial land–atmosphere fluxes to the globe is a key method for diagnosing the Earth system from a data-driven perspective. We describe the first global products (called X-BASE) from a newly implemented upscaling framework, FLUXCOM-X, representing an advancement from the previous generation of FLUXCOM products in terms of flexibility and technical capabilities. The X-BASE products are comprised of estimates of CO₂ net ecosystem exchange (NEE), gross primary productivity (GPP), evapotranspiration (ET), and for the first time a novel, fully data-driven global transpiration product (ET_T), at high spatial (0.05°) and temporal (hourly) resolution. X-BASE estimates the global NEE at −5.75 ± 0.33 Pg C yr⁻¹ for the period 2001–2020, showing a much higher consistency with independent atmospheric carbon cycle constraints compared to the previous versions of FLUXCOM. The improvement of global NEE was likely only possible thanks to the international effort to increase the precision and consistency of eddy covariance collection and processing pipelines, as well as to the extension of the measurements to more site years resulting in a wider coverage of bioclimatic conditions. However, X-BASE global net ecosystem exchange shows a very low interannual variability, which is common to state-of-the-art data-driven flux products and remains a scientific challenge. With 125 ± 2.1 Pg C yr⁻¹ for the same period, X-BASE GPP is slightly higher than previous FLUXCOM estimates, mostly in temperate and boreal areas. X-BASE evapotranspiration amounts to 74.7×10³ ± 0.9×10³ km³ globally for the years 2001–2020 but exceeds precipitation in many dry areas, likely indicating overestimation in these regions. On average 57 % of evapotranspiration is estimated to be transpiration, in good agreement with isotope-based approaches, but higher than estimates from many land surface models. Despite considerable improvements to the previous upscaling products, many further opportunities for development exist. Pathways of exploration include methodological choices in the selection and processing of eddy covariance and satellite observations, their ingestion into the framework, and the configuration of machine learning methods. For this, the new FLUXCOM-X framework was specifically designed to have the necessary flexibility to experiment, diagnose, and converge to more accurate global flux estimates.

1 Introduction

Energy, water, and carbon exchange between terrestrial surfaces and the atmosphere are key components of the Earth system and impact ecosystems, ecosystem services, weather, climate, and water availability. The exchange (or flux) can be directly observed using eddy covariance (EC) measurement systems (Baldocchi, 2019), which are installed on towers overlooking the ecosystem of interest. The EC stations typically represent an area of a few hundred square meters to a square kilometer. One key advantage of the EC methodology is the ability to provide near-continuous measurements with some records now exceeding 20 years (Pastorello et al., 2020), allowing for examination of flux variations from the order of 30 min to decades. EC systems also provide a unique perspective on the magnitude, temporal variability, and environmental sensitivity of ecosystem CO₂ uptake, water use, and local climate regulation (Baldocchi, 2019; Musavi et al., 2017; Bao et al., 2022). However, while many of the most pressing scientific knowledge gaps surrounding the delicate land carbon balance and the water cycle require spatially and temporally resolved flux patterns at continental to global scales, EC observations are confined to individual locations in space and limited periods in time (Kumar et al., 2016; Papale et al., 2015). Methodologies to transcend the gap between local and global scales are needed to ultimately support societally relevant activities of building greenhouse gas monitoring systems, taking informed climate and land management actions, and verifying the effectiveness of mitigation strategies (Baldocchi and Penuelas, 2019; Bonan et al., 2011; Novick et al., 2022).

Coordinated and consolidated data collections from EC networks are invaluable for the mapping of in situ fluxes to regional and global scales. For example, EC measurements aid both the parameterization (Huang et al., 2021) and the validation (Turner et al., 2006; Heinsch et al., 2006) of mechanistic models of ecosystem productivity and land surface processes. The latter generate widely used reference data sets for terrestrial carbon cycle applications (Zhao and Running, 2010; Ukkola et al., 2022). A complementary approach to modeling terrestrial fluxes at continental and global scales is of empirical nature and links observations of explanatory variables at the EC stations, particularly meteorological and remote sensing data, to the EC fluxes via machine learning models. This upscaling concept does not prescribe any mechanistic formulations and assumes that the EC observations cover all complexities of ecosystem functioning. Based on a trained machine learning model and globally gridded input data of the explanatory variables, EC fluxes can be mapped to the global scale.

First implementations of this flux upscaling concept emerged in the early 2000s. They focused on net ecosystem exchange (NEE) of CO₂ and utilized the growing EC networks in Europe (Papale and Valentini, 2003) and North America (Xiao et al., 2008). The release of the FLUXNET La Thuile Synthesis Dataset of harmonized EC data in 2007, as well as methodological improvements in the training of the machine learning models (Jung et al., 2009), led to the first global products of terrestrial CO₂ and water fluxes at a monthly time step and in 0.5° grids in 2011 (Jung et al., 2011). While good agreement of flux estimates derived from complementary process-based models with the upscaled global gross photosynthetic CO₂ uptake (gross primary productivity, GPP) and energy fluxes demonstrated the potential of the approach, important inconsistencies remained, in particular regarding the globally integrated NEE and its year-to-year variability (Jung et al., 2020, 2019, 2011; Zscheischler et al., 2017). A critical component of data-driven methods, including the flux upscaling methodologies described here, is understanding the related uncertainties. The multistage process of going from the actual measured values of high-frequency (10–20 Hz) measurements of three-dimensional wind components, temperature, water vapor, and CO₂ concentration to the standard half-hourly, gap-filled, long-term data records from synthesis data sets such as FLUXNET has multiple associated potential biases and uncertainties. The well-documented issues of gap filling (Vekuri et al., 2023; Soloway et al., 2017; Papale et al., 2006), instrumentation errors (Zhang et al., 2023; Fratini et al., 2014; Vitale et al., 2020; Rannik et al., 2016), and energy balance non-closure (Mauder et al., 2020; Leuning et al., 2012; Stoy et al., 2013; McGloin et al., 2018; Franssen et al., 2010) will likely have an effect on the resulting gridded global data sets.

In an effort to better understand the uncertainties associated with mapping of EC fluxes to larger scales, the FLUXCOM intercomparison initiative built an ensemble of flux estimates as a type of factorial experiment (Tramontana et al., 2016; Jung et al., 2019, 2020). The ensemble consisted of multiple machine learning algorithms, meteorological forcing data, and combinations of predictor variables, resulting in 120 individual upscaled estimates per flux. These were summarized in two overall ensemble configurations, which differed in the set of predictors and spatial–temporal resolution. Apart from creating a large ensemble, the FLUXCOM evaluation included a consistent site-level cross-validation analysis as well as cross-consistency checks with terrestrial flux estimates from independent approaches, such as complementary modeling concepts or observational surrogates. From a methodological point of view, the key lessons learned from FLUXCOM were the following:

1. The overall approach seems to be primarily limited by the input information given to the machine learning algorithms rather than by the ability of the algorithm to extract the information.

2. The largest qualitative differences among flux products were related to the set of the predictor variables rather than to the choice of the machine learning method or meteorological forcing.

3. The cross-consistency checks with global independent data are essential for supplementing site-level cross-validation.

4. The largest qualitative discrepancy with independent data was a very high (strongly negative) tropical NEE that was shared among all ensemble members.

In addition to the systematic intercomparisons in FLUXCOM, the empirical upscaling concept has been implemented for a series of regional- and global-scale applications, each of them adopting disparate and individual methodological choices (e.g., Ichii et al., 2017; Yao et al., 2018; Joiner and Yoshida, 2020; Virkkala et al., 2021; Dannenberg et al., 2023; Burton et al., 2023; Zhu et al., 2024). These potentially important choices relate to data treatment (quality control, gap filling, processing pathways), ingestion (sampling, as well as matching EC and space-borne observations), and methodological configurations (machine learning methods and their training configuration, choice of predictor variables, resolution). The disparity between the setups and implementations of the empirical upscaling approach impedes any direct comparisons among products and valid conclusions as to where potential differences originate. Hence, flexibility to explore the large methodological space, as well as the ability to diagnose and evaluate global products in parallel to site-level cross-validation, is required to understand the importance of individual methodological choices and to make progress in empirical upscaling of EC fluxes.

We are developing a modeling framework that allows experimenting with many of these methodological choices. We coin this extended upscaling environment FLUXCOM-X. Compared to the FLUXCOM intercomparison initiative and the products produced therein, it follows the same overall principles. The key innovation of FLUXCOM-X is that it represents a flexibly adjustable upscaling environment to systematically explore different methodological avenues. Lessons learned from FLUXCOM and other upscaling exercises further imply striving for the enhancement of the information content of the training data with aspects related to coverage and quality of EC measurements as well as complementarity, completeness, and quality of predictor variables. The flexibility of the way FLUXCOM-X is implemented will considerably reduce the latency with which innovations in the related fields of machine learning and space-based Earth observations as well as novel EC data can find their way to empirical flux upscaling. This in turn allows faster progress towards more accurate and fit-for-purpose global biogenic flux estimates. Here, we introduce and evaluate the initial “basic” set of products from this flexible framework, which we refer to as FLUXCOM-X-BASE products (or X-BASE for short; see Appendix A for an overview on the naming conventions).

X-BASE products were comparably generated to the original FLUXCOM ensemble using qualitatively similar predictor variables, i.e., remotely sensed vegetation indices and land surface temperatures from the Moderate Resolution Imaging Spectroradiometer (MODIS) along with meteorological variables. In contrast to FLUXCOM, the remotely sensed and the meteorological predictors are combined in the X-BASE products, and the subdaily variability is resolved. We furthermore made efforts to provide more and improved information to the machine learning models by enhancing coverage and quality of the training data and by further developing the processing of satellite predictor variables (Walther and Besnard et al., 2022). In this article, we show results for X-BASE NEE, GPP, evapotranspiration (ET), and for the first time transpiration (ET_T), for the period 2001–2020 at 0.05° spatial and hourly temporal resolution. We focus here on the evaluation and cross-consistency checks of X-BASE with previous FLUXCOM products and independent data streams. Our specific objectives are the following:

1. to describe the production of X-BASE products

2. to evaluate the X-BASE setup using site-level cross-validation

3. to assess qualitative differences of global patterns compared to previous FLUXCOM products with reference to independent flux estimates where possible

4. to synthesize lessons learned from this basic exercise to guide future FLUXCOM-X developments.

X-BASE products are freely available and serve as a baseline for future FLUXCOM-X developments (see data availability statement). Any future product releases originating from FLUXCOM-X will follow the naming convention X-[specific name].

2 Data and methods

The following section gives an overview on the essential methodological implementations and data choices adopted in the generation of X-BASE products.

2.1 Eddy covariance data

Eddy covariance data consisted of 294 sites from around the world though skewed towards a higher representation of temperate forests from North America and Europe. All EC data were collected, processed, and analyzed for quality by the station teams, before being processed using state-of-the-art approaches in the ONEFlux data processing pipeline (Pastorello et al., 2020). ONEFlux provides consistent quality checks, gap filling, and carbon flux partitioning. The data included were collected between 2001–2020 and available with a CC BY 4.0 license. Based on this criterion, data for each site came from one of five different sources based on most recent availability: FLUXNET 2015 (Pastorello et al., 2020), ICOS Drought 2018 (Team and Centre, 2020), ICOS Warm Winter 2020 (Team and Centre, 2022), or the most recent AmeriFlux or ICOS release as of December 2022. Table 1 lists all sites included as well as the associated digital object identifier specific to the associated release.

(Garcia et al., 2016)(Kutzbach, 2021)(Posse et al., 2016)(Wohlfahrt et al., 2016)(Cleverly and Eamus, 2016b)(Beringer and Hutley, 2016c)(Meyer et al., 2016)(Pendall and Griebel, 2016)(Beringer and Hutley, 2016b)(Beringer and Hutley, 2016f)(Beringer and Hutley, 2016e)(Schroder et al., 2016)(Beringer and Hutley, 2016a)(Macfarlane et al., 2016)(Beringer and Hutley, 2016d)(Liddell, 2016)(Cleverly and Eamus, 2016a)(Woodgate et al., 2016)(Beringer et al., 2016b)(Beringer et al., 2016a)(Arndt et al., 2016)(Beringer and Walker, 2016)(Team and Centre, 2022)(Team and Centre, 2022)(RI, 2021)(Team and Centre, 2022)(Team and Centre, 2022)(Team and Centre, 2022)(Vourlitis et al., 2022)(Saleska, 2016)(Goulden, 2016d)(Staebler, 2022)(Knox, 2022)(Christen and Knox, 2022)(Wagner-Riddle, 2021)(McCaughey, 2016)(Black, 2021)(Amiro, 2016b)(Goulden, 2016a)(Goulden, 2016b)(Goulden, 2016c)(Goulden, 2016g)(Goulden, 2016e)(Goulden, 2016f)(Black, 2016b)(Black, 2016a)(Margolis, 2016)(Amiro, 2016c)(Amiro, 2016a)(Amiro, 2016d)(Arain, 2016b)(Arain, 2016a)(Arain, 2022b)(Arain, 2016c)(Arain, 2022a)(Nouvellon, 2016)(Team and Centre, 2022)(Team and Centre, 2022)(Team and Centre, 2022)(Team and Centre, 2022)(Team and Centre, 2022)(Ammann, 2016)(Team and Centre, 2022)(Zhang and Han, 2016)(Dong, 2016)(Shi et al., 2016)(Zhou and Yan, 2016)(Chen, 2016k)(Shao, 2016b)(Tang et al., 2016)(Wang and Fu, 2016)(Shao, 2016a)(Team and Centre, 2022)(Sigut et al., 2016)(Team and Centre, 2022)(Team and Centre, 2022)(Team and Centre, 2022)(Team and Centre, 2022)(Team and Centre, 2022)(Team and Centre, 2022)(Team and Centre, 2022)(Team and Centre, 2022)(Team and Centre, 2022)(Team and Centre, 2022)(Team and Centre, 2020)(Team and Centre, 2022)(Team and Centre, 2022)(Lindauer et al., 2016)(Knohl et al., 2016)(Team and Centre, 2022)(RI, 2022)(Team and Centre, 2022)(Team and Centre, 2022)(Schneider and Schmidt, 2016)(Klatt et al., 2016)(Bernhofer et al., 2016)(Team and Centre, 2022)(Sachs et al., 2016)(Pilegaard and Ibrom, 2016)(Olesen, 2016)(RI, 2022)(Team and Centre, 2022)(Team and Centre, 2022)(Team and Centre, 2022)(Poveda et al., 2016)(Team and Centre, 2022)(Team and Centre, 2022)(Team and Centre, 2022)(Team and Centre, 2022)(Reverter et al., 2016b)(Reverter et al., 2016a)(Team and Centre, 2022)(Lohila et al., 2016)(Team and Centre, 2022)(Team and Centre, 2022)(Aurela et al., 2016a)(Team and Centre, 2022)(Team and Centre, 2022)(Aurela et al., 2016b)(RI, 2022)(Team and Centre, 2022)(Team and Centre, 2022)(RI, 2022)(Team and Centre, 2022)(Team and Centre, 2022)(Team and Centre, 2022)(Team and Centre, 2022)(Berbigier and Loustau, 2016)(RI, 2022)(Team and Centre, 2022)(Ourcival, 2016)(RI, 2022)(Team and Centre, 2022)(Valentini et al., 2016b)(RI, 2022)(Hansen, 2016)(Lund et al., 2016b)(Lund et al., 2016a)(Team and Centre, 2022)(Team and Centre, 2022)(Team and Centre, 2022)(RI, 2022)(Sabbatini et al., 2016c)(Sabbatini et al., 2016a)(Sabbatini et al., 2016b)(Matteucci, 2016)(Team and Centre, 2022)(Valentini et al., 2016a)(Gruening et al., 2016b)(Cescatti et al., 2016)(Team and Centre, 2022)(RI, 2022)(Team and Centre, 2022)(Spano et al., 2016)(Manca and Goded, 2016)(Team and Centre, 2022)(Valentini et al., 2016c)(Papale et al., 2016)(Team and Centre, 2022)(Gruening et al., 2016a)(Team and Centre, 2022)(Kotani, 2016b)(Kotani, 2016a)(Yepez and Garatuza, 2021)(Kosugi and Takanashi, 2016)(Dolman et al., 2016a)(Team and Centre, 2020)(Wolf et al., 2016b)(Wolf et al., 2016a)(Griffis and Roman, 2021)(Merbold et al., 2016)(Dolman et al., 2016b)(Team and Centre, 2022)(Team and Centre, 2022)(Belelli et al., 2016)(Ardö et al., 2016)(Team and Centre, 2022)(Team and Centre, 2022)(Team and Centre, 2020)(Team and Centre, 2022)(Team and Centre, 2022)(Team and Centre, 2022)(Christensen, 2016)(Boike et al., 2016)(Tagesson et al., 2016)(Billesbach et al., 2022)(Billesbach et al., 2016b)(Billesbach et al., 2016a)(Biraud et al., 2022)(Torn, 2016b)(Torn, 2016a)(Zona and Oechel, 2016a)(Euskirchen, 2022b)(Euskirchen, 2022c)(Euskirchen, 2022d)(Euskirchen, 2022a)(Rey-Sanchez et al., 2022b)(Rey-Sanchez et al., 2022a)(Goldstein, 2016)(Huggins, 2021)(Huggins, 2022c)(Huggins, 2022a)(Huggins, 2022b)(Chen and Chu, 2016b)(Desai, 2022a)(Desai, 2022c)(Desai, 2022d)(Desai, 2022b)(Bowling, 2016)(Oikawa, 2021)(Massman, 2016)(Massman, 2022)(Meyers, 2016b)(Forsythe et al., 2021)(Goslee, 2022)(Munger, 2016)(Liu et al., 2022)(Hollinger, 2022)(Matamala, 2016)(Euskirchen et al., 2022a)(Euskirchen et al., 2022b)(Zona and Oechel, 2016b)(Vivoni and Perez-Ruiz, 2022)(Brunsell, 2022a)(Brunsell, 2022b)(Drake and Hinkle, 2016a)(Drake and Hinkle, 2016b)(Hinkle, 2022)(Meyers, 2016a)(Fares, 2016)(Desai, 2016c)(Novick and Phillips, 2022)(Wood and Gu, 2022)(Law, 2016c)(Law, 2022)(Law, 2016a)(Law, 2016e)(Law, 2016d)(Law, 2016b)(Litvak, 2021)(Sturtevant et al., 2016)(Torn and Dengel, 2021)(Blanken et al., 2022)(Suyker, 2022)(Suyker, 2016b)(Suyker, 2016a)(Silveira, 2021)(Bohrer, 2021)(Bohrer and Kerns, 2022)(Chen et al., 2016)(Desai, 2016d)(Kobayashi and Suzuki, 2016)(Flerchinger, 2022c)(Baker et al., 2022)(Baker and Griffis, 2022a)(Baker and Griffis, 2021)(Baker and Griffis, 2022b)(Flerchinger and Reba, 2022)(Flerchinger, 2022a)(Flerchinger, 2022b)(Kurc, 2022)(Scott, 2016a)(Scott, 2016b)(Shortt et al., 2022)(Kusak et al., 2022)(Ewers and Pendall, 2016)(Desai, 2016b)(Baldocchi and Ma, 2016)(Valach et al., 2021)(Sturtevant et al., 2022)(Chamberlain et al., 2022)(Sanchez et al., 2016)(Valach et al., 2022)(Baldocchi, 2016)(Bohrer, 2022)(Gough et al., 2016)(Gough et al., 2022)(Baldocchi et al., 2016)(Desai, 2016a)(Chen and Chu, 2016a)(Scott, 2016d)(Chen, 2016g)(Chen, 2016e)(Chen, 2016j)(Chen, 2016b)(Chen, 2016d)(Chen, 2016a)(Chen, 2016h)(Chen, 2016i)(Chen, 2016c)(Chen, 2016f)(Litvak, 2022)(Scott, 2016c)(Network, 2022)

Table 1 Citation data for the 294 sites used in the X-BASE products.

Meteorological data measured at each site consisted of incoming shortwave radiation, air temperature, and vapor pressure deficit, of which all data were gap-filled using the marginal distribution sampling method (Reichstein et al., 2005), as well as the computed potential shortwave incoming radiation (top-of-atmosphere theoretical maximum radiation) for every hour. Carbon dioxide flux data consisted of gap-filled net ecosystem exchange (NEE, variable USTAR threshold 50th percentile i.e., NEE_VUT_50) and the corresponding gross primary productivity (GPP, nighttime partitioning method after Reichstein et al., 2005). Water flux data consisted of evapotranspiration (ET, no energy balance correction), which was converted from the latent energy and transpiration estimates based on the transpiration estimation algorithm (TEA; Nelson et al., 2018; Nelson, 2021). The TEA estimates transpiration based on the relationship between GPP and ET under conditions where surface evaporation and soil evaporation are expected to be minimal, accounting for residual non-transpiration evaporation. Note that estimates of GPP and transpiration are not direct measurements but instead are based on statistical relationships localized to each station and time period and thus contain their own assumptions and uncertainties. All data were aggregated to a common hourly time resolution, an overview of which can be found in Table 2.

Table 2 Fluxes to be predicted and predictor variables used in X-BASE. The units of the fluxes correspond to the native hourly resolution. Upon temporal aggregation as in some analyses in the presented results, the units may change.

Data from the EC data set that ultimately were used for training the models varied 12¹⁴×10⁶ site hours depending on the target variable (i.e., GPP, NEE, ET, or ET_T). Training of the machine learning algorithms was only conducted on hours where all input variables passed quality control. The quality control procedure consisted of two levels, with the first being that each hour must have at least one value of good quality measured or gap-filled with confidence (i.e., at least one half hour was either 0 or 1 based on the ONEFlux _QC flags). Second, a set of consistency tests were performed on each used variable to check the consistency both among variables and across sites. As the consistency flags were based on daily aggregates of the meteorological and flux data, entire days were removed if the test indicated inconsistencies among related variables. The consistency flag also checked the relationship between variables across sites, ensuring that the relationships found across the data are coherent. A detailed explanation of these consistency flags can be found in Jung et al. (2024).

2.2 Global meteorology

For the generation of global flux maps we used hourly meteorological data from ERA5 global reanalysis products at 0.25° (Hersbach et al., 2020). Variables included air temperature at 2 m height, incoming shortwave radiation at the surface, and vapor pressure deficit (computed from relative humidity, air temperature, and surface pressure). Units were converted to correspond to the site-level measurements, which were used for training the machine learning model, and the data were regridded to a 0.05° resolution using bilinear interpolation for every hour.

2.3 Satellite Earth observation

The X-BASE products are based on measurements of the MODerate Imaging Spectroradiometer (MODIS) of surface reflectance and land surface temperature from collection 006 at daily resolution. Missing records were gap-filled consistently in both the average time series per EC station and in the global gridded data following the procedures of the FluxnetEO data version 2 (Walther and Besnard et al., 2022; Walther, 2023).

2.3.1 Spectral vegetation indices

At site level we used surface reflectance in the first seven MODIS spectral bands from the MCD43A4 c006 reflectance data set (500 m and daily, where each daily value is inverted from all valid observations within a 16 d window; Schaaf and Wang, 2015b). The spectral vegetation indices computed from the reflectance data were the enhanced vegetation index (EVI; Huete et al., 2002), the spectral reflectance of vegetation in the near-infrared (NIRv; Badgley et al., 2017), and the normalized difference water index (NDWI) with MODIS band 7 as reference (Gao, 1996). We followed the procedure of the FluxnetEO data set version 2 (Walther and Besnard et al., 2022) for data acquisition from Google Earth Engine for all pixels in a cutout of 4×4 km² around each EC station, as well as for quality checks in terms of cloud, snow, and land cover; index values outside the defined ranges; and outliers. An iterative approach then determined both the strictness of the inversion quality of the bidirectional reflectance distribution function (BRDF, based on the MCD43A2 data, Schaaf and Wang, 2015a) and the set of pixels in a cutout that shall represent a given EC station. Appendix B1 outlines all technical details of the dynamic procedure.

Global data of BRDF-corrected surface reflectance stem from the MCD43C4 c006 data (Schaaf and Wang, 2015b), available in a climate modeling grid of 0.05° with the same temporal sampling and subject to the same removal of snow and water pixels and outlier values as at site level. The BRDF quality control of the global data followed the same dynamic approach (see Appendix B1), which maximized data availability especially in tropical regions.

2.3.2 Land surface temperature

Satellite observations of land surface temperature (LST) were based on the MODIS c006 TERRA observations, which are available every day at 1 km resolution (Wan et al., 2015). We selected the 1 km² pixel containing a specific tower and treated the two MODIS LST data streams as independent predictor variables, which represent clear-sky LST at a specific time of the day (namely around 10:30 and 22:30 local time). Quality checks and gap filling followed the procedure described in FluxnetEO version 2 (Walther and Besnard et al., 2022).

For the global spatialization of the flux estimates, we relied on climate modeling grid LST from the MODIS TERRA data sets (Wan et al., 2015) and apply consistent quality control and imputation of missing values like at site level.

2.3.3 Land cover

Land cover information used the IGBP global vegetation classification. Site-level classification was as reported by the principal investigators. Global data were based on the yearly resolved MODIS MCD12C1 c006 product (Friedl and Sulla-Menashe, 2015). In order to ease the transition between site and global land cover classifications, an intermediate classification scheme was utilized which translated each classification into characteristics (e.g., trees, crops, needleleaf, deciduous) based on whether the classification has a specific feature (value = 1.0), might have a specific feature (value = 0.5), does not have a specific feature (value = 0.0), or is unknown (value = −1.0). A full description of this intermediate classification system can be found in Appendix B2.

2.4 Machine learning method

All X-BASE products are based on gradient-boosted regression trees using the XGBoost library (Chen and Guestrin, 2016). XGBoost is known as a robust algorithm that is able to handle a variety of variable types (numeric, Boolean, categorical). Training was conducted using a two-thirds training subsampling ratio and a 0.05 learning rate. Hyperparameters were chosen based on preliminary testing to predict eddy covariance data, and a full list of hyperparameters can be found in Appendix B3. Boosting was stopped when no model improvement (based on mean squared error of validation data) was observed for 10 consecutive rounds, and the best-performing model was stored to generate predictions. In all cases, the model reached the stopping criteria relatively quickly, with the final number of boosting rounds between 80–230, depending on the flux.

2.5 Cross-validation

All cross-validation was performed using a 10-fold, leave-site-fold-out scheme, where each fold was constructed by randomly assigning each site to a fold. Leave-site-fold-out is distinct from traditional leave-one-out cross-validation as each fold is comprised of independent sites rather than fully randomized individual points. For each round of cross-validation, eight folds were used for training, one for validation, and the remaining one as the test fold for which the actual predictions were made. The leave-site-fold-out scheme ensures that no data from the sites in the test fold were ever seen by the algorithm during training and in turn iterated such that each site was in the test set once. As eddy covariance sites are sometimes clustered in the same location (e.g., as different treatments) and can therefore be both physically closely located and not truly independent, sites are assigned to the same fold if they are less than 0.05° apart to reduce overfitting. We evaluate the accuracy of the cross-validation models by computing the Nash–Sutcliffe modeling efficiency (NSE, Nash and Sutcliffe, 1970), where a negative NSE indicates a model accuracy that is worse than a mean prediction, while a value close to 1 indicates high model accuracy. We compute the NSE for each site and for a range of temporal scales from hourly to interannual.

2.6 Upscaling

The final step to train a model to use in the final global prediction step was identical to the training in the cross-validation, with the exception that, because no test fold was required, we used 9 of the 10 folds for the training, and validation was done on the remaining fold. The final trained models (one trained model for each target flux) were then used to predict fluxes at the global scales using the associated globally gridded input variables that correspond to those used at site level, as outlined in Table 2.

2.7 Previous FLUXCOM and independent global flux estimates

We compare X-BASE with upscaling results from FLUXCOM (Jung et al., 2019, 2020). As mentioned earlier, FLUXCOM comprised an ensemble of upscaling experiments that differed in the choice of machine learning method and meteorological forcing data. The FLUXCOM ensemble was summarized in two groups of setups that shared the same predictor variables and spatiotemporal resolution: the “remote-sensing-only” setup (RS) mostly used spaceborne observations of MODIS as explanatory variables and produced flux estimates every 8 d at 0.083° resolution, while the “remote sensing plus meteorology setup” (RS + METEO) produced daily flux estimates at half-degree resolution from meteorological predictor variables and an average seasonal cycle of satellite observations (Tramontana et al., 2016; Jung et al., 2019, 2020). Comparisons to FLUXCOM RS + METEO data sets always refer to the ensemble over multiple machine learning methods for all realizations driven by the ERA5 meteorology (Hersbach et al., 2020). RS + METEO uses average seasonal cycles of MODIS c005 observations. For the FLUXCOM RS setup we use the ensemble over all machine learning methods. Please note that both the previous RS runs and the X-BASE runs presented here are driven by data from MODIS c006, but the processing has changed in some aspects such as quality control and gap filling. For clarification, an overview on FLUXCOM(-X) naming conventions is given in Appendix A.

For evaluating X-BASE NEE globally, in particular its seasonal cycle and for different regions, we used two different atmospheric inversion model products: the Orbiting Carbon Observatory-2 (OCO-2) v10 model intercomparison project (Byrne et al., 2023) and the CarboScope inversion (Rödenbeck et al., 2018) version s99oc_v2022 (Roedenbeck and Heimann, 2022). Estimates from the OCO-2 came from the LNLGIS experiment, which combines satellite-based column-averaged CO₂ (XCO₂) retrievals and in situ CO₂ measurements as observational constraints in the assimilation and consists of 13 different ensemble members covering the period 2015–2020 with a monthly frequency and 1° spatial resolution (https://gml.noaa.gov/ccgg/OCO2_v10mip/index.php, last access: 6 May 2022). The CarboScope product consisted of a single inversion output at the same spatial resolution as OCO-2 but with a longer temporal period from 2001–2020. In each case, as the inversion products estimate net biome exchange, we subtracted from the inversions data fire emissions as estimated by the Global Fire Emissions Database, Version 4.1 (Randerson et al., 2017).

We compared temporal patterns of X-BASE GPP with the patterns in global retrievals of sun-induced chlorophyll fluorescence (SIF) from the Sentinel-5P TROPOMI instrument (Köhler et al., 2018), which under most conditions approximates the variability in GPP. For the comparison we used estimates of daily mean SIF applying a correction factor to instantaneous observations (Zhang et al., 2018) and averaged both X-BASE GPP and TROPOMI SIF to a temporal resolution of 16 d and 0.5° spatial grids for the common period April 2018 to December 2020.

X-BASE ET and ET_T were cross-compared with transpiration estimates from the Global Land Evaporation Amsterdam Model (GLEAM) v3.6a (Martens et al., 2017; Miralles et al., 2011). GLEAM also utilizes satellite and reanalysis data sets but in a more physically constrained way, relying on semiempirical models such as the Priestley–Taylor (Priestley and Taylor, 1972) and Gash models (Gash, 1979). Further comparisons were made to precipitation data from GPCC (Schneider et al., 2022).

2.8 Units across scales and conversion

Given the broad range of spatial and temporal scales reported here, we utilized different units depending on each use case. All base measurements and predictions are at hourly timescale (in units of µmol CO₂ m⁻² s⁻¹ for carbon fluxes and mm h⁻¹ for water fluxes; as reported in Table 2). In the case of carbon fluxes, all timescales which are not subdaily are reported in grams of carbon per unit time and area. Globally integrated values are reported in petagrams of carbon or square kilometers of water and are aggregated using a common land area estimate per grid cell.

3 Results

3.1 Cross-validation and data space

One important innovation in FLUXCOM-X compared to the previous FLUXCOM ensemble was the training database, which was larger due to an increase in both number of sites and years. Furthermore, the EC methodology has changed considerably in many aspects ranging from collection and processing to quality filtering in the last 15 years. We show here one illustrative example of the changes in the environmental space that is represented in the training samples for daily NEE: between daily vapor pressure deficit (VPD) and daily incoming shortwave radiation, the distribution of training samples was considerably broader in X-BASE compared to the RS + METEO ensemble (Fig. 1). Furthermore, the number of unique sites contributing to a certain VPD-radiation bin has increased (Fig. C1); i.e., the number of ecosystems sampled in each climatic condition has also increased. The increases were seen particularly at the margins of the distribution, i.e., for days with high VPD along the full radiation spectrum, and vice versa for days with high-radiation conditions along the full VPD spectrum. Of particular note, Figs. 1 and C1 demonstrate improved sampling across more sites with conditions of high VPD and radiation, i.e., hot and dry conditions, for X-BASE compared to RS + METEO. Overall for NEE, the number of sampled site days increased over 3-fold (552 878 to 183 216 for X-BASE and RS + METEO, respectively). Note, however, that X-BASE is modeled at an hourly instead of daily resolution, and thus the number of sampled site days should not be considered a metric of how well the feature space is sampled.

Figure 1 Cross-validation sampling in meteorological space: number of site days contributing to sampling for NEE for the previous FLUXCOM RS + METEO ensemble (a) compared to the sampling of FLUXCOM-X-BASE (b) in environmental space of daily aggregated incoming shortwave radiation and VPD. Color corresponds to number of site days per bin in log scale. Only bins with at least 20 site days are shown.

The results from the 10-fold cross-validation showed an overall high performance with most fluxes and scales of variability having an NSE above 0.6 (Fig. 2). In terms of scales of variability across all fluxes, the monthly mean diel cycle (“diel”) and the daily median seasonal cycle (“seasonal”) were very regular patterns that the trained models reproduced best. Also, among-site changes (“spatial”, except for NEE) and monthly aggregated fluxes (“monthly”) were reliably predicted. Deviations from the median daily seasonality (“anom”) were only moderately reliable with NSE between 0.25 and 0.5. The XGBoost models did not succeed in accurately reproducing interannual variability of all fluxes and between-site patterns in NEE. Consistently across all scales, the net fluxes which are directly calculated (i.e., ET and even more so NEE) showed lower performance than their respective modeled gross fluxes (i.e., GPP and ET_T). Note that the cross-validation results from Fig. 2 cannot be quantitatively compared to previous cross-validation results in FLUXCOM as the training data are not the same. However, qualitatively the accuracy gradient among fluxes, as well as along scales of variability, corresponded to patterns identified in FLUXCOM and in comparable empirical modeling activities (Jung et al., 2011; Tramontana et al., 2016; Virkkala et al., 2021; Dannenberg et al., 2023).

Figure 2 FLUXCOM-X-BASE site-level accuracy of predicted fluxes in 10-fold leave-site-fold-out cross-validation in terms of NSE computed per site for a range of scales of variability. Scales of variability include the hourly timescale (“hourly”), daily (“daily”) and monthly (“monthly”) aggregated fluxes, between-site changes (“spatial”), monthly mean diel cycle (“diel”), daily median seasonal cycle (“seasonal”), deviations from the median daily seasonality (“anom.”), and interannual variability (“i.a.v.”). Boxes denote the range from the 25th to the 75th percentile of sites; whiskers extend 1.5 times the interquartile range from the 25th and 75th percentile of NSE across sites.

3.2 Global flux estimates

One asset of FLUXCOM-X is flexibility in the spatiotemporal resolution of the flux estimates. We produce X-BASE products at 0.05° spatial and hourly temporal resolution globally. Figure 3 illustrates the increase in spatial and temporal detail in X-BASE compared to RS (0.083°, 8-daily) and RS + METEO (0.5°, daily) using the example of NEE.

Figure 3 Resolution improvements for the X-BASE products compared to RS and RS + METEO: average NEE for an 8 d period in Europe in 2010 as estimated from the RS, RS + METEO and X-BASE setups (top panel), as well as snapshots of temporal trajectories of NEE in pixels closest to selected EC station locations (A: UK-Tad, B: DE-Hai, C: ES-LM1, D: FR-Pue). Negative values of NEE denote a CO₂ flux from the atmosphere to the land.

3.2.1 Net ecosystem exchange (NEE)

The X-BASE product estimates the global terrestrial NEE to be −5.75 ± 0.33 Pg C yr⁻¹ (2001–2020), with hotspots of strong CO₂ flux into the ecosystems in the tropical regions and temperate regions of North America and Europe (Fig. 4). In contrast to both RS and RS + METEO, India and some regions in central Sahel show prominent patterns of a mean CO₂ flux from the ecosystems to the atmosphere in X-BASE, corresponding mostly to crop-designated areas (Fig. C2). However, X-BASE global terrestrial NEE (−5.63 Pg C yr⁻¹) agrees well with the inversion estimates of OCO-2 (−4.12 Pg C yr⁻¹) and CarboScope (−3.46 Pg C yr⁻¹) over the common period (2015–2020). Here, GFED 4.1 (Randerson et al., 2017) emission estimates have been used to correct the inversion estimates for fire emissions.

Figure 4 Comparison of annually integrated NEE from X-BASE, RS + METEO with ERA5 forcing, and RS averaged over the period 2001–2020. The difference maps show the difference of the averages over 2001–2020.

Comparison with OCO-2 and CarboScope inversions also indicates a substantial improvement of the global mean seasonal cycle of NEE (Fig. 5) in X-BASE compared to RS and RS + METEO. The systematic bias present in RS and RS + METEO has essentially disappeared in X-BASE. The shape, and in particular the amplitude, of the global NEE seasonal cycle of X-BASE is more consistent with the inversions. The larger and more realistic seasonal cycle amplitude of global NEE in X-BASE originates primarily from improved and increased amplitudes in boreal regions. Interestingly, X-BASE suggests slightly larger NEE seasonal cycle amplitudes in temperate regions compared to the inversions. In seasonally dry regions, the timing of maximum uptake is consistent between X-BASE and inversions, while the peak of maximum net release is larger and delayed in the inversions. In Australia, the peak of CO₂ release to the atmosphere at the end of the year present in both inversions is not evident in X-BASE, which instead shows a relatively consistent CO₂ flux to the atmosphere throughout the year. In tropical regions, the patterns of seasonal variations are qualitatively consistent between X-BASE and the previous RS and RS + METEO products. The seasonal patterns in tropical regions are relatively weak overall and seem inconsistent between the inversions and X-BASE as well as among the inversions.

Figure 5 Seasonal and interannual variability of global NEE. Comparison of mean seasonal cycles (calculated over the common time period, 2015–2020) and interannual variability (2001–2020) of NEE estimated from CarboScope and OCO-2 inversions as well as FLUXCOM-X-BASE, FLUXCOM RS + METEO, and RS outputs. All products were integrated with a common mask that removes sparsely vegetated arid regions not predicted by RS and RS + METEO.

As shown in Fig. 5, the X-BASE product shows the same large underestimation of globally integrated NEE interannual variance as the previous RS and RS + METEO products. In terms of temporal trends, the X-BASE products show almost no change in annual NEE in time, which is in contrast to the RS + METEO (slight positive trend) and RS (slight negative trend) and more consistent with the CarboScope inversions (Table C2). However, as interannual variability was poorly reproduced even in the cross-validation (Fig. 2), trends in the X-BASE products should be taken with caution and interpreted with careful scrutiny.

3.2.2 Gross primary productivity (GPP)

X-BASE estimates the globally integrated GPP at 124.7 ± 2.1 Pg C yr⁻¹ on average in the time period 2001–2020. Globally integrated GPP over vegetated areas (RS and RS + METEO do not have estimates for non-vegetated areas) was approximately equal for X-BASE (121.9 ± 2.1 Pg C yr⁻¹) and RS + METEO (121.6 ± 0.4 Pg C yr⁻¹) but considerably higher than RS (113.2 ± 1.8 Pg C yr⁻¹) over the same period. In terms of regional patterns, X-BASE GPP consistently exceeds both RS + METEO and RS in temperate, boreal, and most subtropical ecosystems but is lower in sparsely vegetated (semi-)arid regions like southwestern North America as well as southeastern Asian croplands (Fig. 6). This qualitatively consistent pattern is only broken in the humid tropics, where X-BASE GPP is higher than RS but lower than RS + METEO.

Figure 6 Comparison of annually integrated GPP from X-BASE, RS + METEO with ERA5 forcing, and RS averaged over the period 2001–2020. The difference maps show the difference of the averages over 2001–2020.

Comparing the estimated trend over the last 2 decades, X-BASE GPP has a clear increasing linear trend of 0.34 Pg C yr⁻¹, which is slightly higher than the trend in RS (0.25 Pg C yr⁻¹; Table C2). In contrast, the RS + METEO product shows nearly no trend in annual GPP. The increases in both the X-BASE and RS products may be related to increases in surface greenness coming from variability in the remote sensing forcing data which are interannually dynamic in both products, whereas the remote sensing data were not interannually dynamic in the RS + METEO product, which instead used only the mean seasonal cycle of the remote sensing data. The magnitude of between-year changes in globally integrated X-BASE GPP is 0.575 Pg C yr⁻¹ over the years 2001–2020, which is about twice as large as RS + METEO (0.248 Pg C yr⁻¹) but only half the magnitude estimates in the RS setup (1.02 Pg C yr⁻¹; Table C2).

We further compared the temporal trajectory in GPP estimates against TROPOMI SIF as an independent proxy for GPP dynamics (Fig. 7) at a temporal resolution of 16 d. The temporal variability of X-BASE GPP strongly agrees with that in TROPOMI SIF, with squared Spearman correlation values (denoted as R²) of the time series above 0.85 across most of the vegetated land surface (Fig. 7a). The only exceptions are regions with no or very small variability in both GPP and SIF, such as in evergreen tropical ecosystems in South America, Africa, and southeastern Asia; sparsely to non-vegetated areas due to aridity (e.g., Mexican and African deserts); or cold conditions (e.g., Canadian and Siberian subpolar regions). In inner Australia, despite being sparsely vegetated, variability between years is expected in GPP due to precipitation increases during La Niña years, which is however not reflected in the squared correlations. R² values for the deviations from the average seasonality (again computed with a temporal resolution of 16 d) show the same qualitative spatial patterns (Fig. 7b) but are overall lower with R² values between 0.55 and 0.8. Anomalies of X-BASE GPP and SIF agree best in eastern European temperate forests as well as in grassy and shrub ecosystems in eastern South America.

Comparison of the level of agreement of SIF and X-BASE with that of SIF and RS and RS + METEO, respectively, illustrates that X-BASE and RS GPP estimates have comparable consistency both for the time series (global area weighted mean R² values of 0.72 and 0.73, respectively) and anomalies (global mean R² values of 0.64 and 0.66, respectively). In contrast, the R² between RS + METEO and SIF is lower in both cases (R² values of 0.66 for the time series and 0.58 for anomalies). X-BASE GPP shows a higher agreement with SIF than RS both in terms of the actual trajectory and anomalies in evergreen tropical forests with no or only a very short dry season in the Amazon and Africa, as well as in fully humid parts of southeastern Asia (Fig. 7c, d). Improvements in X-BASE GPP compared to RS are also consistent in the very continental and polar tundra areas in eastern Siberia, northern Canada, and Alaska. Conversely, in arid steppe climates globally, X-BASE GPP variability agrees less with SIF than does RS GPP. X-BASE GPP variability is consistently and widespread much more similar to the variability in TROPOMI SIF than RS-METEO GPP. Increases in R² for X-BASE compared to RS + METEO are most pronounced in arid to semiarid ecosystems (large parts of the Caatinga and Gran Chaco regions in South America, steppe regions in Mexico, southern and eastern Africa, Australia, and central Siberia) as well as in global crop regions, especially for the deviations from the seasonality (albeit with the magnitude of R² change being quite variable between regions; Fig. 7e, f).

Figure 7 Similarity of temporal patterns between GPP estimates and TROPOMI SIF observations: R² (computed as the square of the Spearman correlation) between X-BASE GPP and TROPOMI SIF (Köhler et al., 2018) for the actual time series at a temporal resolution of 16 d (a) and anomalies from the median seasonality in both variables (b). Panels (c), (d), (e), and (f) compare the agreement between X-BASE GPP and TROPOMI SIF to the agreement between FLUXCOM GPP and TROPOMI SIF, where (c) and (d) refer to TROPOMI SIF and GPP from the RS + METEO setup, and (e) and (f) refer to the RS setup. All comparisons are done for time series with a resolution of 16 d for the common time period April 2018 to December 2020. SIF observations have been applied a correction factor to estimate daily average SIF before aggregation. Semitransparent areas mark pixels in which the correlation of at least one of the data sets is negative.

3.2.3 Water vapor fluxes

Globally integrated ET amounts to 74.7×10³ ± 0.9×10³ km³ yr⁻¹ for 2001–2020 (Table C1) for X-BASE, with the highest rates in the tropics (Fig. 8). Comparison of global totals for vegetated areas only (where all products give outputs) shows similar values for X-BASE (68.9×10³ ± 0.9×10³ km³ yr⁻¹), GLEAM (70.9×10³ ± 0.9×10³ km³ yr⁻¹), and RS + METEO (68.3×10³ ± 0.3×10³ km³ yr⁻¹) ET estimates, while the RS ET is more than 11 % higher (78.5×10³ ± 0.5×10³ km³ yr⁻¹; Table C1). Particularly in evergreen tropical ecosystems, X-BASE estimates a considerably lower ET than GLEAM, RS + METEO, and RS (Fig. 8). Furthermore, in the temperate and high latitudes of the Northern Hemisphere, annually integrated X-BASE ET is consistently lower than the other estimates, though the magnitude of the bias is smaller than in the tropical regions. The pattern is only reversed with higher X-BASE ET in the semiarid and arid ecosystems of the lower and middle latitudes, especially with respect to annual ET in RS + METEO and GLEAM.

Comparison to precipitation estimates shows that X-BASE ET greatly exceeds precipitation inputs over large areas, indicating a strong overestimation of X-BASE ET in many arid regions with sparse vegetation (e.g., the Sahara region; Fig. C3). While transport of water both laterally (i.e., runoff from other areas) and from deeper groundwater could cause ET to exceed precipitation inputs in some areas, the extent of area where ET exceeds precipitation (e.g., the entire Sahara region) and the magnitude of the excess ET (over 3 times precipitation inputs) indicates a major bias in these areas and is likely due to a lack of EC data in similar ecosystems. As a rough estimate, constraining the X-BASE estimates with precipitation (see Appendix C5) suggests about 4–6×10³ km³ yr⁻¹ of water is overestimated globally.

Figure 8 Comparison of evaporative flux estimates ET, ET_T, and ET_T $/$ ET from X-BASE and its difference with RS, RS + METEO, and GLEAM. ET_T is compared in the case of GLEAM but is unavailable in the previous FLUXCOM ensembles.

The globally integrated ET_T amounts to 42.6×10³ ± 1.0×10³ km³ yr⁻¹ (2001–2020) in X-BASE, resulting in an average ratio of transpiration to total evaporation of 57.0 % ± 0.6 % (Table C1). In contrast to ET, the ET_T estimates from X-BASE do not commonly exceed precipitation estimates (Fig. C3), which could indicate that because the water vapor flux is more tightly coupled with vegetation, the model is able to distinguish that no vegetation corresponds with no transpiration, which is not generally the case for non-transpiration evaporation. The RS and RS + METEO products did not produce ET_T estimates, so the comparison is limited to GLEAM (50.7×10³ ± 0.6×10³ km³ yr⁻¹), which estimates ET_T on average 17 % higher than X-BASE, with strong contributions from the evergreen tropics. Only in single semiarid regions, such as the northernmost Sahel as well as large parts of the South American Caatinga and Chaco regions, is this pattern reversed (Fig. 8).

Spatially, X-BASE-estimated ET_T $/$ ET exceeds 50 % in the majority of areas, with the highest values seen in the higher-latitude regions of Europe and Asia, as well as in subtropical ecosystems (Fig. 8). Arid regions with sparse vegetation show the lowest ET_T $/$ ET overall, with values generally below 20 %. With 71.4 % ± 0.6 % over global vegetated surfaces, GLEAM attributes about 10 % more of its ET to ET_T than does X-BASE (Table C1). Regionally, this difference can even reach up to 40 %, with the only exception being boreal forests and very dry ecosystems in the Sahel, the Arabian Peninsula, and central Asia (Fig. 8).

Trends in ET, ET_T, and ET_T $/$ ET are positive and exceed the trends seen in all other estimates over the years 2001–2020. Conversely, the magnitude of interannual changes in X-BASE ET, ET_T, and ET_T $/$ ET is mostly less than half the variability in GLEAM (Table C2). Low interannual changes are common to the RS and RS + METEO ET as well.

Figure 9 shows the temporal correlation at 16-daily temporal scale using GLEAM as a reference, showing overall high values of squared correlation between X-BASE and GLEAM ET and ET_T (top and bottom left). Notable exceptions with low correlations are areas with low variability in ET such as the arc of deforestation, very dry areas, and tropical evergreen ecosystems in Africa. Compared to RS + METEO and RS (middle panels, left column in Fig. 9), X-BASE ET temporal patterns are more similar to GLEAM ET in many areas, especially so in areas north of the arc of deforestation and parts of tropical evergreen areas in central Africa and southeastern Asia. Conversely, X-BASE ET agrees less well with GLEAM than RS or RS + METEO in the arc of deforestation itself, the eastern parts of the Amazon basin, and dry areas. The deviations from the mean annual cycle in ET and ET_T (right column) show overall lower correlations than the actual time series, with the highest agreement between GLEAM and X-BASE in large parts of the Amazon forest and central European ecosystems. X-BASE ET anomalies are much more strongly correlated with GLEAM ET than either RS or RS + METEO everywhere except for most (semi-)arid regions.

Figure 9 Similarity of temporal patterns of X-BASE and GLEAM terrestrial water fluxes in comparison to those from FLUXCOM: R² (computed as the square of the Pearson correlation) between X-BASE ET and GLEAM ET for the actual time series (left column) and anomalies from the median seasonality (right column). The middle panels compare the agreement between X-BASE ET and GLEAM ET to the agreement between FLUXCOM ET and GLEAM ET. The bottom panel shows the squared correlations between X-BASE and GLEAM ET_T but no comparisons to FLUXCOM because FLUXCOM did not include ET_T. All comparisons are done for time series with a resolution of 16 d and 0.05° for the years 2001–2020. Semitransparent areas mark pixels in which the correlation of at least one of the data sets is negative.

4 Discussion

We have presented qualitative and quantitative comparisons of the X-BASE products with previous FLUXCOM estimates and independent data sets where possible. Statements on the main causes of (dis-)agreements in the magnitude and spatiotemporal variability in the flux estimates remain speculative as detailed analyses would require further investigation in the form of experiments, which is beyond the scope of this article. In the following we discuss the prominent changes in NEE compared to FLUXCOM as a key scientific highlight in X-BASE and outline potential ways forward to tackle persistent challenges in empirical upscaling.

4.1 Higher consistency of NEE with atmospheric carbon cycle constraints

Although FLUXCOM-X follows the same fundamental approach as FLUXCOM, we find a substantial improvement of the magnitude of the annually integrated NEE of FLUXCOM-X-BASE over previous FLUXCOM products (Jung et al., 2020) when compared to independent estimates from atmospheric inversions. The mean global X-BASE NEE of −5.75 Pg C yr⁻¹ is slightly smaller than the inferred NEE of −3.92 Pg C yr⁻¹ (corrected for fire emissions based on GFED 4.1) from CarboScope. The remaining difference could easily be explained by carbon sources such as aquatic evasion and volatile organic compounds that are included in the atmospherically based estimate but not in eddy-covariance-based FLUXCOM (see Jung et al., 2020, and Zscheischler et al., 2017, for further discussion).

The improved global NEE of FLUXCOM-X-BASE originates most likely from enhanced quality of eddy covariance measurements in the training resulting from years of consistent improvements in instrumentation, method development, and quality assurance from the point of measurement to the synthesis data sets. Previous upscaling-based NEE products of Jung et al. (2011, 2020); Bodesheim et al. (2018) – all based on the La Thuile FLUXNET data set but varying with respect to machine learning methods, predictor variables, and temporal resolution – consistently estimated a nearly 3-fold larger global terrestrial carbon uptake compared to X-BASE. As discussed and speculated in Jung et al. (2020), La Thuile likely contained biased NEE measurements, in particular for some tropical sites (Fu et al., 2018), and together with the sparsity of data in the tropics, these biases were propagated to unrealistic tropical and global NEE estimates. The fact that we can now reconcile bottom-up global eddy-covariance-based NEE and estimates from top-down atmospheric inversions is a major achievement of the FLUXNET community. For context, 1 Pg C yr⁻¹ over the global vegetated area (145×10⁶ km²) corresponds to ∼ 7 g C m² yr⁻¹, which marks a challenge for achieving such an accuracy of mean NEE at any one flux tower site. The lesson learned here emphasizes once more that it is crucial to control for and minimize systematic biases of in situ eddy covariance measurements (Moncrieff et al., 1996).

The improved seasonality of X-BASE NEE, in particular for boreal regions, likely also results from enhanced information in the training data due to the hourly resolution. Similar improvements were observed by Bodesheim et al. (2018), who extended RS + METEO by training on half-hourly flux observations. The hourly resolution improves the seasonal high-latitude NEE likely due to better capturing the responses to light when day length varies strongly.

4.2 New opportunities by X-BASE products

The improvements of NEE make X-BASE attractive as a data-driven biogenic prior for atmospheric inversions (Munassar et al., 2022). Moreover, its hourly resolution facilitates better integration in inversion systems due to the accounting of diurnal flux and atmospheric transport variations, while its high spatial variations can provide patterns of flux variations that cannot be resolved by atmospheric constraints alone.

For the first time, X-BASE includes a global data-driven product of ecosystem transpiration. The estimated global ET_T $/$ ET ratio of 57 % is consistent with independent top-down assessments from isotope base methods (Good et al., 2015; Coenders-Gerrits et al., 2014) and past upscaling estimates (Wei et al., 2017; Schlesinger and Jasechko, 2014). The spatially and temporally high-resolution data-driven X-BASE ET_T product provides a valuable complementary perspective to simulations from process-based land surface models, which show large disagreements and often indicate global ET_T $/$ ET below 50 % (Berg and Sheffield, 2019; Miralles et al., 2016). This advancement opens new opportunities for large-scale studies of carbon–water relations on a diurnal timescale. The generation of the X-BASE ET_T product was facilitated by the development of site-level evapotranspiration partitioning methods (Nelson et al., 2018, 2020) underlining once more the importance of advances by the FLUXNET community for Earth system science.

5 Tackling persistent challenges

In addition to the improvements and opportunities that X-BASE brings through the enhancement in site-level training data quality and quantity as well as the higher spatiotemporal resolution, we find that some key issues previously identified in FLUXCOM (Tramontana et al., 2016; Jung et al., 2020; Bodesheim et al., 2018) persist in X-BASE. These include (1) the insufficient representation of water-related effects, (2) the limited predictability of the spatial patterns of mean NEE, and (3) severe limitations with respect to the variability between years and over decades. In the following, we outline what we consider the most promising ways forward to tackle these (partly interrelated) challenges in future developments beyond X-BASE: novel additional predictor variables, more in situ observations (especially for some key ecosystems), and innovative approaches to constraining machine learning models.

The persistent challenge and importance of capturing water effects on land–atmosphere flux variations is illustrated by the overestimation of mean ET in very dry, sparsely vegetated areas (Fig. C3), as well as by the poorer consistency of NEE seasonality with inversions in water-limited regions (Fig. 5). For GPP temporal patterns we find that X-BASE shows improved agreement with SIF in water-limited regions compared to RS + METEO (Fig. 7), which is likely because X-BASE uses concomitantly changing remote sensing observations opposed to a mean seasonal cycle only in RS + METEO. However, X-BASE shows deteriorated agreement with SIF when compared to RS, even though X-BASE was trained on hourly flux observations with improved coverage of dry conditions (Fig. 1). This decrease in performance indicates clearly the importance and uncertainty related to the predictor variable set for capturing water-related effects. Uncertainties due to water limitation effects are also seen in the divergence of ET and ET_T estimates in semiarid regions such as the Sahel and South American Caatinga (Figs. 8 and 9). Thus, there is considerable potential for advancements by including remote-sensing-based predictors on soil moisture, subdaily varying land surface temperature from geostationary satellites, SIF, and vegetation optical depth. Here, key challenges reside in achieving a sensible integration of flux observations with footprints that are much smaller than corresponding Earth observation products on the one hand and an integration of disparate Earth observation products with each other on the other hand. This is one research priority for further FLUXCOM-X developments.

Missing predictor variables is likely also a main reason for the limited skill of predicting between-site variability of mean NEE (Fig. 2), which can depend on legacy effects of disturbances and management that are not accounted for. Novel and complementary Earth observation products that characterize ecosystem structure and states related to biomass and canopy heights from SAR and lidar should help raise the accuracy of FLUXCOM-X-based mean NEE in future efforts. Also high-resolution Earth observation products such as from Landsat may help to better resolve spatial details. X-BASE shows a prominent pattern of carbon flux to the atmosphere in subtropical and crop-dominated regions of India and the Sahel, which emphasizes the need to improve in situ data coverage for agricultural systems, especially outside the temperate zone, and to include important metadata to better characterize these ecosystems and their site history. Despite the greatly reduced overall bias of mean NEE, we emphasize that X-BASE products are premature for diagnosing spatial variations of mean NEE.

The representation of longer-term dynamics remains an area with opportunity for improvements in X-BASE. Interannual variability is still poorly reproduced in cross-validation (Fig. 2), particularly for NEE, which is likely not only due to the complexity of processes shaping interannual variations but also due to temporal discontinuities in flux tower time series related to changes in instrumentation and factors like management (Jung et al., 2024) that are not accounted for. The complexities of relying on field-deployed instrumentation, together with the uncertainties related to linking satellite and flux data, cause poor signal-to-noise ratios and may impede good cross-validation results for interannual variability. Globally, comparisons of X-BASE with inversions reveal an underestimated interannual variance and a poor correlation for global NEE interannual variability (Fig. 5 and Table C2). Interestingly, X-BASE GPP shows improved correspondence with SIF anomalies compared to RS + METEO, especially in water-limited regions (Fig. 7), while no such improvement is evident for global NEE, which is likely due to the compensatory water effects in the global NEE signal (Jung et al., 2017). That a comparison of RS + METEO runs with different meteorological forcing data showed the weakest correspondence with inversion interannual variability when using ERA5 (Jung et al., 2020) explains the substantially better correlations of RS + METEO with inversions for NEE interannual variability in earlier studies (Jung et al., 2017, 2020), and it may also explain the poor correlation with X-BASE. Thus, testing whether alternative meteorological forcing data can improve global NEE interannual variability for X-BASE is an important next step. It remains unclear at this point whether accurate interannual variations at site level and globally can be achieved by the FLUXCOM approach in the near future. Additional constraints beyond FLUXNET such as atmospheric CO₂ measurements (Upton et al., 2023) or theoretical considerations in the form of hybrid (Reichstein et al., 2019) or deep learning models (Camps‐Valls et al., 2021) are promising, and we will foster such endeavors.

6 Conclusions

We presented X-BASE, a new set of global high-resolution data-driven products of land–atmosphere fluxes from the FLUXCOM approach. This represents a cornerstone of our developments of the FLUXCOM-X framework designed to explore and mitigate current limitations to upscaling from site to global scale. Improvements of the eddy covariance data facilitated reconciling estimates of global terrestrial net carbon exchange from X-BASE with top-down atmospheric inversions and allowed for the first time the generation of a global data-driven estimate of ecosystem transpiration. Beyond fostering all activities to enhance quality and coverage of available flux tower observations, most promise for future advancements by FLUXCOM-X relates to the synergistic exploitation of complementary satellite data streams to better capture effects related to water, site history, and management. This will be challenging as it requires developing strategies and methodologies to better integrate in situ flux observations and spaceborne Earth observations with very heterogeneous acquisition properties and with spatial resolutions that are often very coarse compared to flux tower footprints. The recent de-orbiting of the TERRA spacecraft requires employing alternative satellite missions where practical issues of data acquisition and conceptual issues related to temporal consistency and reduced overlap with FLUXNET records pose imminent challenges. With FLUXCOM-X we have prepared the ground for tackling these challenges, which can facilitate up-to-date and accurate flux estimates and thereby contribute to increased understanding of the Earth system in the future.

Appendix A: Naming conventions

Table A1 Names of products and modeling environments may create confusion, especially given the legacy of FLUXCOM. We therefore will follow the following naming conventions and require potential users of products and collaborators to also follow this convention.

Appendix B: Details on processing of Earth observation data

B1 Dynamic quality control and cutout size

The conditions in the pixels around a given EC station should best represent the conditions of the land surface in the area where the actual fluxes originate from. Given that the actual flux footprints are not generally available or computable for lack of critical information, we assume that the pixel containing the actual EC station (the “tower pixel”) is most representative of the dynamics of the area of influence on a tower. However, data availability and quality in the tower pixel are often insufficient. An iterative approach therefore selects both the cutout size and the strictness of the BRDF inversion quality from within defined bounds in a way that maximizes data availability and that ensures representativeness of the spatially averaged time series for the given site at the same time. In more detail, we start with a strict criterion for BRDF inversion quality (BRDF_Albedo_Band_Quality_Bandx flag in MCD43A2 ≤2, meaning only full inversions). Then, three options regarding the cutout size are considered:

1. A cutout option only considers the tower pixel, i.e., the single pixel closest to the tower location.

2. Those 20 % of pixels within 4×4 km² around a tower that are best correlated with the tower pixel are linearly regressed against the tower pixel and subsequently spatially averaged.

3. The 25 % of pixels within a 4×4 km² area that are closest to the tower are averaged with the inverse of the distance to the tower as weight.

Table B1 Land cover intermediary classification encoding for MCD12C1 classifications where the classification has a specific feature (value = 1.0), might have a specific feature (value = 0.5), does not have a specific feature (value = 0.0), or is unknown (value = −1.0).

The criteria for selection between options A–C are based on the number of available good-quality observations n in the resulting spatial average time series per site as follows:

if (n_A >= 60 %) and (n_B <=  70 %),
    select A;
elif (n_A >=  60 %) and (n_B >=  70 %),
    select B;
elif (n_A < 60 %) and (n_A > 15 %),
    select B;
else
    select C.

If after the previous steps still less than 40 % of good-quality observations outside of snow-covered times are available in the resulting average time series for a given site and index, the BRDF inversion quality threshold is relaxed to also allow magnitude inversions (MCD43A2 BRDF inversion quality flag ≤3), and the procedure to select the pixels contributing to the average described above is repeated. Consequently, the size of the area that a MODIS reflectance time series represents varies between sites and so does the BRDF inversion quality.

For the global gridded MODIS data, the BRDF inversion quality is consistently selected as ≤2 or ≤3 based on the number of available good-quality observations in a pixel.

B2 Details on the treatment of land cover information

Land cover information was passed through an intermediary classification system both to act as an encoding mechanism and to allow for arbitrary links between classification schemes. Rather than simple true/false classification for each category, different attributes are classified based on whether the classification has a specific feature (value = 1.0), might have a specific feature (value = 0.5), does not have a specific feature (value = 0.0), or is unknown (value = −1.0). In the specific case of the MCD12C1 classification scheme, the conversion is as shown in Table B1.

B3 Overview of XGBoost hyperparameters used

Table B2 Hyperparameters used for the XGBoost models.

Appendix C: Additional results

C1 Additional cross-validation results

Figure C1 Cross-validation sampling in meteorological space: number of unique sites contributing to sampling for NEE for FLUXCOM RS + METEO (a) compared to the sampling in the X-BASE setup (b). Color corresponds to number of unique sites per bin in log scale.

C2 Large carbon release in tropical croplands

Figure C2 Large carbon release in tropical croplands.

C3 Global magnitude of all fluxes

Figure C3 Potential ET overestimation based on the ratio of estimated ET to precipitation from the Global Precipitation Climatology Centre (GPCC; Schneider et al., 2022).

Table C1 Global magnitude of all fluxes. The first column is the globally integrated flux for all areas including sparsely vegetated dry areas from 2001–2020. The second column includes a common mask, which removes sparsely vegetated areas that are not computed for the RS and RS + METEO products. Values reported after the ± correspond to the standard deviation across years.

C4 Linear trends and interannual variability for all fluxes

Table C2 Long-term variability of fluxes. The first column is the linear slope of annually integrated fluxes over the years 2001–2020. The second column is computed as the standard deviation of annually integrated fluxes after the trend is removed.

C5 Potential overestimation of ET in dryland areas

Maps in Fig. C3 show the extent to which ET and ET_T exceed precipitation as the ratio of the total of each flux to total precipitation (from GPCC Schneider et al., 2022). Overall, X-BASE ET largely exceeds precipitation in most dry, sparsely vegetated areas, indicating overestimation. In contrast, ET_T does not show such extensive overestimation, limited instead to only smaller regions of the Sahara.

The amount of overestimation of X-BASE ET can be roughly estimated by replacing areas where annual ET exceeds precipitation inputs with the corresponding annual precipitation inputs for each grid cell, i.e., replacing areas where the ET $/$ precipitation ratio is more than a threshold with the precipitation rather than the estimated ET. Using thresholds from 1.25 to 2.5 gives an excess of ET (i.e., original ET minus precipitation corrected) from 3.9×10³ to 6.1×10³ km³ yr⁻¹.