Development of a model framework for terrestrial carbon flux prediction: the Regional Carbon and Climate Analytics Tool (RCCAT) applied to non-tidal wetlands

Wetlands play a pivotal role in carbon sequestration but emit methane (CH₄), creating uncertainty in their net climate impact. Although process-based models offer mechanistic insights into wetland dynamics, they require extensive site-specific parameterisation (e.g. soil carbon profiles, pore-water chemistry, vegetation-specific model parameters), as well as high-resolution hydrological and meteorological inputs that are often difficult to obtain outside of well-instrumented research sites, which makes regional upscaling challenging. In contrast, data-driven models provide a scalable alternative by leveraging available datasets to identify patterns and relationships, making them more adaptable for large-scale applications. However, their performance can vary significantly depending on the quality and representativeness of the data, as well as the model design, which raises questions about their reliability and generalizability in complex wetland systems. To address these issues, we present a data-driven framework for upscaling wetland CO₂ and CH₄ emissions, across a range of machine learning models that vary in complexity, validated against an extensive observational dataset from the Sacramento-San Joaquin Delta. We show that artificial intelligence (AI) approaches, including Random Forests, gradient boosting methods (XGBoost, LightGBM), Support Vector Machines (SVM) and Recurrent Neural Networks (GRU, LSTM), outperform linear regression models, with RNNs standing out, achieving an R² of 0.73 for daily CO₂ flux predictions compared to 0.64 for linear regression, and an R² of 0.53 for CH₄ flux predictions compared to 0.47 for linear regression. Interestingly, linear regression performed better than random forest for methane flux, which highlights the necessity for comparison. The, interannual variability is less well captured, with annual mean absolute error of 176 gC m⁻² yr⁻¹ for CO₂ fluxes and 9 gC-CH₄ m⁻² yr⁻¹ for CH₄ fluxes. By integrating vertically-resolved atmospheric, subsurface, and spectral reflectance information from readily available sources, the model identifies key drivers of wetland CO₂ and CH₄ emissions and enables regional upscaling. These findings demonstrate the potential of AI methods for upscaling, providing practical tools for wetland management and restoration planning to support climate mitigation efforts.

1 Introduction

Wetlands provide a wide array of ecological, economic, and environmental benefits (Costanza et al., 2014). They play a crucial role in biodiversity conservation, water purification, flood control, and climate regulation (Grande et al., 2023; Sharma and Singh, 2021). Significant attention has been recently given to wetland restoration due to their ability to sequester carbon from the atmosphere (Lolu et al., 2020; Upadhyay et al., 2020). These ecosystems are highly effective at storing carbon in their soils because the anaerobic conditions in waterlogged soils suppress organic matter decomposition, allowing carbon to accumulate over time (Mitsch and Gosselink, 2015a). However, wetlands can also be significant sources of CH₄, a potent greenhouse gas (Brix et al., 2001), leading to potentially net positive effects of wetlands on climate warming. The most accurate way to determine the carbon balance in natural ecosystems is through direct and continuous measurements of carbon and GHG sources and sinks (Baldocchi et al., 2001). This involves monitoring carbon dynamics using techniques such as eddy covariance (EC) towers (Aubinet et al., 2012), soil carbon stock assessments (Harrison et al., 2011), and lateral carbon transport measurements (Ciais et al., 2008). However, these measurements are time-consuming to carry out, costly, and require specialized instruments and expertise, limiting their application to a few representative sites globally (Hill et al., 2017; Kumar et al., 2017). The Ameriflux network offers roughly 500 EC sites comprising about 3600 site years of data, monitoring carbon fluxes across various ecosystems such as forests, grasslands, and wetlands (Pastorello et al., 2020). Eddy-covariance site footprints range in scale and are typically determined by the sensor height and atmospheric turbulence (Chu et al., 2021). Data from these Ameriflux sites could potentially be upscaled and used for estimating fluxes from non-monitored sites to obtain regional assessments of carbon balance for various ecosystem types, including wetlands.

In this study, we focus on nontidal wetlands due to the presence of a cluster of EC towers in a small region located in the Sacramento-San Joaquin Delta, including three sites, each with over a decade of continuous data. Reported sequestration rates in wetlands vary widely, influenced by factors such as climate, vegetation, and management. For instance, reported sequestration rates range from as low as 26 gC m⁻² yr⁻¹ in boreal rain-fed bogs (Villa and Bernal, 2018) to as high as 797 gC m⁻² yr⁻¹ in constructed wetlands with emergent Phragmites in the Netherlands (de Klein and van der Werf, 2014). Similarly, temperate wetlands in central Ohio exhibit a wide range of carbon sequestration rates depending on vegetation: forested depressional wetlands dominated by Quercus palustris sequester up to 473 gC m⁻² yr⁻¹, while marshes dominated by Typha sequester around 210 gC m⁻² yr⁻¹ (Bernal and Mitsch, 2012). In Victoria, Australia, freshwater marshes show varying sequestration rates from 91 gC m⁻² yr⁻¹ in shallow marshes to 230 gC m⁻² yr⁻¹ in permanent open freshwater wetlands (Carnell et al., 2018). More relevant to this study, in the San Francisco Bay-Delta region, nontidal managed wetlands dominated by Schoenoplectus and Typha species sequester carbon at rates of approximately 355 ± 249 gC-CO₂ m⁻² yr⁻¹. This estimate is based on direct calculations using Ameriflux data from sites with over a decade of observations (US-Myb, US-Tw1, and US-Tw4). For this calculation we used full-year annual averages and their corresponding standard deviation to the annual mean, to highlight the significant inter-annual variability, with the standard deviation close to the mean. The unit reported for these Delta sites is in gC-CO2 m⁻² yr⁻¹, as the EC tower directly detects CO₂ exchange, which is convenient for GHG assessment purposes. It is worth noting that, at these sites, some years were a net CO₂ source, due to site-specific disturbances such as caterpillar infestations, drought, or when vegetation cover was fully established (Anderson et al., 2018; Knox et al., 2017; Rey-Sanchez et al., 2021). See Table S1 in the Supplement for more detailed information and references therein.

Although CO₂ balance (photosynthesis minus community respiration) is an important component of carbon sequestration, in many wetland systems sequestration benefits are counterbalanced by CH₄ emissions, a potent greenhouse gas, with a warming potential 27 times higher than CO₂ (Lee et al., 2023) that can often offset climate mitigation efforts. CH₄ emission rates also vary substantially over time and across wetlands, from as low as 0.23 gC-CH₄ m⁻² yr⁻¹ in saltwater zones of estuarine environments (Abril and Iversen, 2002) to as high as 270 gC-CH₄ m⁻² yr⁻¹ in certain freshwater wetlands (Knox et al., 2021). For example, restored freshwater wetlands in Maryland dominated by grasses and sedges emit around 142 gC-CH₄ m⁻² yr⁻¹ (Stewart et al., 2024). Tropical wetlands in Costa Rica exhibit some of the highest emissions, with isolated and floodplain wetlands releasing between 220 and 263 gC-CH₄ m⁻² yr⁻¹ (Mitsch et al., 2013). The San Francisco Bay-Delta wetlands that have high carbon sequestration rates also release CH₄ at rates of 35 ± 13 gC-CH₄ m⁻² yr⁻¹(direct measurements from the eddy covariance tower data, Arias-Ortiz et al., 2021). See Table S2 for further information and reference therein. This dual role of wetlands in both sequestering carbon and emitting CH₄ reveals the complex effect they have on the global greenhouse gas balance. Therefore, integrating CO₂ and CH₄ emissions is critical to assess the net climate benefits of wetland conservation and restoration initiatives.

To evaluate how wetlands contribute to the atmospheric radiation budget at larger scales, it is essential to quantify both GHG emissions and carbon sequestration, especially at sites where direct measurements are unavailable (Moomaw et al., 2018). Upscaling models serve this purpose by allowing estimation of sequestration and emission rates across larger spatial scales than those covered by the original data sources (Villa and Bernal, 2018) which provide GHG accounting and net climate benefit assessments for specific wetland sites (Nahlik and Fennessy, 2016). Moreover, it aids in targeting wetland restoration efforts that aim to optimize sequestration by identifying locations with the greatest potential for net carbon uptake.

Process-based models have traditionally been used to estimate sequestration and emissions (Mack et al., 2023; Zhang et al., 2002). Models such as DNDC (Li, 1996), DayCent (Parton et al., 1998), and Ecosys (Grant et al., 2017) have been applied to simulate biogeochemical processes in terrestrial ecosystems, including modeling CH₄ emissions, carbon balances, and soil carbon and nitrogen cycling (Grant and Roulet, 2002; Weiler et al., 2018; Zhang et al., 2002). While these models can elucidate the processes that play a role in carbon dynamics, they require extensive mechanistic parameterization to accurately represent the interactions in various ecosystems (Pastorello et al., 2020; Yin et al., 2023). This approach often necessitates site-specific information and data collection, making implementation over vast areas challenging (Saunois et al., 2025; Xu and Trugman, 2021). The extensive data needs associated with these process-rich models showcase the need for alternative approaches that can effectively upscale wetland emissions without such intensive resource demands.

Artificial Intelligence (AI) methods, such as machine learning and deep learning, have been widely applied in ecological modeling in recent years, alongside long-term, large-scale data collection efforts (Perry et al., 2022). Recent deep learning applications have demonstrated success in capturing the complex dynamics of carbon and methane fluxes in these systems (Ouyang et al., 2023; Yuan et al., 2022, 2024; Zou et al., 2024). The availability of open-source modeling platforms like TensorFlow and PyTorch has made advanced computational techniques, such as neural networks, more accessible, enabling the rapid development and deployment of a range of specialized modeling tasks (Xu et al., 2021). Despite several recent studies demonstrating the potential of machine learning for large-scale carbon cycling in wetland ecosystems, this remains a relatively young field. Moreover, carbon dynamics in wetland ecosystems are temporally variable and inherently nonlinear, making them particularly well-suited for testing machine learning approaches (Arora et al., 2019, 2022). We therefore emphasize the importance of evaluating and comparing various approaches within this domain and their potential for large-scale assessment.

A pervasive challenge in model development is the ability to balance complexity with generalizability. While more complex models can capture nonlinear relationships, they also increase the risk of overfitting, where the model performs well in the testing, but poorly on new conditions (Hastie, 2009; Tashman, 2000). Furthermore, it is also important to use a robust validation framework. For the application of upscaling, it is important that the model is able to extrapolate spatially. For this purpose, a leave-one-site-out (LOSO) validation approach is typically carried out, whereby the models are trained on data that excludes a single site, with the excluded site data saved for model testing (Bodesheim et al., 2018; Tramontana et al., 2016). It is also important to avoid data leakage, where information from the training set inadvertently appears in the testing set (Kaufman et al., 2012), a risk posed when splitting temporally adjacent data points that are close in value, potentially inflating performance statistics (Bergmeir and Benítez, 2012; Kaufman et al., 2012). For example, daily rates of change relative to a system where seasonal dynamics dominate, such as emissions of CH₄ emissions in vegetated wetlands (Knox et al., 2021).

In this study, we introduce a model framework for coastal nontidal wetland CO₂ and CH₄ emissions using several “off-the-shelf” models. These models are trained and validated against observational data, and results are compared to find the most predictive model. The top performing model is then used to upscale carbon sequestration and CH₄ emissions in nontidal wetlands at regional scale. The San Francisco Bay-Delta serves as the area of interest, due to its network of EC towers that have been operating for a relatively long time and relevance to future wetland restoration efforts. We employ a suite of models, ranging widely in complexity: (1) linear regression; (2) Random Forests (Breiman, 2001), an ensemble method that constructs multiple decision trees to reduce overfitting; (3) gradient boosting techniques such as LightGBM (Ke et al., 2017) and XGBoost (Chen and Guestrin, 2016), which are scalable tree boosting systems able to handle complex nonlinear relationships and variable interactions; (4) Support Vector Machines (SVM) (Cortes, 1995), a kernel-based technique that can approximate nonlinear boundaries between data points and (5) the Recurrent Neural Network (RNN) such as the Long Short-Term Memory (LSTM) neural network (Hochreiter, 1997), an advanced model designed to process sequential data and capture non-linear interactions over long-term dependencies. We also test a model with similar but simpler architecture, the Gated Recurrent Unit (GRU) (Chung et al., 2014), which uses fewer parameters. Linear regressions serve as a baseline to assess the applicability of the more sophisticated methods. Random Forests have been used to upscale northern wetland methane emissions (Peltola et al., 2019), gradient boosting methods have demonstrated success in ecological modeling (Ding, 2024; Räsänen et al., 2021; Zou et al., 2024), and LSTM neural networks have been successfully applied to model CO₂ and CH₄ fluxes in ecosystems (Yuan et al., 2022, 2024; Zou et al., 2024). Our proposed framework is designed to provide transparency, easy determination of model practicality and applicability, and contextualisation to model performances by comparing to a baseline model (i.e. linear regression).

2 Methods

Our ultimate aim is to establish a robust modeling framework for estimating wetland carbon fluxes in sites that are not monitored. To achieve this, we compare a range of models, from simple linear regression to advanced recurrent machine learning neural networks. Since the goal is to predict unseen sites, we emphasize cross-site predictability by validating and testing the models at sites not included in training. Doing so ensures predictions are applicable beyond the training sites and addresses challenges often associated with model generalizability (Meyer and Pebesma, 2022). This strategy serves several purposes:

1. Performance Contextualization. Starting with the simplest type of model provides a baseline for performance and helps evaluate the advantage (or lack thereof) for using more complex models.

2. Practicality and Transparency. Advanced models may offer better performance but often require significant effort to set up and may lack interpretability. By comparing models of varying complexity using the same input data, we assess whether the added complexity is justified.

3. Feature Evaluation. Training with different combinations of relevant features helps us to understand which features are dominating control, and the limitations of the data in terms of predictive capacity.

2.1 Model targets

The model targets two key variables: CO₂ (FCO2) and CH₄ (FCH4) surface emissions. Both variables follow a sign convention where positive values indicate emissions to the atmosphere (source) and negative values indicate sequestration (sink). Both variables are available at half-hourly resolution through the Ameriflux database.

The models we developed all operate on a daily time scale, requiring target variables to be aggregated to the daily time scale. This approach assumes that sub-daily variations have a negligible non-linear contribution to longer time scales, an assumption supported by the dominant seasonal signal typically observed in flux data from these systems (Knox et al., 2021).

These target variables could then be used to calculate annual NECB (Net Ecosystem Carbon Balance; gC m⁻² yr⁻¹) and annual wetland net atmospheric radiative effect (FCO2e (CO₂-equivalent flux) gCO₂e m⁻² yr⁻¹).The global warming potential (GWP) of non-fossil CH₄ is 27.2 as per the latest IPCC assessment (Lee et al., 2023). For this study, we neglect contributions of lateral fluxes due to data limitations, and that lateral transport at these sites is assumed to be negligible due to the limited outflow from these specific non-tidal wetland sites (Miller et al., 2008). FCO2e is defined as annually averaged CO₂ and CH₄ emissions, adjusted for the global warming potential (GWP) of each gas. A positive FCO2e indicates that the ecosystem is contributing positively to atmospheric warming, and vice versa. Here we consider CO₂ and CH₄ emissions but neglect contributions from N₂O due to data limitations and because N₂O emissions are considered negligible in Delta wetlands (Windham-Myers et al., 2018).

2.2 Region of interest

The Sacramento-San Joaquin Delta was selected for this study due to its high density of EC towers and extensive long-term data. We selected sites for model training and validation where data was collected for at least a decade to capture interannual variability. Hence three restored wetland sites, US-Myb (Matthes et al., 2025), US-Tw1 (Valach et al., 2024), and US-Tw4 (Eichelmann et al., 2025) are selected in this study. While data from two other sites (i.e., US-Sne and US-Tw5) are available, the lack of sufficient temporal coverage and, in the case of US-Sne, not fully established vegetation cover, makes them less representative of a stable ecosystem. Focusing on sites with over a decade of continuous data allows for capturing long-term dynamics more effectively and provides sufficient time for the wetlands to reach a stable state. The dataset encompasses 35 full site-years of observations across the three sites within the Delta (Novick et al., 2018) (Table 2, Fig. 1), with detailed mapping data sourced from the Ecoatlas Database (Workgroup, 2019) which provides land use and vegetation surveys across wetlands in California.

Table 1 Model training sites.

Table 2 An overview of the models that are applied to wetland fluxes.

The sites are dominated by Tules (Schoenoplectus), Cattails (Typha), and invasive species such as Phragmites, which are perennial emergent plants well suited to wetland environments (López et al., 2016). The Delta itself is host to the largest estuarine system on the US Pacific coast, spanning approximately 3000 km², and contains a diverse network of wetland systems. Historically, much of the area was drained and converted for agriculture (Laćan and Resh, 2016; Lund et al., 2010), but recent restoration efforts have reclaimed select portions of the landscape for environmental benefits.

Figure 1 Map of the Sacramento-San Joaquin Delta's wetland system. The Eddy-covariance tower site locations outlined in Table 2 are shown in the red and purple boxes. Satellite image: © Google Earth, accessed 2025.

2.3 Model features

The application of this work focuses on upscaling carbon fluxes from similar wetlands at a regional scale. To achieve this, we aim to predict fluxes at unmonitored sites using widely available data that are expected to be key drivers of FCO2 and FCH4. Since site-level measurements from EC towers are not available at a larger spatial scale, we focus on ecosystem drivers that can be accessed across broader spatial extents.

The models utilize a comprehensive set of features from two readily accessible datasets: (i) the Western Land Data Assimilation System (WLDAS) (Erlingis et al., 2021) and (ii) Landsat surface-reflectance products (Landsat, 2020). A list of features can be found in Table S3. Initially surface reflectance products were derived from MODIS (Justice et al., 2002), but we found better model performance with Landsat features. WLDAS provides hydrological and meteorological data at 1 km spatial and daily temporal resolution; we bilinearly interpolate these fields to each tower coordinate (no additional smoothing). Landsat offers 30 m pixels at a nominal 16 d revisit, although temporal resolution increases with time as more satellites are added; we average a 3 × 3 pixel window centred on the tower, linearly interpolate the series to daily resolution, and apply a centred 17 d running mean to improve data continuity.

2.4 Model suite

To evaluate ML model performance in calculating FCO2 and FCH4, we implemented a suite of seven models ranging from simple linear methods to more complex neural networks. These models have been used in various ecosystems to study fluxes and collectively represent a broad spectrum of methodological complexity. Table 3 summarizes the core characteristics and advantages of each approach.

Table 3 Feed-forward feature selection process.

These models act to demonstrate a spectrum of model complexity and how that can be leveraged to improve flux prediction.

After performing simple grid searches we found that all models were largely insensitive to hyperparameter tuning, so we kept almost everything at the package defaults with some minor exceptions. Model hyper-parameter choices can be found in Brereton (2025).

2.5 Validation framework

To evaluate the models' ability to generalize across sites, we employed a Leave-One-Site-Out (LOSO) cross-validation strategy. In LOSO, we train the models on data from all but one site, and test the models on the excluded site. This approach is repeated for each site in the dataset and then aggregated, ensuring that there are no spatio-temporal connections between the training and testing data. While few models are immune to overfitting, this approach minimizes the risk of doing so.

An integral part of our modeling approach is the strategic selection of input features to optimize the model's performance. We perform this selection by first selecting features that are expected to be important, guided by mechanistic considerations of wetland processes gained from fieldwork and insights from mechanistic models (Table S3). Since the total number of possible feature combinations is too large for an exhaustive search, we adopt a feed-forward selection (FFS) strategy. This method begins with a single feature and iteratively adds features that most improves the model's performance based on a chosen statistic. At each step, we evaluate the model's performance with each potential new feature and select the one that provides the greatest improvement. This process continues until adding additional features no longer significantly enhances the model's performance. By using this approach, we efficiently identify the most influential predictors without the computational burden of testing all possible combinations.

2.6 Validation

As suggested above, each model was trained using data from two wetland sites and then validated on the third. Although the number of sites was limited, each site offered over a decade of observations accumulated to a daily time step, ensuring exposure to a range of environmental conditions representative of the wetland type and regional climate. For each excluded site, the model's predictions were compared against measured FCO2 and FCH4 and we calculated R² , Pearson's r, and RMSE for that site. We then pooled all held-out predictions from the three sites into one combined set and recomputed R² (as well as r and RMSE) on the full array to give an overall cross-validation score. This process was paired with the FFS method optimized to maximize R².

After selecting LSTM as the model of choice, it was retrained using all available data from the three sites for upscaling. The Sacramento-San Joaquin Delta contains roughly 700 km² of wetland area, including tidal and nontidal regions. The upscaling domain encompasses approximately 25 km² of nontidal wetlands in the region, dominated by vegetation types relevant to the training sites, specifically Tules, Cattails, and Phragmites. The assumption is that the training sites used in this study are representative of the broader conditions in the Delta, but we acknowledge that local variability in carbon dynamics, such as those caused by microclimates prevalent in the area, may not be fully captured during the ML model training. Improvements to the model might be achieved if additional site data covering a wider range of environmental conditions were incorporated. The feature data used to optimize the model were spatially interpolated onto the regional model grid and the model applied to yield flux estimations. Although relatively modest in spatial extent, these wetlands are of particular interest given their role in carbon sequestration and potential climate mitigation and as targets for conservation and restoration.

3 Results

3.1 Model Validation

We tested six modeling techniques of varying complexities (Table 3). Model performance scores for daily predictions are shown in Fig. 2, demonstrating that nearly all machine learning models outperformed the linear regression baseline (R² = 0.64 for FCO2 and R² = 0.47 for FCH4). For FCO2, LSTM and GRU achieved the highest R² values (0.73 and 0.71, respectively), outperforming other methods. A similar result was found for FCH4, with LSTM and GRU both scoring R² of 0.53. These results suggest that deep learning models can provide tangible benefits over linear regression methods for upscaling flux predictions. The LSTM model was selected for upscaling in this study as it scored highest consistently, though other ML models scored comparably, so we do not assert it as definitively the best model.

Figure 2 Bar plot showing best model performance for each type of machine learning model based on R² score (though other metrics are in agreement, see Pearson r correlation and RMSE).

The feature selection process had access to 26 environmental features from WLDAS and 7 features derived from LANDSAT spectral bands (see Table S3 for full details). These variables encompass a wide range of atmospheric, soil, and vegetation characteristics, such as precipitation, temperature, soil moisture, and spectral indices, key environmental drivers known to influence carbon and methane flux dynamics (Mitsch and Gosselink, 2015).

The feature selection routine converged on variables that map directly onto the three main controls of wetland carbon cycling – vegetation productivity, surface energy-water balance, and microbial temperature sensitivity, see Table 4. For FCO2, the features selected were the Soil-Adjusted Vegetation Index (SAVI) and the upwards sensible heat flux, which are proxies for gross primary production and the surface energy water balance (Anderson et al., 2016; Huete, 1988). For FCH4, the features selected were canopy temperature, soil temperature and Greenness Difference Vegetation Index (GNDVI), which are proxies for short-term thermal forcing and vegetation water status, the anaerobic root-zone temperature that governs methanogenesis, and the supply of photosynthetically derived substrates for microbes, respectively (Bubier et al., 1993; Knox et al., 2021; Whiting and Chanton, 1993; Yvon-Durocher et al., 2014).

Figure 3 shows both FCO2 and FCH4 results, including time series and scatter plots comparing predictions to observations. Overall, the predicted values track the observations reasonably well. For FCO2, predictions tended to regress toward the mean, underestimating peak emissions at local maxima and overestimating at local minima, although reasonable interannual variability was observed. The ML models also displayed less interannual variability than the observations, common in machine learning approaches (Ouyang et al., 2023). For wetlands, this is likely due to limited subsurface process information included in the machine learning models. Still, the scatter plot shows strong performance for FCO2 (r = 0.86, R² = 0.73, RMSE = 1.46 gC-CO₂ m⁻² d⁻¹), despite a noticeable spread around the 1:1 line.

FCH4 predictions exhibited similar behavior, with low interannual variability than the observations. At the US-Myb site, for example, observed FCH4 were initially high (aside from the first year, when vegetation cover had yet to be fully established) but declined over time, stabilizing at lower values. The ML models captured this shift to some extent, predicting higher fluxes early in the time series and then modulating to lower levels later on. However, predictions did not fully replicate the magnitude of the observed downward annual trend, introducing bias into the scatter plots at higher and lower extreme values. This phenomenon is known as regression to the mean, observed in similar machine learning studies (Ouyang et al., 2023). Consequently, the FCH4 model performance was weaker than the FCO2 model (R² = 0.53, r = 0.74, RMSE = 0.05 g C-CH₄ m⁻² d⁻¹), indicating that the processes controlling FCH4 in younger wetlands like US-Myb may require more detailed subsurface information (such as soil organic C, oxygen, or redox information) to be accurately modeled. Restored Delta wetlands are often net GHG sources for 1–3 years after flooding, before vegetation is fully established. Eddy-covariance measurements show positive NEE of +201 ± 101 g C-CO₂ m⁻² yr⁻¹ and elevated CH₄ emissions in the initial period, switching to sinks of between −400 to −700 g C-CO₂ m⁻² yr⁻¹ thereafter (Hemes et al., 2019). A larger synthesis found that this can persist decades in nontidal marshes because CH₄ radiative forcing outweighs CO₂ burial (Arias-Ortiz et al., 2021). Similar contrasts between 2 and 15-year-old wetlands (Knox et al., 2015).

The annual bar plots presented in Fig. 4 highlight the model's difficulty in capturing the interannual variability of carbon fluxes across the study sites. While the average FCO2 and FCH4 predictions are generally aligned with observed average values with small overall mean bias, the model struggles to reproduce the observed year-to-year variability. Although direct subsurface measurements are available at certain sites, at the regional scale their limited spatial and temporal coverage currently limits integration into models designed for regional upscaling over inter-annual timescale. For example, while spatial maps of wetland soil organic carbon exist (Uhran et al., 2022), using only three sites for training purposes would provide just three corresponding data points, limiting model training. The LOSO validation approach revealed that deep learning models, particularly LSTM and GRU, consistently outperformed traditional linear regression and other machine learning methods for both FCO2 and FCH4 predictions. While nonlinear models demonstrated clear advantages, the magnitude of improvement was relatively modest, reflecting the inherent challenges of capturing site-specific inter-annual dynamics of wetland emissions. To improve model performance, additional techniques such as feature transformations or attention mechanisms could be implemented. However, the primary goal of this model suite is to ensure reproducible results with “off-the-shelf” models, which serves as a foundation for more advanced, nuanced approaches.

Figure 3 Time-series plots (left) of observed (blue) and predicted (orange) FCO2 and FCH4 fluxes for US-Myb, US-Tw1, and US-Tw4. The scatter plot (right) compares observed vs. predicted values across all sites, with a 1:1 reference line with overall and site-only performance metrics (R² ).

Figure 4 Annual Observed and Predicted FCO2 and FCH4 Across Three Wetland Sites. Aggregated statistics for all sites are as follows: For FCO2, the Mean Absolute Error (MAE) is 176 gC m⁻² yr⁻¹ and the Mean Bias Error (MBE) is 17 gC m⁻² yr⁻¹. For FCH4, MAE is 9 gC-CH₄ m⁻² yr⁻¹ and the MBE is 1 gC-CH₄ m⁻² yr⁻¹.

3.2 Model Application: Upscaling

Upscaling was repeated 10 times for 10 separately trained ML models using the same data, and the ensemble mean is the value reported, as training includes stochasticity. Figure 5 displays spatial maps of annual flux estimates of Net Ecosystem Carbon Balance (NECB), and the CO₂ equivalent flux rate (FCO2e) in the study domain, including zoom-in subplots highlighting areas with more data. The results show that carbon sequestration, indicated by negative NECB (green) values, are typically dominant throughout the domain, although the northern regions show more carbon sources. In contrast, the FCO2e distribution shows variability across the region, with sources and sinks found CO₂e sink throughout .

Figure 5 Mean annual Net Ecosystem Carbon Balance (NECB, left) and CO₂-equivalent radiative forcing (FCO2e, right) averaged over all model years. Main maps show the Delta area; dashed rectangles (1–3) correspond to zoom-in panels. Tidal wetlands are shaded dark blue, non-tidal light blue. Positive values (red) indicate net carbon loss; negative values (green) net uptake. See Fig. 1 for reference to training sites.

Figure S1 plots the coefficient of variation (CV $=\mathit{\sigma }/\mathit{\mu }$) of the inter-model ensemble for both NECB and FCO2e. Higher CV indicates locations where environmental conditions are poorly represented in the training data – effectively a proxy to determine model confidence. Across the study domain the vast majority of pixels show low dispersion: ≈ 85 % of the mapped area has a CV < 0.5 for NECB, and 69 % falls below that same threshold for FCO2e.

Figure 6 shows averaged fluxes in the upscaling domain over the full study period. The results highlight the Delta as an overall carbon sink, with NECB averaging approximately −450 gC m⁻² yr⁻¹, indicating persistent sequestration across multiple years. CH₄ fluxes average 31 gC-CH₄ m⁻² yr⁻¹, and shows little spatial variability. Values are consistent with those previously reported in the region (Arias-Ortiz et al., 2021). Integrating these fluxes into a CO₂-equivalent metric, this regional wetland system remains a net sink of CO₂e, with approximately 600 gCO₂e m⁻² yr⁻¹ sequestered on average in the upscaling domain, with an increasing trend with time.

Figure 6 Bar plots and box plots of annual NECB, FCH4, and FCO2e fluxes, which have been spatially integrated over the study region, a total of 25 km² total land area vegetated primarily by Tules, but also Cattails and Phragmites. The left column shows annual fluxes for each year, with negative fluxes in green and positive fluxes in orange. Daily fluxes, aggregated to annual totals, are overlaid as grey lines. The right column shows box plots summarizing the distribution of annual fluxes, highlighting the range, median (blue line), and spread of values. Each row represents a different flux variable: (a) NECB, (b) FCH4, and (c) FCO2e.

4 Discussion

This study demonstrates the development and evaluation of a data-driven framework to upscale terrestrial CO₂ and CH₄ flux estimates for non-tidal wetlands in the Sacramento-San Joaquin Delta. By systematically comparing models of varying complexity, including linear regression, ensemble methods, gradient boosting algorithms, and recurrent neural networks (RNNs), we presented a transparent assessment of model performance. The goals were to identify the model that best predicts CO₂ and CH₄ fluxes and critically appraise whether incremental complexity is justified by improvements in predictive capacity. Relevant cited works have included many different machine learning approaches for predicting emissions. This work aims to unify modelling efforts by establishing a standard framework for developing robust data-driven models, particularly for upscaling purposes.

Our results indicate that non-linear and more advanced models generally outperformed simple linear regression approaches. Among all tested models, the Long Short-Term Memory (LSTM) and Gated Recurrent Unit (GRU) neural networks provided the highest overall skill in predicting both CO₂ and CH₄ fluxes at daily timescales. This improvement was marginal but consistent, supporting the notion that time-series models, which inherently capture temporal dependencies and non-linearities, can provide tangible benefits over linear methods and traditional machine learning algorithms.

However, while these deep learning models performed best, the performance gains were not as large as might be expected given their significantly higher complexity and computational demands. Similar outcomes have been noted in other ecological modeling applications, where advanced machine learning methods yield improvements that are statistically significant yet modest in terms of performance gains relative to linear models (Oh et al., 2022; Wood, 2022).

The deep learning models provided reasonable estimates of daily fluxes but struggled to replicate the full range of interannual variability observed in the field measurements, which is a common issue for data-driven models in this field (Nelson et al., 2024). This limited ability to capture long-term trends and extremes mirrors common challenges in machine learning-based modeling, where the absence of explicit mechanistic understanding limits extrapolation beyond the conditions represented in the training data. The difficulty in reproducing interannual fluctuations was particularly evident for CH₄ fluxes, an outcome consistent with the high spatial and temporal complexity of CH₄ cycling in wetland environments and the limited availability of subsurface parameters (e.g., oxygen concentration, redox conditions, substrate availability) that drive CH₄ production. This may not be surprising as the number of annual cycles available in the training set was only 35 years.

The observed regression to the mean and the reduced dynamic range in model predictions may reflect insufficient representation of key environmental drivers in the feature set or inadequate temporal coverage and variability in the training data. While publicly available datasets such as WLDAS and LANDSAT were effective at providing spatially and temporally comprehensive inputs, the lack of direct subsurface and soil biogeochemical measurements likely limited the model's ability to capture critical internal processes that are likely causing the observed differences between years. Although the feed-forward selection process for the model features had access to an extensive pool of relevant features, results indicated that only a small subset of features was necessary to maximise performance. This suggests that, while there are many features that control CO₂ and CH₄, their contribution to predictive accuracy may be redundant or captured indirectly by other variables. The exclusion of particular features, such as the water table depth for FCH4, illustrates the trade-off between mechanistic intuition and data-driven optimization. Strong correlations between features with limited independent variability can lead to features being left out that would typically be considered ecologically relevant.

The feed-forward selection converges on a compact set of features that are mechanistically expected to drive CO₂ and CH₄ gas exchange in natural wetlands, adding qualitative confidence to the model; see Table 4 for the list of features selected. For FCO2; the model selected SAVI (Soil-Adjusted Vegetation Index), which represents vegetation state and photosynthetic capacity, and upward sensible heat flux (Q_h), which is related to the rate of evapotranspiration. For FCH4, canopy temperature and GNDVI (Greenness index) indicate plant activity and substrate supply via plant-mediated transport, while 10–40 cm soil temperature reflects anaerobic microbial production (methanogenesis) in the root zone.

After applying the chosen model (LSTM) to calculate CO₂ and CH₄ fluxes, we estimated NECB and CO2-equivalent fluxes for similar wetland settings across the Delta region (Fig. 5). The results show spatial heterogeneity and pinpoint regions that act as stronger net carbon sinks, as well as areas where CH₄ emissions may offset climate benefits of net carbon sequestration. Such insights support targeted conservation and restoration strategies aimed at maximizing net carbon sequestration benefits, facilitating ongoing efforts to restore and manage wetlands to contribute to net-zero emission goals.

A key advantage of the chosen approach is its reliance on readily available, open-source data streams and standard computational resources. The framework can be deployed efficiently without specialized hardware, making it accessible to resource-limited organizations, practitioners, and researchers.

The primary objectives of this study were to identify a suitable model, contextualize model performance by comparing to a baseline linear regression, and highlight trade-offs between complexity, interpretability, and accuracy. By explicitly testing multiple models ranging from simple linear regressions to advanced recurrent neural networks, we demonstrated that complexity alone does not guarantee a substantial increase in predictive power. Instead, complexity should be adopted judiciously, based on the magnitude of performance gains, the cost of model implementation, and the level of interpretability.

We suggest that future modeling efforts should focus on deriving mechanistically relevant predictors (Ouyang et al., 2023), and incorporating hybrid modeling approaches (Yao et al., 2023) that combine the strengths of process-based and machine learning methods. Attention mechanisms (Yuan et al., 2022), advanced architectures (e.g., Transformers, Vaswani, 2017), or physics-informed machine learning (Raissi et al., 2019) may also help address model performance limitations.

5 Conclusions

This study provides a transparent, methodical demonstration of an artificial intelligence approach to modeling wetland carbon dioxide (CO₂) and methane (CH₄) emissions, using a suite of “off-the-shelf” tools and establishing a standardized benchmarking protocol for model performance evaluation. In the study region (the Sacramento–San Joaquin Delta), inter-model comparisons revealed modest but appreciable performance differences when comparing advanced models with a linear regression baseline. While there are tangible benefits to employing machine learning for these purposes, it is likely that the gap between simpler models and more sophisticated models will widen as data quantity and quality continues to increase. Ultimately, this study lays the groundwork for regional scale model benchmark testing, facilitating the development of more advanced modeling approaches that can guide wetland management, restoration planning, and climate mitigation strategies.

Appendix A

Figure A1 Coefficient of variation ($\mathit{\sigma }/\left|\mathit{\mu }\right|$) of NECB (left) and FCO2e (right) across the 10 runs, showing relative inter-model uncertainty in the predictions. Green areas show high model confidence and red areas shower either lower model confidence (or division by small μ). This can be interpreted as a proxy for the confidence of the spatial upscaling.