Machine Learning-Based Prediction of Ecosystem-Scale CO₂ Flux Measurements

1. Introduction

Tower-based, ecosystem-scale CO₂ flux measurements quantify the exchange of turbulence flux of CO₂ (FCO₂, measured in $\mu$molm⁻²s⁻²) between the land surface and the atmosphere. Plainly, FCO₂ measures how much CO₂ is moving into or out of an ecosystem, per unit area and per unit time. During daytime hours, most ecosystems are a strong sink for CO₂ (negative FCO₂, following the micrometeorological sign convention), as they remove CO₂ from the atmosphere through the process of photosynthesis. By comparison, during the night, ecosystems are generally a moderate source of CO₂ (positive FCO₂), as they release CO₂ back into the atmosphere through the process of respiration. FCO₂ is measured using a method known as eddy covariance (EC) [3]. Eddy covariance measurements are continuous in time (24 hours a day, 7 days a week, 365 days a year) and are generally reported at an hourly or half-hourly temporal resolution. Global networks of eddy covariance flux towers collect in situ carbon flux measurements, providing information on photosynthesis dynamics across different ecosystems and under various environmental conditions.

Currently, FCO₂ is measured at hundreds of research sites across the United States, with 385 of these sites being members of the AmeriFlux Network [9,34]. While these measurements run continuously at high frequency (e.g. at 5 Hz), practical limitations such as technical failures, instrument malfunction, and the necessity for filtering out data with low turbulent conditions have lead to extended gaps in the collected data. Moreover, there is no attempt to standardize the measurements across sites within the AmeriFlux network, meaning that when measurements are available, they may be more or less reliable than another site. Together, these compromise the validity of the fluxes that are measured and reported without issue.

These are the issues we tackle in this paper: developing a consistent, robust, and explainable method for quantifying FCO₂. To do this we experiment with seven machine learning algorithms and 37 explanatory variables (environmental `drivers’) to make predictions about the half-hourly, daily, and annual FCO₂ in 19 different ecosystems across the continental United States. In particular, our contributions can be summarized as follows:

• A wide-ranging comparison of many common machine learning methods for predicting tower-based FCO₂;
• The discovery of a generalizable machine learning-based model that can predict FCO₂ to within 1.81$\mu$molm⁻²s⁻¹ of tower-based measurements;
• An open source gap-filled FCO₂ dataset covering 44 unique sites for free use by other researchers in the climate science community; and
• An open source code repository for reproducibility and wider implementation.

The remainder of the paper is organized as follows: Section 2 describes the background and purpose of the AmeriFlux Network, and provides the reader an explanation of the importance of quantifying carbon fluxes to the theory of Natural Climate Solutions; Section 3 discusses the state of existing work using machine learning for modeling CO₂ and other flux measurements; Section 4 details the data, algorithms, and structure of our experiments, while the results are presented in Section 5 and they are analyzed further in Section 6; finally, we draw conclusions, discuss limitations, and suggest future directions in Section 7.

2. Background Information

2.1. The AmeriFlux network

The driving motivation for the establishment of AmeriFlux almost three decades ago was to measure the carbon balance of different ecosystems, and more specifically, to better understand the distribution of CO₂ sinks and sources across the continent [49]. While this is an ambitious goal, from this perspective, the sampling provided by AmeriFlux is woefully inadequate–assuming that all 385 AmeriFlux sites are currently active (which they often are not), this coverage equates to approximately 1 flux measurement site every 25,000 km².

Therefore, extrapolation and upscaling from individual sites to fine resolutions and regional and continental scales must be done using either process-based or statistically-based models. The former approach is attractive because these simulation models are based on state-of-the-art understanding of how the carbon cycle works. However, parameterization and initial conditions remain outstanding challenges, and past model validation efforts have highlighted serious model errors. By comparison, the latter approach is unattractive because many of these statistical approaches are essentially black boxes from which it is impossible to verify process-level representation. Standardization of inputs for statistical models is also a challenge, and, to the best of our knowledge, validation of model predictions has generally not been conducted against independent datasets.

An extensive model-data comparison project of over 20 ecosystem models conducted under the North American Carbon Program found that process-based models generally performed poorly in representing site-level carbon flux dynamics across sites with varying land cover. Specifically, substantial model errors in representing FCO₂ were found at annual, seasonal, and diurnal time scales [10,44]; models misrepresented the inter-annual variability in observed CO₂ uptake [25]; models did not properly represent phenological transitions in Spring or Fall [37]; and models could not predict photosynthetic uptake within the uncertainty of observations [42]. These results lead to valid questions about the viability of using process-based models to evaluate natural climate solution strategies (discussed later in Section 6).

Statistically-based upscaling of FCO₂ began about two decades ago with the pioneering work of Papale et al. [35]. They used an artificial neural network, trained with CO₂ flux data from 16 measurement sites in Europe to calibrate a simulation model to predict CO₂ fluxes of European forests at 1 km resolution. Several years later, Xiao et al. [53] calibrated a modified regression tree model to FCO₂ measurements across the AmeriFlux network, using satellite observed greenness indicators, such as vegetation indices, leaf area index, and fraction of observed photosynthetically active radiation [24]. The sophistication of these kinds of upscaling efforts has matured over the last 15 years. The current state of the art is probably defined by the FLUXCOM project [23], which uses satellite remote sensing and gridded meteorological products to calibrate a model trained on FCO₂ measurements from sites around the world.

However, a challenge with past efforts to upscale site-level measurements is the lack of standardization in measurement protocols across sites. For example, across the AmeriFlux network, the choice of instrument setup and configuration, and even the details of flux data processing and corrections (which are critically important), may be different for each site. Furthermore, key instrumentation principles (e.g., open vs. closed path gas analyzer or sonic anemometer geometry), installation protocols (e.g., depth profiles of soil temperature and moisture measurements), measured and calibrated quantities (gravimetric vs. volumetric soil water content vs. soil water potential), and even units (hPa vs. kPa for vapor pressure deficit – easily converted, but also easily incorrectly reported or interpreted) are not consistent across sites. In particular, this lack of consistency of site variables across sites is a major barrier for any predictive modeling methods that use machine learning techniques.

The aforementioned inconsistency and variation in the AmeriFlux network’s data largely stems from its design as a “coalition of the willing”, where sites are setup and monitored by a large number of researchers (and consequently, research interests). Fortunately, FCO₂ data from the 47 long-term research sites operated within the National Ecological Observatory Network (NEON) are also contributed to the AmeriFlux data archive. NEON was specifically established to “collect long-term open access ecological data to better understand how U.S. ecosystems are changing” [4], and implicit in this mission statement is the need for standardization of measurement protocols and techniques across sites. This standardization opens up the possibility to use a machine learning algorithm to predict site-level FCO₂ without relying on gridded or reanalysis products as is necessary when using sites from AmeriFlux as a whole. Thus, the network of NEON sites represents an opportunity to train models on observational data across numerous sites which might be viewed as analogous to a model emulator [16]. The key difference being that this model is trained on real observations rather than the output of a simulation model.

NEON sites are strategically located, following a clustering algorithm to identify and group distinct regions of vegetation, landforms, and ecosystem dynamics into 20 different domains, as shown in Figure 1. Within each domain, at one or more monitoring sites, standardized measurements of environmental drivers (weather, solar radiation, etc.) are conducted along with ecosystem-level measurements of FCO₂ and other quantities measured by eddy covariance (e.g. sensible and latent heat fluxes).

As discussed in Section 4, our experiments leverage the reliability, consistency, and high quality of the data collected at NEON-based AmeriFlux sites to train machine learning models.

2.2. Natural Climate Solutions

Accurately quantifying the terrestrial-atmospheric exchange of carbon is vital to assessing the impact of environmental management projects and policies at all scales. Hemes et al. [19] argue that ecosystem-scale CO₂ flux measurements can play an important role in developing and evaluating climate mitigation strategies at the global level, while Hollinger et al. [20] noted the value of CO₂ flux measurements for quantifying the magnitude of carbon storage, on an annual basis, by a single evergreen forest. Both of these articles also highlight the value of accurate CO₂ flux measurements in the context of a theory known as Natural Climate Solutions (NCS).

NCS is a framework for adapting existing theory and knowledge of ecosystem science to mitigate the impact of anthropogenic climate change. It focuses on deliberate actions to manage, restore, and otherwise conserve ecosystems, to increase the quantity of CO₂ they remove from the atmosphere and store in slow-turnover carbon reservoirs, such as soil or woody biomass. While the role of terrestrial ecosystems in the global carbon cycle has been relatively well-understood for decades [43,52], the theory of NCS was first defined in 2017, in a presentation by a group of scientists and practitioners at the Proceedings of the National Academy of Sciences [17]. Since this time, support for natural climate solutions has gradually gained momentum [5,15], (potentially because efforts to reduce fossil fuel emissions have not yet been successful).

A recent work [13] defined the five foundational principles of NCS, with Principle 4 reading: `There are multiple potential NCS actions that can occur in a given landscape and quantifying the overall magnitude of opportunity can help to focus efforts on the actions that can offer the largest mitigation returns. However, appropriate accounting is required to ensure that NCS potential is consistently and clearly quantified.’ The authors argue that accurate carbon dioxide flux estimation is essential to the implementation and adoption of NCS, as it will help optimize the actions taking by governments and other land managers. The machine learning-based predictive modeling conducted in this paper aims to contribute to this important need.

3. Literature Review

The prediction of CO₂ is an important task in environmental sciences as rising levels of atmospheric CO₂ are the primary cause of climate change [26]. However, it is also a difficult predictive modeling problem, and relatively few studies have been conducted in this space, and even fewer using modern machine learning methods.

The majority of existing work using machine learning to predict CO₂ is focussed on estimating the atmospheric CO₂ concentration at various geographical scales (as opposed to our focus of predicting CO₂ flux). Alomar et al. [1] used extreme learning machines (a variation of feed-forward neural networks) to accurately predict CO₂ concentration based on a single site in Hawaii, while Hou et al. [21] used XGBoost to predict emissions in geographical regions of China. Conversely, Fang et al. [14] uses Gaussian processes, and Mardani et al. [31] uses a multi-stage neural network technique to predict national-level carbon emissions. A policy-based need for this modeling is evident in Zhang et al. [54], where a genetic algorithm was used to assess the impact of China’s ecological zones on its CO₂ emissions. Finally, Baareh et al. [2] take a time series forecasting approach to modeling CO₂ emissions with a neural network.

Model choice, hyperparameter tuning, and variable choice are always difficult in machine learning-based work, which is why Hamrani et al. [18] compared nine different machine learning algorithms to predict CO₂ for their specific agricultural sites, and Durmanov et al. [12] looked at the key enablers of greenhouse gases.

When considering the prediction of eddy covariance carbon dioxide flux with machine learning, the existing research either produced a low accuracy model, or a model that cannot generalize beyond the experimental site/s. For example, Tramontana et al. [48] predicted carbon dioxide and energy fluxes across global FLUXNET sites with 4 different algorithms, but the R² for net ecosystem exchange of CO₂ was less than 0.5. Alternatively, multiple works [11,41,50,55] all demonstrate high accuracy results using machine learning models to predict carbon dioxide on single or few experimental sites.

While most gap-filling techniques are process-based, Zhou et al.[56] used a variation of a random forest model to gap-fill extra long periods of missing values in carbon, heat, and energy fluxes. For further reading, a survey on using machine learning to predict various air pollutants including CO₂ is presented in [29].

4. Methods

Our experiments compared the performance of seven machine learning algorithms to predict half-hourly FCO₂ measurements collected between January 1st, 2016 and June 30, 2022. The experimental details are all provided in the following section and the code for the experiments is available at: https://github.com/jsl339/AmeriFlux.

4.1. Data

There are 47 NEON core terrestrial sites located across the U.S. and Puerto Rico, which strategically represent a range of vegetation, climate, and ecosystems divided into 20 different ecological domains as shown in Figure 1. Our experiments used data collected at 44 sites, as three sites–Marvin Klemme Range Research Station (OAES), Mountain Lake Biological Station (MLBS), and Puu Makaala Natural Area Reserve (PUUM)–were removed from the analysis due to inconsistencies in predictor variables, missing flux measurements, and errors arising during preprocessing. With these sites removed, our 44 sites represented 19 out of the 20 ecological domains (see [32] for general information about the data product).

We preprocessed the data using the R package REddyproc, as is the standard approach for gap-filling and u* filtering of carbon flux values. We used the U50 threshold to filter our u* values.

Table 1 shows a general explanation and summary statistics for the environmental drivers that we used as feature variables to learn our models. The data were sourced from 3 locations: AmeriFlux, the Phenocam Network, and MODIS satellite imagery [39,45,49].

Each site was assigned both a primary and secondary vegetation type from the following categories:

• Agricultural (AG)
• Deciduous Broadleaf (DB)
• Evergreen Broadleaf (EB)
• Evergreen Needleleaf (EN)
• Grassland (GR)
• Shrub (SH)
• Tundra (TN)

After preprocessing, our final dataset consisted of 961,340 observations unevenly divided among the 44 NEON sites.

4.2. Experimental Design

We compared the predictive performance of seven machine learning algorithms (explained in Section 4.3 below) in two experimental scenarios. The experiments employ cross validation, a common tool in machine learning experiments to ensure generalizability of the results. It does this by `holding out’ some data (called a `fold’), training the model without it, and testing on the held-out data. This gives an unbiased estimate of how the model would have performed on unseen data. This process is then repeated and averaged for robustness. When the data are split at random into k folds, this is referred to as k-fold cross validation.

In the first experimental scenario, we performed 10-fold cross validation on the data. This means that the data were randomly divided into 10 folds, with each containing approximately 10% of the data. The models were then trained using 9 folds (90% of the available data) and tested on the remaining fold. This process was repeated so that each fold was used in the training set 9 times and appeared as the test set once (see Figure 2a for an illustrated explanation). The performance of each algorithm was reported as the average across the 10 different runs. We note that the data were divided into the same 10 folds for each predictive algorithm.

K-fold cross-validation is a common technique in the testing and comparison of machine learning algorithms as it removes selection bias (whether deliberate or not), and demonstrates the ability of the models to generalize to unseen data [40].

In the second experimental scenario, which we will refer to as Leave-one-site-out cross-validation (L1SO CV), we began by partitioning the data by site, resulting in 44 uneven groups of data. We then employed a similar process to scenario one, where the models were trained on all-but-one group and tested on the remaining group (an example is shown in Figure 2b). This was repeated so that each site was used as the test data once, and therefore the stated performance metrics are the average of the 44 models fitted and tested.

The L1SO CV experiments present an inherently more difficult problem than the prior scenario as a predictive model significantly benefits from learning from data belonging to the test site. These experiments were included to replicate a situation where a site has no prior carbon flux recordings, i.e. it could be a new site or the instrumentation might not be functioning correctly. In addition, this experimental setup also tests whether we might be able to make a minimal set of measurements at a site with lower standardization in measurement protocols in order to predict the FCO₂. This would be helpful for carbon accounting purposes and nature-based carbon solutions, and also to enable a benchmark for land surface model simulations and checking existing datasets.

The performance of each model was assessed using 2 evaluation metrics—Root Mean Squared Error (RMSE) and the Coefficient of Determination ($R^2$). The RMSE is the square root of the average of the squared prediction errors over all of the data in the test set. Specifically,

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

Where y is the measured (true) value and $\widehat{y}$ is the predicted value for a test set of size N. Due to the squared component of the metric, the RMSE is sensitive to large errors in any of the individual predictions.

The $R^2$ evaluation metric is a measure of the goodness-of-fit of the linear model found by regressing the predicted values against the true values. It is calculated as:

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

In contrast to RMSE, $R^2$ is not sensitive to large errors in any of the individual predictions as it measures the amount of total variance accounted for by the predictions. When used together, these metrics complement each other and provide a more comprehensive picture of the performance of the algorithms. Each metric can then be analyzed on half-hourly, daily, and annual timescales. Ensuring accurate model predictions on an annual scale is important for reliable carbon accounting. However, it is also critical to evaluate model performance at finer temporal resolutions, such as half-hourly and daily scales, to ensure that our models produce accurate annual predictions for scientifically sound reasons. In order to produce meaningful predictions of annual sums of FCO₂ for each test site, we must first use models optimized in the 10-fold experimental setting to fill in missing FCO₂ values for each site before making predictions per site in the L1SO experimental setting.

4.3. Machine Learning Models

We compared the performance of 6 different machine learning models on predicting carbon dioxide flux. They are:

• Linear Regression (all predictors): This is a linear model including all of the variables using the maximum likelihood estimates for the coefficients. Linear regression assumes a linear relationship between the predictors and the response variable, which is unlikely in complex modeling problems, but does provide a baseline for the comparison of the performance of other models.
• Stepwise Linear Regression: This model began by testing for the most significant single variable in a linear regression model, and then iteratively added variables and tested for greatest improvement. A threshold number of selection variables was set to 15 for this forward selection technique. In this way, we simplify the basic linear regression model to find feature variables with greater importance for linear prediction.
• Decision Tree: A decision tree is a model based on recursively splitting the data on values of variables to maximize the difference between observations. Decision trees are most effective on problems where there is a non-linear relationship between the predictors and response variable [33,51]. The optimal tree depth was found to be 10 which was found through cross-validation.
• Random Forest: A random forest model [7] is a bagged ensemble of decision trees. The algorithm creates an uncorrelated forest of decision trees by using random subsets of features in each tree. When predicting a regression variable with a random forest model, the overall prediction is the average of the results of each of its constituent trees.
• Extreme Gradient Boosting (XGBoost): The XGBoost model [8] is a boosted ensemble of n underfit decision tree models. In practice, a decision tree is fit the to data and the errors in prediction are measured. Next, a second decision tree is used to fit the errors of the first tree. Then a third decision tree is fit to the errors of the second tree, and we continue until we have n trees in our ensemble. The optimal number of trees in our ensemble was found to be 2000. We also set the number of rounds for early stopping to be 50, and we used a learning rate of 0.05, max depth of 10, subsample ratio of 0.5, and subsample ratio of columns for each node of 0.45. Finally we used the histogram-optimized approximate greedy algorithm for tree construction to optimize our XGBoost model. All hyperparameters were optimized through 10-fold cross-validation using an exhaustive grid search.
• Neural Network (single-layer): A neural network is the sum of weighted non-linear functions of the predictor variables. This model is a single-layer neural network, with 256 neurons in the hidden layer, and uses a feed-forward architecture with ReLU activation. Early stopping was implemented to prevent model over-fitting, and training was performed with a data loader with a batch size of 128. The learning rate was set to 0.0003, and the best performance was achieved with no weight decay using the Adam optimizer. For more information on the mathematics of neural networks, see: [22,30].
• Deep Neural Network: The model uses the same mathematical structure as the single-layer neural network, but increases the number of hidden layers to 3, each consisting of 256 neurons. Compared to the single-layer neural network, the increased depth of the model increases the number of parameters to learn, meaning the model is capable of modeling more complex relationships, but also takes longer to learn from the data.

5. Results

5.1. 10-Fold Cross-Validation Results

The results for fitting each model and testing on each fold of the 10-fold cross-validaiton experiments are shown in Table 2 (RMSE). The XGBoost model and deep Neural Network were the only two models with a RMSE less than 2$\mu$molm⁻²s⁻¹. The strength of these models suggests that there are non-linearities in the the relationships between the environmental drivers and FCO₂. It is important to note that the XGBoost model outperformed our deep Neural Network at each stage throughout model development, and in addition the XGboost model requires significantly less training time than either neural network.

After determining the optimal algorithm, we used the trained XGBoost model to gap-fill all of the missing values for each of the 44 sites. The resulting dataset, consisting of 4,068,459 observations, is freely available at https://zenodo.org/records/10719776 for use by other researchers in the climate science community.

5.2. L1SO Cross-Validation Results

The results for fitting each model and testing on each site of the leave-one-site-out cross-validation experiments are shown in Table 3 (RMSE) and Table 4 (R²). Again, the XGBoost model was superior to all others with a mean prediction RMSE of 2.45$\mu$molm⁻²s⁻¹. This is 35% greater than the RMSE of same model in the 10-fold cross-validation experiments, demonstrating the substantial information the model gains from seeing data from the test site in the training set (as is the case in the 10-fold experiments).

The results also varied greatly between test sites–from a RMSE of 0.66$\mu$molm⁻²s⁻¹ up to 6.22$\mu$molm⁻²s⁻¹. The model performed best on Toolik (TOOL), as well as other sites with Tundra as the primary vegetation–Barrow Environmental Observatory (BARR), Healy (HEAL), and Niwot Ridge Mountain Research Station (NIWO)–suggesting that the environmental drivers for these sites are highly similar. Another justification for lower model RMSE across sites with Tundra primary vegetation is that these sites in general experience smaller magnitude fluxes. Random errors scale with flux magnitude, so it’s almost inevitable that sites with higher magnitude fluxes will have somewhat larger model-data mismatch.

The model performed worst on Lajas Experimental Station (LAJA), which is one of two sites in Puerto Rico, and together these two represent the only two sites with an evergreen broadleaf primary vegetation type. While we cannot separate the domain and primary vegetation effects here, we can say that our training data, which is mostly from the United States mainland, does not generalize well when predicting FCO₂ in vastly different climates and ecosystems.

A map of the average RMSE per domain is shown in Figure 3.

5.3. XGBoost Feature Importance

The XGBoost algorithm has a built-in method for calculating the importance of each feature variable based on the amount that each feature’s split point improves model performance. A plot of the twenty most important features for prediction is shown in Figure 4.

There are two input features that are noticeably more important to the model than others–EVI and net radiation. This is interesting as these are not measurements taken through site-level instrumentation, which suggests that we can learn a lot about the FCO₂ of a site just by knowing the vegetation greenness and the radiation environment of a site. Furthermore, six of the ten most important variables are continuous measurement variables, as opposed to the domain or vegetation categorical variables, meaning the model should generalize easier to any new sites of interest.

6. Discussion

6.1. Comparison of 10-fold and L1SO experimental results

By making predictions on each site in both 10-fold and L1SO contexts, we are able to gain greater understanding of model performance across the 44 NEON sites. We partitioned our model’s 10-fold RMSE by site, and treated a site’s average RMSE value as the irreducible error–that is, error that can be attributed to variability in the dataset, measurement errors, and the error inherent in using a model to predict a biological process. From there, we compare this irreducible error to the average RMSE values for each site obtained through L1SO CV experiments and therefore obtain an estimate of the amount of error attributable to testing on an `unseen’ sight, which we call the L1SO Remainder. A visualization of the baseline error and its corresponding L1SO remainder for each site, ordered by ecological domain, is shown in Figure 5.

The L1SO remainder gives us a reasonable way to identify sites that are difficult to predict without having that site’s data available in the training set. We identified five NEON terrestrial sites that had a L1SO remainder greater than 0.85. These sites are Guanica Forest (GUAN), Lajas Experimental Station (LAJA), LBJ National Grassland (CLBJ), Disney Wilderness Preserve (DSNY), and Wind River Experimental Forest (WREF). There are several reasons that can justify why these sites in particular may be difficult for a model in a L1SO scenario. Firstly, Guanica Forest and Lajas Experimental Station are the only two sites in Puerto Rico and in ecological domain 4. In addition, these two sites are the only two whose primary vegetation type is evergreen broadleaf (EB).

Wind River Experimental Forest is a site in Washington state, located in an old growth forest with very tall trees with a real summer dry-down that restricts FCO₂. Overall, Wind River Experimental Forest is a very unusual site in comparison with the other NEON terrestrial sites.

For each of these five sites, we created time series of predicted FCO₂ values and actual FCO₂values reported both in half hourly increments, and aggregated as an average for each day of the year, as well as a scatter plot of predicted FCO₂ vs. actual FCO₂ for analysis. We then compared these results to sites with the same primary vegetation types for which our model had superior performance. In the case of Guanica Forest and Lajas Experimental Station, since there were no other sites with the same primary vegetation type, both sites are included in Figure 6.

Steigerwaldt Land Services (STEI), Dakota Coteau Field School (DCFS), and Delta Junction (DEJU) were used as comparison for our other three sites representing primary vegetation types of DB, GR and EN respectively. These comparisons are found in Figure 7, Figure 8, and Figure 9.

Note that in most cases we observed large systematic errors in model performance for our 5 sites with the greatest L1SO Remainder values. For example, when considering scatterplots of predicted vs observed FCO₂ for LAJA and WREF, the slope of predicted vs observed FCO₂ is less than 1. At CLBJ, the magnitude of summertime uptake is under-predicted. At DSNY, the seasonality is represented well but there is a consistent offset of several $\mu$molm⁻²s⁻¹, with predicted values higher than measured values. By comparison, at STEI, DCFS, and DEJU, the magnitude and timing of predicted FCO₂ is much better.

What is interesting from this analysis is that even on sites with relatively high L1SO remainder, our model seems to do a good job on average predicting patterns and dips in daily average FCO₂. It appears that most of the errors associated with sites with large L1SO remainder can be attributed to the model being too conservative in its predictions, that is predicting values closer to zero than the true measured flux values. As seen by the right column of plots in Figure 7, Figure 8 and Figure 9, sites of the same primary vegetation type where our model had stronger performance seem to generally have less large positive and negative flux values. This makes sense, since our model learns to minimize prediction error, and since each error ends up being squared, predicting very large positive or negative values in general would be more heavily penalized. A good example of this fact can be seen in the half-hourly time series for Lajas Experimental Station in Figure 6. This site has a mix of large positive and negative observed flux values, and our model rarely made large positive or negative predictions. Compare this to a site like Delta Junction in Figure 8. Here, there are a number of large negative observed flux values, but not as many large positive values. Spikes in the negative direction are less erratic, and model predictions, as a result, more closely represented measured flux values. When looking at scatter plots of predicted flux vs. observed flux, one can see that the results for sites in the right hand column are more tightly clustered about the 1:1 line, resulting in higher $R^2$ scores.

6.2. Relevance of L1SO predictions for unseen sites

Historically, process-based models have been considered the “gold standard” for predicting ecosystem CO₂ fluxes. However, past model-data evaluation studies have shown that although process-based models can often predict daily- or sub-daily fluxes that agree reasonably well with measured values, model performance on longer time scales (seasonal, annual and inter-annual) is often quite poor [10,25,47]. Models that cannot accurately predict ecosystem carbon budgets on annual and inter-annual time scales are not likely to be useful for carbon accounting purposes or for developing strategies for nature-based climate solutions. This suggests that alternatives to process-based models are needed. While machine learning-based models have been used for flux upscaling for almost two decades [23,35,53], these analyses have generally attempted to extrapolate from individual sites to regions and continents using only remotely-sensed variables as drivers. While this strategy is intuitively appealing, it is unable to leverage the site-level characteristics that are undoubtedly relevant for making fine-scale predictions. Indeed, basic ecosystem theory suggests that without accounting for these site-level characteristics such as disturbance and land use history, it is impossible to predict ecosystem carbon balance. Notably, we found that site characteristics related to vegetation type, as well as to soils, were identified as among the most important features for predicting FCO₂. However, remotely-sensed variables from MODIS such as EVI and NDVI were found to be more important than site-level vegetation indices (e.g. Gcc, Rcc) derived from PhenoCam imagery. We can hypothesize that while PhenoCam imagery can provide phenological information at a fine spatial and temporal scale, it may be subject to issues related to the mismatch of footprints with eddy covariance flux measurements. In the case of heterogeneous landscapes, MODIS vegetation indices with larger spatial coverage may actually be more representative of seasonal variations in vegetation dynamics within the flux tower footprint.

Finally, we note that although site-level meteorological and environmental drivers (e.g. air temperature, relative humidity or VPD, soil temperature, soil moisture, and precipitation) were not ranked highly in terms of feature importance, this is not to say that these variables do not matter. Rather, it is likely that in the context of variation in FCO₂ from the Arctic to the Tropics, from winter to summer, and from day to night, that the additional information contributed by these variables explains only a small amount of the half-hourly variation in FCO₂, although it may contribute greatly to improved estimates of annual FCO₂.

A persistent challenge in estimating site-level carbon balance via FCO₂ measurements has always been that small but selectively systematic measurement errors in 30 minute data can accumulate to large errors in annual integrals [37]. In our machine learning approach, selectively systematic prediction errors could occur if important meteorological or environmental variables were not accounted for as covariates. Omission of these variables might do little to impact the $R^2$ calculated on 30 minute values but could seriously impact annual flux integrals. Adoption of model optimization criteria that place more weight on reducing selectively systematic bias (which might not even show up when bias is calculated over a multi-year dataset) and improving predictive power on annual and multi-year time scales could be important for further improving the application of machine learning methods to carbon accounting and nature-based climate solutions.

6.3. Leveraging Site-level Data when Standardized Model Inputs are not Available

Our feature importance plot (Figure 4) shows that, in spite of our assertion that site-level data are critical for correctly predicting ecosystem carbon balance, much of the information needed to predict half-hourly FCO₂ actually comes from variables that are already available from gridded land cover maps (i.e. vegetation type classifications), satellite data products characterizing phenology (i.e. EVI, NDVI), and basic energy balance data that are also widely available as satellite data products (e.g., net radiation). This suggests that there is the potential for leveraging the much greater abundance of AmeriFlux towers, (for which site-level measurements are not standardized, but still useful), together with key remotely sensed data products to generate an initial map of ecosystem carbon balance. This initial map, when fused with elements of the analysis presented here, could lead to a hybrid data product that leverages the sampling intensity of AmeriFlux and the standardized sampling of NEON. Development of a data fusion platform such as we describe here is beyond the scope of the present analysis, but it is potentially an exciting direction to be pursued in future research.

6.4. Annual carbon sums

For most sites, we managed to obtain low RMSE and high R² for predicting the measured half-hourly FCO₂, even in the L1SO analysis (Table 3 and Table 4). However, in the context of carbon accounting and nature-based climate solutions, it more important to know the overall carbon balance on an annual time scale. That is, we want to answer the question of how much carbon (if any) the ecosystem is removing from the atmosphere and putting into biomass and soil carbon on an annual basis. This carbon balance reflects the balance between plant photosynthesis (carbon uptake, or negative flux) and ecosystem respiration (carbon release, or positive flux). It is a challenge for models, either process-based or data-driven, to get the overall carbon balance correct because of the opposing nature of these processes on different timescales. For example, in most ecosystems there is a strong seasonal pattern of carbon uptake during the growing season and release during the dormant season. During the growing season there is also a diurnal pattern of carbon uptake during the day and release during the night. Annually, the difference between photosynthesis and respiration is much smaller (0-30%) than the flux associated with either of these two key processes.

A model that predicts the annual carbon balance for an unknown site would be extremely valuable if it successfully estimated the multi-year mean carbon balance. The model would be even more useful it if successfully represented the inter-annual variability in carbon balance. State-of-the-art process-based models have generally failed to meet either of these targets [25]. Our results show that across all vegetation types, annual sums predicted in the L1SO analysis did a surprisingly good job at hitting the first target (see Table 5). For 29 out of 44 sites (66%), the L1SO-predicted multi-year mean carbon balance was within ± 50gCm⁻²y⁻¹ of the “true” value estimated by gap-filling missing values in the CV analysis. This is quite remarkable given that the total uncertainty on the annual carbon balance, derived from gap-filled FCO₂ measurements, is typically estimated to be about ± 50gCm⁻²y⁻¹ [38]. However, for 7 of 44 sites (16%), the deviation between the L1SO-predicted multi-year mean and the “true” value was greater than 150gCm⁻²y⁻¹. Three of these were deciduous broadleaf forest sites, one was an evergreen needleleaf forest site, and one was a grassland site. We expect that there may be land use history, disturbance, or similar factors that might explain these deviations, but were not included in our model.

Annual sums predicted in the L1SO analysis also did a reasonably good job of representing the “true” inter-annual variability estimated from gap-filled time series. At more than a quarter of sites (12 of 44, 27%), the correlation of L1SO-predicted annual sums and the gap-filled annual sums was greater than 0.75, while for almost half of sites (21 of 44, 47%) the correlation was greater than 0.50. While these results are based on at most 5 years of data per site, they point to the enormous potential of machine learning to predict not only the long-term carbon balance of an unknown site, but even the inter-annual variation in that carbon balance. By comparison, it has been known for more than a decade that even the most sophisticated process-based models are unable to capture this inter-annual variability [6,36,46], despite accurately capturing the dynamics of “fast” processes operating on timescales of hours to days.

7. Conclusions

In this paper we showed the ability of machine learning-based models to make skillful predictions of tower-based CO₂ flux measurements. Specifically, we found that an XGBoost model trained on 37 environmental drivers, from 44 AmeriFlux sites, can predict FCO₂ at an unseen site to within an average error of 2.45$\mu$molm⁻²s⁻¹. Furthermore, this error reduces significantly–down to as little as 0.66$\mu$molm⁻²s⁻¹–when a site in the training data has similar ecological characteristics to the unseen sites. This suggests that, with strategic placement of instrumentation to record future data, there is potential to predict most locations of interest with high accuracy. Our research underscores the importance of integrating advanced modeling techniques into carbon accounting frameworks, enabling more accurate quantification of carbon sequestration potential and guiding the implementation of effective natural climate mitigation strategies.

While our results are a significant step forward for quantifying carbon fluxes, we note that this work, like all machine learning-based modeling, is limited by the quality of the training data–the tower-based flux measurements. Both NEON flux measurements, and those from AmeriFlux more generally, have substantial uncertainties [38], but for the most part it is believed by the CO₂ flux measurement community that theoretically-based corrections largely eliminate the systematic biases in measurements, and that the random errors then to average out over long time scales. Unfortunately, validating the accuracy of these measurements is a challenge because of all of the possible pathways by which CO₂ can be removed from, stored, and returned to the atmosphere. Despite the lack of complete standardization across AmeriFlux sites, it is widely believed that tower-based measurements of FCO₂ provide the most accurate and informative estimates of ecosystem carbon uptake and storage. Importantly, interpretation of annual FCO₂ is also possible in a cross-site context, whereas it is not so straightforward to compare biometric forest inventory measurements with estimates of grassland productivity based on biomass clipping, or estimates of agricultural productivity based on crop yield. For this reason, the ability of our ML-based model to successfully predict across-site variation in annual FCO₂ integrals, and within-site inter annual variation in annual FCO₂ integrals, represents an important step forward in coast-to-coast mapping of ecosystem carbon balances, at fine spatial resolution, and the application of these carbon balance estimates in implementing natural climate solutions. Further improvements to our machine learning model (such as domain adaptation techniques [27,28]), using a wider range of data from across the AmeriFlux network, would undoubtedly be of value in advancing this socially-relevant field.