Globally Scalable Approach to Estimate Net Ecosystem Exchange Based on Remote Sensing, Meteorological Data, and Direct Measurements of Eddy Covariance Sites

Abstract

Despite a rapid development of Nature-Based Solutions (NBS) for carbon removal in recent years, the methods for evaluating NBS still have certain gaps. We propose an approach based on a combination of remote sensing data and meteorological variables to reconstruct the spatiotemporal variation of net ecosystem exchange from eddy-covariance stations. A Lagrangian particle dispersion model was used for upscaling satellite images and flux towers. We trained data-driven models based on kernel methods separately for each selected land-cover class. The results suggest that the proposed approach to quantifying carbon exchange on a medium-to-large scale by blending eddy covariance flux data with moderate resolution satellite and weather data provides a set of key advantages over previously deployed methods: (1) scalability, achieved via the validation design based on a separate set of eddy covariance stations; (2) high spatial and temporal resolution thanks to the use of Landsat imagery; and (3) robust and accurate predictions due to improved data quality control, advanced machine learning techniques, and rigorous validation. The machine learning models yielded high cross-validation results. Stratification that uses separate Fluxnet stations for each fold of validation ensures that the models are accurate across the area covered by the Fluxnet sites. Overall, we present here a globally scaled technology for the land sector based on high resolution remote sensing imagery, meteorological variables, and direct carbon flux measurements of eddy covariance flux stations.

1. Introduction

Nature-based solutions (NBS) are increasingly recognized as an essential tool in combating global warming as they remove greenhouse gases (GHG) from the atmosphere [1], reduce GHG emissions through sustainable land use practices, and store carbon in the soil and biomass [2]. While NBS are becoming quite advanced and continuously going through significant developments, the methods for accurately quantifying their resulting benefits are still lagging [3]. A vast majority of the existing projects use biometric soil or tree survey data to estimate carbon stock changes [4]. These methods provide solid estimates of the sequestered carbon but come with multiple important limitations. They typically describe either soil layer or partial canopy layer (e.g., above ground), with time resolution limited to multiple months and years and spatial resolution limited to a particular plot where such measurements have been performed. As a result of such limitations, these methods may not be best suited for large-scale ongoing NBS quantification [2]. Moreover, the measurements required by these methods for a meaningful time-resolved, large-scale assessment are expensive and unaffordable for most users.

A pure modeling approach employing either Earth system models or terrestrial biosphere models has been used as an alternative to direct measurements on a large scale [5,6], but the spatial resolution of these models is too coarse, while the uncertainty remains high due to their inability to respond timely to actual changes in the ecosystem. Likewise, satellites that measure carbon dioxide (CO₂) concentrations and are used in inverse modeling to map GHG emission and sequestration [7,8] have coarse spatial resolution and hence can hardly be used for NBS assessment at the farm or forest stand level. Optical, radar, and lidar satellites provide a scalable solution for mapping soil carbon [9] and carbon fluxes at a much finer spatial resolution [10,11] and do have the potential to help large-scale space- and time-resolved NBS assessments if anchored in accurate ground measurements. The current remote-sensing-based studies typically use allometric/biometric data for machine learning (ML) model training [2]. The amount of ground data required for capturing the landscape heterogeneity and temporal variation is challenging to sustain and scale. The models thus inherit the disadvantages of the biometric measurements.

There is, however, another potential solution that combines many advantages of the above methods while overcoming some of their limitations and can result in a high spatial resolution, scalability to the regional level and practical affordability. This can be achieved by blending the direct highly time-resolved small-scale observations from the numerous eddy covariance flux towers [2] and the less direct but large-scale remote sensing methods. Considerable efforts have been made in this area lately to improve the quantification of the spatiotemporal CO₂ patterns, employing multiple different approaches vs. existing estimates to help resolve the majority of issues with current carbon accounting methods in the land sector. For example, [12] evaluate how seven light-use-efficiency (LUE) models captured the spatiotemporal gross primary productivity (GPP) variations from the LaThuile FLUXNET dataset. This study used data from 157 eddy covariance (EC) towers located in the six major terrestrial biomes (evergreen broadleaf forest (EBF), deciduous broadleaf forest (DBF), mixed forest (MF), evergreen needle leaf forest (ENF), shrubland (SHR), and grassland (GRA)). None of the seven models matched well with estimated GPP based on overall EC measurements. The models explained 41–57% GPP variations in the study sites. Reference [13] evaluated four LUE models for GPP estimation and compared them with data obtained by EC in 51 sites with diverse vegetation types. The average daily $R^2$ for all models and all sites was about 66%. Reference [14] used three different models to estimate GPP in a local deciduous forest in the US in combination with flux measurements and remote sensing data. Comparisons were made yearly, but on average, the correlation coefficient was much better, at 0.91. Reference [15] used two enhanced vegetation index (EVI)-based models to estimate the GPP in 15 sites across the USA. The model’s results vary according to the vegetation; in some cases, coefficients of determination can vary from 0.27 up to 0.94 for specific sites. Using remote sensing data, disturbance information, and a support vector regression model, [16] upscaled the CO₂ fluxes from a network of 21 EC towers to estimate the regional-scale Alaskan CO₂ budget from 2000 to 2011 and compared the results with an inverse model. The predicted values were consistent with the field observations ($R^2$ = 0.79 for NEE) and the magnitude of the inverse model. However, the interannual variations of the results achieved by the remote-sensing-based model and those of the inverse model were inconsistent and therefore require additional validation. Reference [17] used three remote sensing data-driven models to evaluate and understand GPP estimations in 119 EC sites in the northern hemisphere. Considering all the study sites, the mean $R^2$ was 0.63. However, the values can vary significantly depending on the model used and the vegetation, and in some cases, they can reach 0.78. Usually, the highest values were found for deciduous broadleaf and mixed forest. Reference [18] mapped the soil organic carbon (SOC) in Italy using ground observations and ancillary remote sensing data. Compared with all observations, the geographically weighted regression using MODIS normalized difference vegetation index (NDVI) obtained an $R^2$ of 0.45 for SOC. Reference [19] conducted comprehensive work generalizing and mapping NEE, GPP, and respiration data for all types of classification as a collaborative FLUXCOM initiative. The overall coefficients of determination for GPP and respiration were more than 0.71 and 0.64, respectively, while about 0.48 for NEE.

As seen in the examples above, despite the significant advances in the use of remote sensing for CO₂ flux estimations, the results point to mixed success and lack of consistency, so there is room for major improvements. The key areas of such improvements include increasing spatial resolution, capturing temporal variability, improving the data quality (primarily by reducing uncertainties and increasing consistency) and developing better algorithms to enable the simultaneous use of multiple relevant parameters and different spatial and temporal scales. In order to achieve such improvements, we propose to rely on the significant advances which have been made in recent years in the fields of machine learning, artificial intelligence, and neural networks to extract patterns and insights from remote sensing data, processing large volumes of remotely sensed data with minor human interference [20]. Here, we present the successful use of one such approach for globally scaled remote-sensing-based carbon accounting models anchored by the direct carbon flux measurements using eddy covariance flux stations. We use input data on the net ecosystem exchange (NEE, [21,22,23,24]) as targets for the training process and then combine the remote sensing and meteorological variables [19] as predictors to reproduce variations of NEE on a daily and monthly scale to demonstrate the increased predictive power of the proposed approach over more traditional methods.

An overarching goal of this work is to bridge the domains of science and technology by integrating multi-scale remote sensing and eddy covariance measurements to support projects aimed at avoiding GHG emissions or to remove CO₂ from the atmosphere. For this purpose, we target the following methodological improvements: (1) increasing spatial resolution, (2) capturing temporal variability, (3) improving the data quality (primarily by reducing uncertainties and increasing consistency), and (4) developing better algorithms to enable the simultaneous use of multiple relevant parameters and different spatial and temporal scales.

2. Data and Methods

2.1. Data

2.1.1. Eddy Covariance Data

The eddy covariance method is a standard method to directly measure fluxes between ecosystems and the atmosphere by computing the covariance between the vertical wind velocity and the entity of interest [21,25,26]. Such measurements are made routinely on all continents [27,28,29], with many sites linked across a global network called FLUXNET, dispersed across the world’s climate zones and representative biomes [29,30,31,32,33]. FLUXNET teams standardize the data format, perform uniform quality checks, and use highly vetted gap-filling and flux partitioning, with the newest product version being the FLUXNET2015 dataset (we used data for 2000–2014), which includes over 212 sites [32] (see Figure 1).

We used data from 171 Tier 1 sites with a data policy consistent with the Attribution 4.0 International (CC-BY-4.0) data usage license. Daily net ecosystem exchange (NEE) was filtered using a variable friction velocity threshold for each year and quality flags indicating the measured and good quality percentage (data higher than 0.8 were used, indicating percentage of measured and good quality gap-fill data). The following land-cover classes based on the International Geosphere-Biosphere Programme (IGBP) classification were used: cropland (CRO); grassland (GRA), which includes grasslands, closed shrublands (CSH), and open shrublands (OSH) [34]; evergreen needleleaf forest (ENF); deciduous broadleaf forest (DBF); evergreen broadleaf forest (EBF); mixed forest (MF); wetland (WET); and savannas and woody savannas (SAV).

2.1.2. Remote Sensing Data

We used Landsat surface reflectance (SR) and land surface temperature (LST) as well as a number of indices derived from the former as the main predicting features for our machine learning model. Landsat 5, 7, and 8 (TM, ETM+, and OLI respectively) cover our study period and have a 30 m spatial resolution. USGS provides open access to their Landsat Collection 2 Level-2 Science Products (L2SP) through Amazon s3 cloud service. L2SP data are geometrically and atmospherically corrected using the LEDAPS algorithm for TM and ETM+ and LaSRC for OLI imagery and thus prepared for scientific analyses. USGS Landsat was accessed from https://registry.opendata.aws/usgs-landsat (accessed on 1 February 2022).

We filtered out imagery with cloud cover above 75% of a scene. To ensure consistently high data quality, we masked out all possible artifacts as marked in the quality assessment bands provided with the product. We also left out high- to low-confidence clouds and cloud shadows detected by the CFMask algorithm [35].

We utilized SR from three visible, near-infrared (NIR), and short-wave infrared (SWIR 1 and SWIR 2) spectral bands (bands 1–5 and 7 for TM and ETM+; and 2–7 for OLI). These bands are well-suited for assessing vegetation and soil state and, consequently, together with LST, can serve as accurate predictors of NEE.

In addition, we selected two spectral indices and two tasseled-cap components [36] to map proxy land surface parameters to predict NEE better. Soil-adjusted vegetation index (SAVI) and tasseled-cap greenness (TCG) reflect vegetation intensity and chlorophyll content. In contrast, normalized-difference moisture index (NDMI), normalized difference vegetation index (NDVI), enhanced vegetation index (EVI), and tasseled-cap wetness (TCW) reflect moisture content in vegetation in our areas of interest [37,38].

$$TCG=-0.1603·blue+0.2819·green-0.4934·red+0.7940·NIR-0.0002·SWIR-0.1446·SWIR2$$

(1)

$$TCW=0.0315·blue+0.2021·green+0.3102·red+0.1594·NIR-0.6806·SWIR-0.6109·SWIR2$$

(2)

$$NDMI=\frac{NIR-SWIR}{NIR+SWIR}$$

(3)

$$NDVI=\frac{NIR-red}{NIR+red}$$

(4)

$$SAVI=\frac{NIR-red}{NIR+red+0.5}·1.5$$

(5)

$$EVI=2.5·\frac{NIR-red}{NIR+6·red-7.5·green+1}$$

(6)

where $blue$, $green$, $red$, $NIR$, $SWIR$, $SWIR2$ are Landsat spectral bands of the corresponding wavelengths.

2.1.3. Meteorological Data

In addition to the remotely sensed data characterizing the vegetation cover, the performance of NEE is influenced by weather parameters such as temperature, dew-point temperature, solar radiation, and precipitation. The global dataset AgERA5 was chosen as the source of this data [39].

This service is based initially on hourly European Centre for Medium-Range Weather Forecasts Reanalysis version 5 (ECMWF ERA5) data at the surface level, which are available at 30-km spatial resolution (about 0.28125°). Data were processed to daily temporal resolution by summing or averaging depending on a particular parameter and corrected to a more refined surface description at a 0.1° spatial resolution. Aggregated data follow a local time zone definition and include several major agronomic parameters. The correction to the 0.1° grid was performed by applying grid and variable-specific regression equations to an ERA5 dataset interpolated at a 0.1° grid. The equations were trained on operational ECMWF’s high-resolution atmospheric model (HRES) data at a 0.1° resolution.

The following set of parameters was chosen for our purposes: 24 h mean of 2 m air temperature ($TemperatureAir$), 24 h mean dew-point temperature ($DewPoint$), precipitation flux (sum of total precipitation per day, $PrecipitationFlux$), and surface solar radiation downwards (sum of solar radiation per day, $Solar$). A description of the additional meteorological features generated and used in this work is provided in Section 2.2.2.

2.2. Simulation Experiment

The eddy covariance method measures fluxes integrated over an upwind area of a specific size, known as a footprint [40,41]. The size of this area depends significantly on wind speed and direction. The NEE measurements from stations have a half-hour temporal resolution, which means that the station is sensitive to the area around at a distance that air mass has passed upstream during a specific half-hour. The linear scale of such an area may range from hundreds of meters to several kilometers.

Thus, in the case of 30 m resolution remote sensing data, we cannot directly associate the surface reflectance value in a band’s pixel (or calculated index) with the station inside and NEE estimation on this station for the specific date. In order to do so, we need to calculate the relevant weighted footprint [42,43] of the impact area around the station for the corresponding time scale (daily in our case, see Figure 2). To calculate this footprint, we used a meteorological model which can reproduce transport processes on a local scale depending on weather conditions on a given day. For this purpose, we used the Lagrangian particle dispersion model (LPDM) FLEXPART. We weighted six Landsat-derived spectral bands and indices and integrated them by a footprint probability density function estimated using the FLEXPART model in backward mode. There is no need to do such a calculation for meteorological features, as their spatial resolution is coarser (linear size of gridcell about 9–10 km or 100 km²).

2.2.1. Footprint Calculation with LPDM

We used the Lagrangian particle dispersion models (LPDM) FLEXPART [44,45] to estimate the footprint of Fluxnet stations. In Lagrangian models, individual particle trajectories (representing air parcel transport) are calculated following a predetermined atmospheric wind field and other meteorological variables. LPDMs can run forward (source to a receptor) or backward (receptor to source) modes.

We used the FLEXPART version 8.2, which utilizes National Center for Environmental Prediction Final (NCEP FNL) operational global meteorological data on a regular latitude–longitude horizontal grid with 0.25–1.0 degree spatial resolution on hybrid sigma-pressure vertical levels and temporal resolution of 3–6 h. The model output grid has a spatial resolution on a latitude–longitude grid with 0.002 degrees (about 150–200 m), and the smoother result, the model synchronization time step, was 15 s.

We ran the model in a backward mode to obtain the weighted footprint for a specific station, setting station parameters such as location coordinates and above-ground height as receptor position. The duration of trajectories was selected corresponding to the original EC system time scale, which is 30 min for which the residence time is calculated. Residence time is a regular output parameter of FLEXPART model in the backward mode. This model setup was run six times a day (i.e., every 4 h) for every station for every date where satellite images were available to obtain a daily footprint.

One of the regular output parameters of the FLEXPART model in backward mode is residence time, which is measured in seconds. That gives an understanding of how long one particle (effective air parcels), or group of particles, has “lived” in a particular grid cell, which allows us to evaluate the receptor’s sensitivity to the surface grid cells. Thus, the resulting residence time map can be interpreted as a weighting function (see Figure 3 right panels), which we need to rely on when processing the remote sensing data.

2.2.2. Feature Generation and Data Augmentation

We selected the resampled, weighted, and aggregated SR and the LST values from the Landsat bands, indices calculated from these bands, and the set of meteorological parameters taken from the reanalysis AgERA5 as predictors of a regression model to predict daily NEE, which is a product from Fluxnet stations (see Section 2.1.1).

In this study, a set of additional features was created to increase the flexibility of the regression model and expand the interactions between the base predictive features and the target variable. First, all possible ratios were created between 6 Landsat bands and LST. Thus, we have five ratios for blue band, four ratios for green band, three ratios for red band, two ratios for NIR, and one ratio for SWIR1 (e.g., blue/green, blue/red, blue/NIR, 15 features in total). When we computed vice versa, we obtained five ratios for SWIR2, four ratios for SWIR1, three ratios for NIR, two ratios for red, and one ratio for green (e.g., green/blue, red/blue, NIR/blue, 15 features in total). Another set of features is a list of ratios between bands and LST (blue/LST, green/LST etc. 12 features in total). As a result, 42 features were added to the predictive dataset. Second, from meteorological variables such as solar radiation, 2 m air temperature, and precipitation flux, additional features were also created. For solar radiation, cumulative and mean radiation over five different periods were calculated: four days, one week, two weeks, three weeks, and four weeks (ten features as a result), which may show delayed phenological effects on sub-seasonal and seasonal scales. The same set of parameters was prepared for 2 m air temperature. For precipitation flux, cumulative precipitation for the last one, two, three, four, and five days was calculated (five features). Thus, 25 more features were added to the base predictive features. In total, 87 features were prepared for further training.

Another issue in the dataset collection process is the availability of clear-sky observations from the Landsat constellation. After merging the dates with available satellite images and good quality NEE observations, fewer training points remain, mainly due to the images’ cloud cover and other artifacts. We decided to apply training sample augmentation techniques to solve this problem.

The procedure for increasing the training sample is based on the assumption that the surface reflectance on the day of satellite measurements and the neighboring days will be similar given similar weather conditions. Weather variables, i.e., precipitation, air temperature, short-wave, and incoming long-wave radiation, were determined using data from local meteorological measurements at each eddy-covariance station. From these four parameters, a weather profile was generated for the whole time series of daily NEE estimations, which was then compared with a similar profile on the day of satellite observations. If the $R^2$ value was greater than 0.9, it was assumed that the day had similar weather conditions, and images of the original day could be applied to such a day. As a result, we used the same Landsat scenes for remote sensing parameters at augmented days from the original observation date, while applying the weighted footprint function for that specific day. The time considered in this study was a week before the date of Landsat observation. As a result, the size of the training sample was increased by about a factor of 3.

2.2.3. Building Regression Model

We used Python as a programming language and performed the calculations on Amazon Web Services (AWS, Seattle, WA, USA). We processed 20,438 Landsat scenes for this study. Apart from satellite imagery processing, calculations associated with the feature selection process (from a few thousand features) were among the most computationally intensive.

Figure 2 demonstrates the data flow for the dataset creation process. The formulas for calculating weighted parameter values and a regression model are as follows:

$${\displaystyle F_{weighted}^j=\sum _{i=1}^N{w}_i\ast F_i^j}$$

(7)

$$NEE=H(F_{weighted}^1,F_{weighted}^2,\cdots ,F_{weighted}^M)$$

(8)

where $w_i$—weight value for the i-th grid cell (see Figure 3); $F_i^j$— value of j-th predictive feature for i-th grid cell; $F_{weighted}^j$—weighted sum for j-th feature; NEE—eddy covariance data; N—number of grid cells; M—number of features; and H—non-linear operator.

Figure 3 shows an example of original features calculated from Landsat data and results of LPDM calculations of the footprint for the up-scaling process. The weighted sum for each parameter may be interpreted as the dot product of the feature and footprint map. To train the model, we match the resulting input dataset (

Table 1. Example of the aggregated dataset for each feature for each date. Note that the values are scaled by 10,000.

Date	Blue	NDMI	SAVI
27 December 2009	344.9	327.6	2314.2
28 December 2009	346.8	345.2	2335.3
29 December 2009	338.3	315.5	2276.5
30 December 2009	342.1	317.5	2289.9
31 December 2009	355.2	355.5	2392.7
1 January 2010	351.3	356.4	2354.5
2 January 2010	349.3	325.1	2346.3
3 January 2010	351.4	361.2	2373.4
4 January 2010	355.4	363.2	2382.4
5 January 2010	352.2	348.9	2367.3

) with eddy covariance data (

Table 2. Example of “target” NEE from eddy covariance stations for each date.

Date	$\mathbf{NEE},\mathbf{gC}/{\mathbf{m}}^2/\mathbf{Day}$
27 December 2009	−5.77
28 December 2009	−7.43
29 December 2009	−10.24
30 December 2009	−5.49
31 December 2009	−7.83
1 January 2010	−5.47
2 January 2010	−5.24
3 January 2010	−3.87
4 January 2010	−5.52
5 January 2010	−3.45
⋯	⋯

).

We selected an algorithm based on the kernel method kernel ridge regression (KRR) [46,47]. This method is similar to the more widely known class of support vector machine (SVM) methods [48]. However, it uses a different loss function. The KRR method uses linear least squares with l2-norm regularization, meaning it has scale factor for squared feature coefficients as an additional term in the loss function. This allows for creating non-linear dependencies between the features according to the selected kernel and then searching for the optimal coefficients of these dependencies with control of the regularization process. The result of solving such a problem is a hyperplane in the multidimensional feature space that describes the dependencies of the target values on the selected features. This method works well with small datasets, which is the case in this study and is also faster than the original support vector regression (SVR).

KRR has three main hyperparameters: $\alpha$—the regularization parameter of the model (the higher the value, the stronger the regularization); kernel—kernel interaction between the features included in the model (see Equations (9) and (10)); $\gamma$—internal kernel parameter (see Equations (9) and (10)).

$$k(x,y)={\mathrm{e}}^{{-\gamma \parallel x-y\parallel }_1}$$

(9)

$$k(x,y)={\mathrm{e}}^{{-\gamma \parallel x-y\parallel }_2}$$

(10)

where e—exponent; x, y—input vectors; ${\parallel x-y\parallel }_1$—Manhattan distance between the input vectors; and ${\parallel x-y\parallel }_2$—Euclidean distance between the input vectors.

Another critical aspect of preparation is the feature selection process. As mentioned in Section 2.2.2, in addition to the primary remote sensing and weather parameters, 87 additional features in total were prepared. In addition, we used products between all meteorological and remote sensing features that led to a few thousand features.

For each IGBP class, a small set of features has to be selected to allow optimal prediction of NEE at the stations. The exact solution for finding the optimal combination of features can be obtained by a complete enumeration of all possible variants. The complexity of this approach is $2^n$, where n is a number of features, so the approach becomes inapplicable for a large number of features. The utilities in the package machine learning extensions (Mlxtend) [49] were used to select the features.

The selection process was divided into two stages. The first stage estimates mutual information for each feature in the list. A standard method from the scikit-learn library was used [50], which allows us to reduce the number of features to select. The second step is a sequential selection of features in the forward mode. That means the model starts to learn from one feature and sequentially adds features to improve the quality of the model. As a result of this selection, the feature sets for each IGBP class have been generated, as shown in

Table 3. Feature sets for predicting NEE for each IGBP class (detailed description in Appendix B).

Land Cover	List of Features
DBF	$TCG$, $Cum\text{\_}Temperature\text{\_}over\text{\_}5\text{\_}3week$, $blue/green$, $DewPoint$, $Solar$, $NIR/SWIR1$, $SWIR2/green$, $Cum\text{\_}precipitation\text{\_}3days$, $Cum\text{\_}Temperature\text{\_}over\text{\_}5\text{\_}2week$, $Cum\text{\_}Temperature\text{\_}over\text{\_}2\text{\_}1week$, $Cum\text{\_}Temperature\text{\_}over\text{\_}3\text{\_}3week$
GRA	$SWIR1/NIR$, $Cum\text{\_}Temperature\text{\_}1week$, $green/red$, $Cum\text{\_}precipitation\text{\_}2days$, $Mean\text{\_}Solar\text{\_}2week$, $Delta\text{\_}T$, $red/NIR$, $LST$, $NDVI$, $red/green$, $green/NIR$, $red/blue$, $SWIR2/SWIR1$, $Cum\text{\_}Solar\text{\_}4days$, $blue/NIR$, $Mean\text{\_}Solar\text{\_}4days$
CRO	$NIR/green$, $Mean\text{\_}Solar\text{\_}4week$, $Cum\text{\_}Temperature\text{\_}over\text{\_}5\text{\_}3week$, $NIR/SWIR2$, $DewPoint$, $blue/green$, $Cum\text{\_}Solar\text{\_}2week$, $Cum\text{\_}Temperature\text{\_}over\text{\_}2\text{\_}3week$, $SWIR2/NIR$, $Mean\text{\_}Temperature\text{\_}1week$
ENF	$Solar$, $SWIR1/NIR·Cum\text{\_}Temperature\text{\_}over\text{\_}5\text{\_}4week$, $Cum\text{\_}Temperature\text{\_}over\text{\_}3\text{\_}3week$, $TemperatureAir$, $NIR/SWIR2·Mean\text{\_}Solar\text{\_}2week$, $Delta\text{\_}T$, $blue/SWIR2$, $Cum\text{\_}Temperature\text{\_}over\text{\_}4\text{\_}3week$, $NIR/SWIR2·Solar$, $NIR·SWIR1$, $Cum\text{\_}Temperature\text{\_}4days$, $red/SWIR1$
MF	$TCG·Solar$, $Delta\text{\_}T$, $DewPoint$, $SWIR2/blue$, $blue/red$, $Cum\text{\_}precipitation\text{\_}5days$, $TCG·DewPoint$, $NIR/LST·Solar$, $NIR·Solar$
EBF	$LST/SWIR2·Solar$, $TemperatureAir$, $NIR/LST$, $Cum\text{\_}precipitation\text{\_}2days$, $LST/SWIR2·Mean\text{\_}Solar\text{\_}1week$, $Delta\text{\_}T$, $red/SWIR2·Mean\text{\_}Solar\text{\_}2week$
WET	$NIR·Mean\text{\_}Solar\text{\_}3week$, $Delta\text{\_}T$, $SWIR1/NIR$, $SWIR1/green$, $Cum\text{\_}precipitation\text{\_}3days$, $green·NIR$, $NDMI$, $NIR/SWIR1$, $NDVI·Mean\text{\_}Solar\text{\_}4week$, $TCG·Cum\text{\_}Temperature\text{\_}2week$, $TCG·Mean\text{\_}Temperature\text{\_}4days$, $TCG·Cum\text{\_}Temperature\text{\_}4days$, $NIR/SWIR2$, $DewPoint$, $TCG·Cum\text{\_}Temperature\text{\_}over\text{\_}4\text{\_}2week$, $TCG·Cum\text{\_}Temperature\text{\_}over\text{\_}3\text{\_}2week$, $green·SWIR1$, $TCG·Cum\text{\_}Temperature\text{\_}over\text{\_}2\text{\_}2week$, $SWIR2/NIR$
SAV	$NDVI·Solar$, $NIR/SWIR2·Solar$, $TemperatureAir$, $green/red$, $red/green$, $Cum\text{\_}Temperature\text{\_}over\text{\_}5\text{\_}2week$, $SWIR1/red$, $SWIR2/SWIR1$, $EVI·DewPoint$, $SWIR1/green$

.

Model training and hyperparameters tuning processes were based on N-fold cross-validation. In our case, N equals 3, 5, or 7 depending on the size of the available dataset per IGBP class. The original dataset was divided into N nearly same-sized folds with a unique set of Fluxnet sites inside. Then, the model was trained on $N-1$ folds and validated on the remaining fold after the procedure was repeated N times with different combinations of folds. Each fold contains a unique set of stations without overlapping with other folds [51], which ensures spatial scalability. Thus, N times independent validation was performed on the sets of stations that were not involved in the training of the model, which ensures spatial scalability of the model predictions. Similarly to [19], hyperparameters for each N fold were tuned separately for a fixed set of features. As a result, we have not just one model per land cover but an ensemble of models with their respective metrics trained on a different dataset configuration.

2.2.4. Metrics for Accuracy Assessment

Several standard metrics were used to assess the level of accuracy for the trained NEE estimation model: mean absolute error (MAE), root mean squared error (RMSE), coefficient of determination ($R^2$), and relative error (RE).

$$RMSE=\sqrt{{\displaystyle \frac{1}{N}\sum _{i=1}^N{(y_i^t-y_i)}^2}}$$

(11)

$$R^2=1-\frac{{\displaystyle \sum _{i=1}^N{(y_i^t-y_i)}^2}}{{\displaystyle \sum _{i=1}^N{(y_i^t-\overline{y})}^2}}$$

(12)

$$RE=\frac{{\displaystyle \frac{1}{N}\sum _{i=1}^N|y_i^t-y_i|}}{{\displaystyle \frac{1}{N}\sum _{i=1}^N|y_i^t-\overline{y}|}}$$

(13)

where N is the number of observation points to estimate; $y_i^t$ and $y_i$ are the i-th true and modeled NEE value respectively; and $\overline{y}$ is the mean NEE value.

3. Results

As a result of N-fold cross-validation for each land cover (Section 2.2.3), we received an ensemble of N different models trained on different independent datasets with the same inputs. Thus, we can use the mean value over the ensemble as a prediction.

First, selected metrics were applied to estimate the model’s accuracy by reproducing the original daily NEE values and then used to estimate the accuracy of the predicted NEE behavior. The results of the training and cross-validation experiments based on independent stations are presented in

Table 4. Cross-validation metrics (and standard deviation of the residuals for each metric) for daily NEE models per land-cover class.

IGBP	${\mathbf{R}}^2$	$\mathbf{RMSE},\mathbf{gC}/{\mathbf{m}}^2/\mathbf{Day}$	$\mathbf{RE}$
CRO	0.73 (0.17)	2.32 (0.67)	0.44 (0.19)
GRA	0.61 (0.04)	0.70 (0.35)	0.60 (0.05)
DBF	0.76 (0.05)	1.64 (0.16)	0.42 (0.04)
ENF	0.62 (0.05)	1.30 (0.22)	0.60 (0.06)
MF	0.55 (0.02)	1.25 (0.18)	0.64 (0.01)
EBF	0.42 (0.06)	1.28 (0.32)	0.73 (0.05)
WET	0.53 (0.06)	1.55 (0.19)	0.64 (0.08)
SAV	0.56 (0.03)	0.68 (0.14)	0.63 (0.02)

. The trained models can adequately reproduce the behavior of NEE at stations with different land cover (overall mean across land-cover NEE, $R^2$ > 0.6 with a range of 0.42 to 0.76). The lowest predicted values are for EBFs because such forests are mostly located in the tropics, where there are few Fluxnet stations, and there could be problems with cloud conditions during some seasons. The largest variation in error (standard deviation) is observed for CRO, possibly due to the difference in crop types and significant local anthropogenic influence on specific fields. An increase or decrease in NEE can be explained by the state of the biomass and the influence of agricultural management. Comparisons of observed values and predicted values are shown in Figure 4 and Figure 5.

To compare monthly NEE, additional daily footprint calculations similar to Section 2.2.1 were carried out. In this case, the calculation period of one week before and one week after the available satellite images were selected. Usually, periods of 8 days are used [19,34], but in our case, we decided to use a similar averaging backward and forwards. Furthermore, remote sensing data from the original image was used for the additional dates, but with meteorological variables and weighting functions corresponding to the actual date. Monthly model NEE data were then calculated for those months where it was possible to predict more or equal to 2/3 days of the month. The resulting monthly NEE was compared with monthly NEE values from stations for which the quality flag was greater than or equal to 0.8. The results of the comparison are presented in

Table 5. Cross-validation metrics for monthly NEE values per land-cover class.

IGBP	${\mathbf{R}}^2$	$\mathbf{RMSE},\mathbf{gC}/{\mathbf{m}}^2/\mathbf{Day}$
CRO	0.83	1.94
GRA	0.60	0.54
DBF	0.84	1.23
ENF	0.78	0.89
MF	0.70	0.75
EBF	0.70	0.80
WET	0.66	1.25
SAV	0.69	0.49

.

4. Conclusions

There has been a drastic rise in NBS development in recent years [3]. Although the methods for evaluating NBS are advancing, there are still certain gaps that need to be addressed [2]. Previous studies found that most of the current methods for quantifying NBS benefits use biometric soil or tree survey data to estimate carbon stock changes. These methods do not enable seamless large-scale assessments, while associated regular measurements make the methods financially unsustainable for most players in the land sector. Modeling approaches based on Earth system models or terrestrial biosphere models, as well as inverse modeling based on CO₂ concentrations data, have spatial resolution too coarse to be used locally. The current remote-sensing-based approaches mainly rely on allometric data as ground truth and therefore involve challenges with capturing landscape heterogeneity and temporal variation that negatively impact scalability. As the result of the analysis from the recent studies, it is fair to say that the approaches based on observations from eddy covariance towers and remote sensing are the best fit to enable affordable, high spatial resolution, globally scalable CO₂ estimates. We analyzed the studies representing a vast amount of scientific efforts targeted on these approaches and identified the significant challenges. Here, we presented an improved approach to quantifying carbon fluxes by blending eddy covariance data with moderate resolution satellite data and weather data. Our method (1) is scalable due to the validation design based on a separate set of eddy covariance stations, (2) features a high spatial and temporal resolution of the Landsat data and set of meteorological variables, and (3) delivers robust and accurate predictions due to improved data quality control, advanced machine learning techniques, and rigorous validation.

Fluxnet towers are expected to estimate NEE values for homogeneous land covers that are stable over time. Therefore, in order to apply our method the regions should be thoroughly controlled for these parameters. In addition, models can only be applied to the areas with known land cover. Application to areas with different land cover may skew the results significantly. At the same time, the eddy covariance system captures the net ecosystem exchange, which accounts for total CO₂ fluxes, including those related to the heterotrophic respiration process. So if the AI remote sensing model was trained with EC NEE, it consequently accounts for heterotrophic respiration.

In terms of global coverage, it is true that the sites we use may not be representative for all possible climate zones. However, since we use several validation folds that are blindly compiled (we do not know which station from which region is in which fold), we can assume that the regression we have constructed (ensemble of models) may be able to describe the underlying behavior processes for NEE. Such estimates then can be applied globally. A exact conclusion can be drawn using new EC data from independent sources and other climate zones.

Further in the Discussion, we suggest the key reasons for these improvements and also list our opinion on the current limitations of the proposed approach. Our methodology employed eddy covariance data from the Fluxnet network. Since the network has a standardized methodology for data acquisition and processing, we collected a robust and consistent training dataset, typically unavailable for more traditional carbon accounting methods based on biomass and soil sampling. The availability of a large amount of consistent and uniform data allowed for setting a fairly conservative threshold of 80% for the NEE data coverage and resulted in substantially reduced noise. On the remote sensing side of the approach, the Landsat archive provided a consistent and robust set of data with great temporal depth and spatial resolution. The FLEXPART model was instrumental in mapping the area that contributed to the NEE values captured by flux stations, helping to match these with remote sensing products better to account for landscape heterogeneity. Utilization of the KRR and its parameterization tailored to the goals of this study, and strict feature selection, enabled high prediction accuracy.

Still, the metrics of our regression estimators were not uniform and varied across the land-cover classes. The $R^2$ for the DBF was the highest, likely due to clear seasonality patterns and relatively low heterogeneity within the class. The same metric for the EBF class was the lowest, likely due to other carbon exchange process drivers and lower data availability due to high cloud cover in tropical regions.

Current results are presented for each land-cover class globally, providing visibility on the unified, scalable calculation approaches and the global accuracy baselines for various land covers (Figure 6). Overall, the achieved accuracy was in line with similar studies but outperformed most of those targeting global geographical coverage (overall land cover cross-validation $R^2$ = 0.73, RMSE = 1.53 gC/m²/day). Daily data accuracy has the $R^2$ range between 0.42 for EBF and 0.76 for DBF (monthly data accuracy is between 0.6 for GRA and 0.84 for DBF). The models have better accuracy in those parts of the world where the number of stations is higher and, conversely, slightly worse accuracy in areas of lower station density.

To resolve or minimize these deficiencies and improve accuracy for specific cases or regions, the current system allows additional training (not a local re-training) in the presence of the additional local flux observations, as suggested by [2]. Another way to improve the confidence of model data is to add other multispectral, hyperspectral, and radar satellite products, e.g., Sentinel-1, 2, or EnMap. It is expected that these new data streams can improve the accuracy of the daily data slightly but can contribute significantly to the monthly and annual data improvement.

Finally, since in this study we have trained N independent models for each land cover (where N = 3, 5, 7). In the future, it may be better to use ensemble predictions to allow additional characteristics such as mean, median, and standard deviation. These and other activities are currently in development but not advanced enough to be included in this manuscript.

Given the urgent need to decarbonize the atmosphere and the critical role of NBS in doing so, a fusion approach of eddy covariance and remote sensing data proved to be an efficient and cost-effective solution for quantifying these efforts. We hope our methods could contribute to mitigating global climate change in a cost-effective and transparent way.