WetCH₄: a machine-learning-based upscaling of methane fluxes of northern wetlands during 2016–2022

Wetlands are the largest natural source of methane (CH₄) emissions globally. Northern wetlands (>45° N), accounting for 42 % of global wetland area, are increasingly vulnerable to carbon loss, especially as CH₄ emissions may accelerate under intensified high-latitude warming. However, the magnitude and spatial patterns of high-latitude CH₄ emissions remain relatively uncertain. Here, we present estimates of daily CH₄ fluxes obtained using a new machine learning-based wetland CH₄ upscaling framework (WetCH₄) that combines the most complete database of eddy-covariance (EC) observations available to date with satellite remote-sensing-informed observations of environmental conditions at 10 km resolution. The most important predictor variables included near-surface soil temperatures (top 40 cm), vegetation spectral reflectance, and soil moisture. Our results, modeled from 138 site years across 26 sites, had relatively strong predictive skill, with a mean R² of 0.51 and 0.70 and a mean absolute error (MAE) of 30 and 27 nmol m⁻² s⁻¹ for daily and monthly fluxes, respectively. Based on the model results, we estimated an annual average of 22.8±2.4 Tg CH₄ yr⁻¹ for the northern wetland region (2016–2022), and total budgets ranged from 15.7 to 51.6 Tg CH₄ yr⁻¹, depending on wetland map extents. Although 88 % of the estimated CH₄ budget occurred during the May–October period, a considerable amount (2.6±0.3 Tg CH₄) occurred during winter. Regionally, the Western Siberian wetlands accounted for a majority (51 %) of the interannual variation in domain CH₄ emissions. Overall, our results provide valuable new high-spatiotemporal-resolution information on the wetland emissions in the high-latitude carbon cycle. However, many key uncertainties remain, including those driven by wetland extent maps and soil moisture products and the incomplete spatial and temporal representativeness in the existing CH₄ flux database; e.g., only 23 % of the sites operate outside of summer months, and flux towers do not exist or are greatly limited in many wetland regions. These uncertainties will need to be addressed by the science community to remove the bottlenecks currently limiting progress in CH₄ detection and monitoring. The dataset can be found at https://doi.org/10.5281/zenodo.10802153 (Ying et al., 2024).

1 Introduction

Methane (CH₄) is the second most important greenhouse gas after carbon dioxide (CO₂) and has contributed to around one-third of anthropogenic warming (IPCC, 2024). Wetlands are the largest natural source of CH₄ emissions. Northern freshwater wetlands (>45° N) account for roughly 40 % of global wetland area (ranging from 1.3×10⁶ to 8.7×10⁶ km²; Zhang et al., 2021), yet the amount of CH₄ emissions from this region is highly uncertain – currently estimated to be 22–49.5 Tg CH₄ yr⁻¹ (Aydin et al., 2011; Baray et al., 2021; Heimann, 2011; Kirschke et al., 2013; Peltola et al., 2019; Saunois et al., 2020; Treat et al., 2018; Watts et al., 2023). The uncertainties in the estimates of wetland CH₄ emissions are primarily attributed to challenges in mapping vegetated wetlands versus open water, leading to double-counting (Thornton et al., 2016); seasonal wetland dynamics; and uncertainties in estimates on flux rates.

Characterized by nutrient, moisture, and hydrodynamic conditions, northern freshwater wetlands are classified as wet tundra in treeless permafrost areas, peat-forming bogs and fens in boreal and temperate biomes, and permafrost bogs (Kuhn et al., 2021; Olefeldt et al., 2021). Bogs were estimated to cover the largest area (1.38×10⁶–2.41×10⁶ km²) in the northern high latitudes, followed by fens (0.76×10⁶–1.14×10⁶ km²) and wet tundra (0.31×10⁶–0.53×10⁶ km²) (Olefeldt et al., 2021). Climate change poses significant threats to these wetlands, affecting their extent and the duration of conditions suitable for wetland formation in permafrost zones (Avis et al., 2011). The rates of CH₄ emissions may increase quickly because of intensified warming at the northern high latitudes (Masson-Delmotte et al., 2021; Rawlins et al., 2010; Rößger et al., 2022; Walsh, 2014; Zhang et al., 2023).

Reflecting the CH₄ response to warming, northern wetlands may account for a high portion (∼78.5 %) of the global surface emission anomaly of CH₄ in 2020 relative to 2019 (6.0±2.3 Tg CH₄ yr⁻¹) (Peng et al., 2022; Zhang et al., 2023). This is concerning, as the responses of high-latitude CH₄ emissions to a warming and possibly wetting climate could produce a positive carbon–climate feedback (McGuire et al., 2009; Natali et al., 2019). However, the ability of models to account for and predict the spatiotemporal variability in high-latitude wetland CH₄ emission rates remains very limited (Treat et al., 2024).

Field observations of gas fluxes typically measure CH₄ exchange between the land and atmosphere at sub-meter to ecosystem (hundreds of meters to kilometers) scales (Bansal et al., 2023; Chu et al., 2021). Tower eddy-covariance (EC) methods provide near-continuous measurements over ecosystem-scale footprints (5×10³–100×10³ m²), the size of which matches fine- to medium-resolution satellite remote sensing. Typical EC measurement system records include carbon, water, and energy fluxes along with environmental conditions at half-hourly intervals. Long-term EC datasets can support the analysis of daily, monthly, seasonal, or interannual patterns and drivers of carbon fluxes (Baldocchi, 2003). Chambers can also measure CH₄ fluxes, although at sub-meter resolution and with low spatial coverage (Bansal et al., 2023; Kuhn et al., 2021). Most chamber studies have a limited temporal sampling period. To avoid footprint disagreement between EC and chamber measurement techniques, we focused on EC-based CH₄ upscaling in this study.

Data-driven upscaling uses empirical models (Bodesheim et al., 2018; Jung et al., 2011), including machine learning (ML) approaches, to compute CH₄ fluxes. It provides independent estimates to those from process-based models and atmospheric inversions (Bergamaschi et al., 2013; Spahni et al., 2011). These approaches have been used to estimate CH₄ fluxes from various ecosystems, such as northern wetlands (Peltola et al., 2019; Virkkala et al., 2024; Yuan et al., 2024), global reservoirs (Johnson et al., 2021), and global aquatic ecosystems (Rosentreter et al., 2021).

Two types of methods are generally used for data-driven upscaling. The first uses a look-up table approach and applies emission rates or emission factors via data synthesis to the corresponding land surface areas, or activity data, over the study region. Emission rates from field observations are associated with environmental drivers that have been spatially characterized and are then applied to the land covers with the same environmental drivers. For example, Rosentreter et al. (2021) collected 2601 CH₄ flux records measured using various methods through a literature review and characterized emission rates over 15 aquatic ecosystem types to upscale global aquatic CH₄ emissions. The study provided estimates of total and per-ecosystem emissions but did not produce spatial distributions and was unable to estimate temporal changes. A similar method was applied for the northern permafrost region, where statistical CH₄ flux rates from the Boreal–Arctic Wetland and Lake CH₄ Dataset (BAWLD-CH₄) were analyzed for emission estimation by wetland type (Kuhn et al., 2021; Ramage et al., 2024). This method favors homogeneous ecosystems and static environments, and the results may be biased for large-scale studies where spatial heterogeneity is prevalent.

Another approach uses ML methods to upscale fluxes (Bodesheim et al., 2018; Tramontana et al., 2016; Yuan et al., 2024). ML models are developed with large training datasets. Generally, ML models can learn from high-dimensional data by optimizing many statistical parameters and identifying variables that are closely associated with spatiotemporally varying CH₄ emissions. The efficient computational cost makes it easier to apply the models over large regions at higher spatial resolutions. Among ML methods, decision-tree-based algorithms have been widely used in upscaling for computational efficiency and predictive performance (Beaulieu et al., 2020; Jung et al., 2020; Virkkala et al., 2021; Zhang et al., 2020). Specifically, random forest (RF) algorithms have been utilized in regional to global wetland CH₄ upscaling (Davidson et al., 2017; Feron et al., 2024; McNicol et al., 2023; Peltola et al., 2019) due to their robustness and the prevention of overfitting to noise in the input data. For example, Peltola et al. (2019) used RF and EC measurements to upscale monthly CH₄ fluxes from the northern wetlands at 0.25–0.5° spatial resolution over the 2013–2014 period.

ML-based upscaling studies usually incorporate information from remote sensing to inform wetland extent, changes in vegetation, and other surface biophysical properties (Davidson et al., 2017; Virkkala et al., 2024; Watts et al., 2014, 2023). For example, recent ML-based large-scale upscaling approaches used the Moderate Resolution Imaging Spectroradiometer (MODIS) land surface temperature (LST) at night, the enhanced vegetation index (EVI), the vegetation canopy height, and ancillary environmental variables from remote-sensing products (McNicol et al., 2023; Ouyang et al., 2023; Peltola et al., 2019). (See Sect. S1 and Table S1 in the Supplement for the detailed predictive variables used in existing ML-based wetland CH₄ upscaling products.) However, soil moisture and soil temperature, two controlling factors of wetland CH₄ fluxes (Knox et al., 2021; Yuan et al., 2022), were missing in previous ML-based regional to global upscaling studies. Soil moisture has been identified as one of the important controlling factors for freshwater wetland CH₄ fluxes (Euskirchen et al., 2024; Voigt et al., 2023). This is the first ML-based study that incorporates remote-sensing constraints from Soil Moisture Active Passive (SMAP) microwave-sensed soil moisture and MODIS Nadir Bidirectional Reflectance Distribution Function (BRDF)-Adjusted Reflectance (MODIS NBAR) data. Surface reflectance provides information on (1) vegetation properties that affect the production and transport of CH₄ to the atmosphere and (2) ecosystem wetness (Alonso et al., 2020; Chen et al., 2013; Entekhabi et al., 2010; Houborg et al., 2007; Murray-Hudson et al., 2015; Wang et al., 2018).

The goal of this study is to develop a scalable framework to upscale daily CH₄ fluxes from EC observations to northern-latitude wetlands (>45° N) using the ensemble RF ML approach with a suite of reanalysis and remote-sensing products representing spatiotemporal environmental conditions. Our specific objectives are as follows:

1. compile an updated EC-based CH₄ flux dataset that extends the temporal and spatial coverage of the FLUXNET-CH₄ database (Delwiche et al., 2021) for the northern high latitudes;

2. build RF models of CH₄ fluxes at the site level based on physical variables that are measured in situ, thereby allowing us to prioritize the selection of gridded variables for upscaling, and then build ensemble RF models at the grid level using gridded reanalysis inputs and constraints from satellite remote sensing; and

3. apply grid-level models to produce a 10 km gridded daily distribution of the CH₄ flux product for the northern high-latitude wetlands using bootstrapped models and their derived uncertainties (Table S1).

Figure 1 Workflow and experimental design: abstract methodological steps are integrated in the dashed box on the left, while a detailed experimental design is described on the right. Colors on the right match the associated step on the left.

2 Materials and methods

2.1 Overview

The scalable framework of upscaling CH₄ fluxes from EC observations for wetlands (hereafter referred to as WetCH₄), which selects physical predictors at the site level and constructs upscaling models at a grid level, is illustrated in Fig. 1. In situ, reanalysis, and remote-sensing products were compiled as candidate predictors for modeling (Fig. 1, purple boxes; see Sect. 2.2 for details). We first ran a feature selection, which uses ensemble RF models to choose important predictors from an extensive list of in situ variables available from the flux tower sites. Gridded versions of selected site variables were taken from version 2 of the Modern-Era Retrospective analysis for Research and Applications (MERRA-2) reanalysis (Gelaro et al., 2017) to model with RF at the grid level. We then added remote-sensing-based products from MODIS NBAR and SMAP soil wetness, as well as topographic data, to strengthen the model and provide a finer delineation of environmental gradients based on the literature and expert knowledge (Poulter et al., 2023; Sturtevant et al., 2012). The predictive performance of grid-level models with input variables at their native spatial resolution (except for MERRA-2 variables, which were interpolated to 10 km resolution) was then evaluated. We also compared model performance with those from two additional ML algorithms: support vector machine (SVM) and an artificial neural network (ANN) (Fig. 1, pink boxes). The RF algorithm modeled on all gridded input variables gained the highest mean R² and lowest daily median errors with respect to model predictive performance and was selected for bootstrap modeling and upscaling the 0.098° (∼10 km along the latitudinal length) gridded time series of daily CH₄ fluxes and ensemble uncertainty estimation (Fig. 1, gray boxes).

2.2 Data

2.2.1 Eddy-covariance measurements

The base of our EC data collection stems from a publicly available global synthesis coordination of FLUXNET-CH₄ (Delwiche et al., 2021; Knox et al., 2019), which includes 79 EC tower sites (42 of which are freshwater wetland sites) and 293 site years of data. FLUXNET-CH₄ represents a first compilation of global CH₄ fluxes measured by EC towers (Delwiche et al., 2021; Knox et al., 2019), although more EC data exist outside of the network. We collected both daily and half-hourly data from 44 sites in the northern high latitudes (>45° N), accounting for 167 site years as our base dataset, to which we added data from 6 new sites (31 site years) and additional data to 9 existing sites (21 site years) contributed by principal investigators (PIs; Table S2). In total, we assembled data from 50 EC tower sites in northern latitudes (219 site years), of which 33 are from wetlands (155 site years), with 13 wet tundra sites, 11 fens, and 9 bogs. Data entries with missing data in gridded predictors were excluded, including five wetland sites (FI-LOM, DE-SFN, RU-SAM, RU-VRK, and SE-ST1) where data were collected before SMAP data were available. Another two sites (CA-BOU and RU-COK) were excluded after quality control revealed an instrument anomaly that affected the measurements. As a result, daily and half-hourly EC data from the 26 wetland sites were compiled for analysis from 22 sites in FLUXNET-CH₄ (of which 8 were sites with updated data to recent years, including US-ATQ, US-BEO, US-BES, US-BRW, US-IVO, US-NGB, US-NGC, and US-UAF) and 4 additional sites using information provided directly by PIs (including CA-ARB, CA-ARF, CA-PB1, and CA-PB2), consisting of 138 site years of data in total and representing the largest high-latitude EC data compilation for CH₄ to date (Table S2; see Sect. S2). The sites were distributed among wetland types, including 9 fens, 7 bogs, and 10 wet tundra sites (Fig. 2). RU-CHE and RU-CH2 were two Chersky sites in Eastern Siberian Russia that were located about 600 m from each other to form a paired disturbance experiment. RU-CH2 was a control tower over an undisturbed wetland, whereas RU-CHE was a tower affected by artificial drainage. The aboveground conditions at the two sites were virtually identical, but soil temperature and moisture were different. Drainage caused lower CH₄ fluxes at RU-CHE compared to those at RU-CH2. However, the SMAP data could not discern the drainage impact on soil moisture at the RU-CHE site due to a coarser spatial resolution; thus, the RU-CHE site was excluded from grid-level modeling.

Figure 2 Eddy-covariance tower sites: distribution (>45° N), class, and data size (site years) used in WetCH₄. Colored circles represent EC tower locations and land cover classes, with wetland sites in cyan (wet tundra), yellow (bog), and orange (fen). The circle sizes represent observation years (n) of available CH₄ fluxes at the site (e.g., 1–3 stands for $\mathrm{1}\le n<\mathrm{3}$). The background image shows the estimated maximum annual fractions of wetland cover in 2011–2020 from WAD2M (Zhang et al., 2021).

Half-hourly fluxes acquired from FLUXNET-CH₄ were already gap-filled (see Sect. S2; Irvin et al., 2021). Additional half-hourly fluxes acquired from site PIs were not gap-filled; therefore, we performed per-site gap filling following the FLUXNET-CH₄ approach (Irvin et al., 2021; Knox et al., 2019). Gap-filled fluxes were temporally consistent and agreed with validation data (mean R²=0.68 and mean RMSE = 6 nmol m⁻² s⁻¹; see Sect. S2).

The mean difference in daily mean fluxes between the gap-filled data and the original data converged to −0.2 nmol m⁻² s⁻¹ when there were more than 11 half-hourly EC tower observations in a day but showed substantial bias and larger differences when including days with less than 11 half-hourly observations (Fig. S1 in the Supplement). Therefore, daily data entries were only kept when the number of half-hourly EC tower observations per day was greater than 11. All data were retained for four sites where only daily, quality-filtered data were provided by site PIs (Table S2). As a result, we identified 12 784 daily data entries for upscaling models (Table S2), spanning from 2015 to 2021 with seasonal observation distributions of 44.0 % in June–July–August (JJA), 29.0 % in March–April–May (MAM), 24.5 % in September–October–November (SON), and 2.5 % in December–January–February (DJF) (Fig. S2).

Site-level candidate predictors were identified by their known influences on CH₄ fluxes at multiday to seasonal scales from field control experiments, in situ flux synthesis, and process-based modeling (Bloom et al., 2010, 2017; Knox et al., 2021; Olefeldt et al., 2013, 2017). Only physical variables measured in situ were considered candidate predictors for site-level modeling. In situ candidate predictors that were gap-filled and available in FLUXNET-CH₄ included daily averages of air temperature, soil temperature, air pressure, vapor pressure deficit, relative humidity, latent heat flux, sensible heat flux, longwave incoming radiation, shortwave incoming radiation, net radiation, wind speed, and daily total precipitation (Fig. 1, site-level model solid blue box). We were unable to include water table depth (WTD) or soil water content (SWC) in our site-level model, as they were not available at many sites. However, we explored ML results that included WTD or SWC for a subset of individual sites (36 % of total) where these variables were available (see Sect. S2 for more details).

2.2.2 Reanalysis data and satellite data products

Reanalysis data were used as the gridded input to replace selected predictors at the site level for training the grid-level models and upscaling. These data provided long-term, continuous estimates of nearly all of the candidate predictors of the variables measured in situ (Fig. 1). MERRA-2 is an atmospheric reanalysis of the modern satellite era produced by NASA's Global Modeling and Assimilation Office (Gelaro et al., 2017). We calculated daily means for air pressure, surface air temperature, latent heat flux, sensible heat flux, downward incoming shortwave radiation, downward incoming longwave radiation, and soil temperature at three depths (9.88, 19.52, and 38.59 cm) (Jiao et al., 2023) as well as the relative humidity using the hourly average of surface flux diagnostics, land surface diagnostics, and land surface forcings. The original 0.5°×0.625° resolution data were resampled to 0.5° grids using a bilinear interpolation method in the NASA MERRA-2 web service tool available from the Goddard Earth Sciences Data and Information Services Center (GES DISC). The MERRA-2 data were further bilinearly interpolated from 0.5 to 0.098° grids weighted by the multiple-error-removed improved-terrain digital elevation model (MERIT-DEM) at 90 m resolution that significantly improves elevation estimates in flat terrain over previous DEM products (Yamazaki et al., 2017). Daily time series of the nearest 0.098° grid to each EC location were extracted for grid-level modeling, whereas daily grids were input for the 10 km upscaling products.

To improve the predictive performance of grid-level models, we added remotely sensed biophysical variables, including SMAP soil wetness, MODIS NBAR bands, and topographic data (Fig. 1, Table 1). All remote-sensing products were extracted at daily time steps and their native spatial resolutions at EC tower sites for modeling and were aggregated to 0.098° grids over the study domain for upscaling using Google Earth Engine. We filtered out data gaps in SMAP and MODIS NBAR time series extracted at the native spatial resolution during model training and validation. Gaps in MODIS NBAR were negligible when aggregated from 500 m to 0.098° grids. Gaps in winter SMAP data were filled with zero values to represent frozen soils for upscaling.

Table 1 Description of input variables for grid-level upscaling model.

The SMAP soil moisture product is generated using passive-microwave-radiometer-measured brightness temperature merged with estimates from the GEOS Catchment Land Surface and Microwave Radiative Transfer Model in a soil moisture data assimilation system, providing global products of surface and root zone soil moisture (Reichle et al., 2017). For soil moisture, we employed Level-4 daily soil wetness products (SPL4SMGP.007) from the SMAP mission as proxies for water table depth in the grid-level model (Reichle et al., 2017). Surface, root zone, and soil profile wetness are dimensionless variables that indicate relative saturation for top-layer soil (0–5 cm), root zone soil (0–100 cm), and total profile soil (0 cm to model bedrock depth), respectively. These three variables are originally 3-hourly data at 9 km resolution and were converted to daily means.

Static topographic variables were added as additional attributes in the grid-level model. We used the topographic slope and indices that represent the water flow from the MERIT-DEM based on Geomorpho90m (Amatulli et al., 2020). Two topographic indices were applied: the compound topographic index (cti), which is considered a proxy for the long-term soil moisture availability, and the stream power index (spi, https://gee-community-catalog.org/projects/geomorpho90/, last access: 20 March 2023), which reflects the erosive power of the flow and the tendency of gravitational forces to move water downstream. We tested the impact of elevation on model performance in explaining the inter-site variability in CH₄ upon the current locations of wetland EC sites (see Sect. S6). However, elevation was not considered an ecologically controlling factor for wetland CH₄ fluxes; hence, it was excluded from the input variable importance analysis that ranked the importance of predictors to the prediction accuracy in RF models.

We included MODIS NBAR (MCD43A4v061) products as predictor variables to represent the vegetation productivity in the grid-level model in order to enhance our model predictive performance in vegetated wetlands. The seven NBAR bands (including red, green, and blue bands; two near-infrared bands; and two shortwave-infrared bands) are developed daily at 500 m spatial resolution, using 16 d of Terra and Aqua data to remove view angle effects, and are temporally weighted to the ninth day as the best local solar-noon reflectance (Schaaf et al., 2002; Wang et al., 2018). We did not explicitly include a vegetation productivity variable, as such information is retained in MODIS NBAR, which is used to produce vegetation indices (e.g., EVI) and gross primary production (GPP). Emergent aerenchymous vegetation is another important component in the plant-mediated pathway of CH₄ transport, although it is less represented in existing upscaling models (Table S1).

2.3 Machine learning model

2.3.1 General model design

We used an RF regression algorithm to train site-level and grid-level ML models (Kim et al., 2020). RF regression builds an assembly of independent trees, each of which is trained from a random subset of input data and tested against the rest of the data (Breiman, 2001). A tree grows two leaves when a random selection of subset features reduces the mean-squared error (MSE) of predictions after splitting at a leaf node. Each tree is trained on a bootstrap sample of input data. Trees constructed in this way are less correlated in the ensemble. The generalization error converges as the forest grows to a limit to avoid overfitting. Compared to other ML algorithms, RF has been shown to have better accuracy and lower uncertainty (Irvin et al., 2021; Kim et al., 2020). This approach has been previously applied to upscaling CH₄ fluxes in wetlands and rice paddies across multiple ecosystems (Davidson et al., 2017; Feron et al., 2024; McNicol et al., 2023; Ouyang et al., 2023; Peltola et al., 2019).

A grid search hyperparameter tuning for daily models was performed before predictor selection. We carried out analyses in Python (version 3.6) with the ensemble RF regressor in the “scikit-learn” package (Pedregosa et al., 2011). With all of the predictors and data, hyperparameters were set after tuning for optimized model performance, including the number of trees (n_estimators = 100), number of variables to consider when looking for the best split (max_features = “sqrt”, meaning the square root of the total number of feature variables), the maximum depth of the tree (max_depth = 10), the minimum number of samples required to split a node (min_sample_split = 10), and the minimum number of samples at a leaf node (min_samples_leaf = 4).

For predictor selection and comparisons between the site-level model using in situ variables and the grid-level model using gridded versions of in situ variables, we built the model across all sites and adopted 5-fold cross-validation and “out-of-bag” scores from ensemble trees to evaluate model performance, because, at this stage, we aimed to find physically reasonable variables from in situ measurements and to compare how the differences in scales and measuring methods between in situ predictors and gridded proxies affect the model-learned temporal variability in CH₄ fluxes. A subset of the data were bagged to train each tree in the RF model, with the remaining out-of-bag data used as independent validation data to evaluate the prediction accuracy of each tree, resulting in the average out-of-bag scores of all of the trees in the model. Cross-validation was applied to daily predictions to select variables that can best predict the daily variability in CH₄ fluxes within sites. The changes in model performance after predictor selection and after switching from site-level variables (in situ measurements) to grid-level proxies (reanalysis data) were assessed, which helped quantify differences in model performance when modeling on in situ predictor variables versus modeling on substitute variables at the grid level. Because the data sources of model input from in situ versus from gridded variables were different, we separated site-level and grid-level modeling to ensure that the importance of input features was comparable within a model. The feature importance reflects the relative importance of each input variable in an RF model. It also pertains to the input data distribution and model structure. Therefore, the feature importance by site models can help us identify controlling physical variables but would not necessarily translate to the same rank in the feature importance of grid models, especially when additional gridded variables from remote-sensing products were added to complement the missing controllers from site models.

A summary of input variables for grid-level modeling is provided in Table 1. Although RF can enhance model robustness when collinearity presents in input variables, the collinearity could affect the interpretation of feature importance measured by the impurity decrease in RF models. Therefore, we first built a baseline grid-level model with independent variables after a pairwise Pearson correlation test (Table S3) to exclude covariates. We grouped significantly correlated variables (p<0.001, r>0.8; white grids except for those on the diagonal line in Table S3), thereby forming three groups: SMAP soil moisture variables in group 1 (we also included surface soil moisture that was significantly correlated with the other two soil moisture variables and r>0.7); air temperature (tas), downward longwave radiation (rsdl), spfh, and soil temperatures (ts1, ts2, and ts3) in group 2; downward shortwave radiation (rsds) and latent heat (le) in group 3. We then selected the one most important variable in each group for the baseline models according to the feature importance of modeling on all predictor variables (Fig. S14). The rest of the variables out of the groups were included in the baseline features. The resulting baseline features included air pressure (pa), latent heat flux (le), sensible heat flux (h), soil temperature (ts2), root zone soil wetness (sm_r_wetness), slope, spi, and cti. Then, we designed four additional different model settings by changing predictor variables, including (1) baseline variables plus covariates, (2) only variables from MODIS NBAR, (3) baseline variables plus NBAR bands, and (4) all predictor variables. In this forward feature selection process, we evaluated the impacts of adding constraint variables from remote-sensing products on model performance.

Model predictive performance evaluates the accuracy of a model to predict at a new site without any prior knowledge. For the spatial predictive performance evaluation of grid-level ML models, we used a nested leave-one-site-out cross-validation scheme (hereafter referred to as LOOCV). Such a scheme selects one site to use as independent validation data to evaluate models trained and tested with data from the remaining sites, repeating the process for all sites. Without any prior knowledge of the validation site added to a model, the LOOCV scheme can (1) assess the predictive ability of the model in a new place and (2) evaluate the uniqueness of a site in the dataset. Similar forms of spatial LOOCV have been used to evaluate upscaling models for global or regional CO₂ and CH₄ (McNicol et al., 2023; Peltola et al., 2019; Virkkala et al., 2021). The validation of the upscaling model was not only performed with respect to daily predictions but also on monthly means. The predictive performance of the upscaling model on the monthly variability in CH₄ fluxes and spatial variability across sites is important for studies that vary with respect to their temporal and spatial scales.

Model predictive performance was assessed using three evaluation metrics: mean absolute error (MAE), root-mean-square error (RMSE), and R² score. Daily modeled CH₄ fluxes were compared to EC observations at each validation site. The evaluation metrics were calculated at daily and monthly scales for each site separately to examine the model performance by general wetland types and for all sites pooled together to evaluate the overall performance and undertake a comparison with existing studies. Squared error metrics are more sensitive to outliers and highly skewed data, which is often the case with CH₄ fluxes. Therefore, we selected both the MAE and RMSE to quantify the errors. The mean error (ME) between model predictions and validation data was calculated, representing the systematic bias in predicted fluxes. The standard deviation of model residuals was also included to measure the spread of the residuals. This matches the RMSE when the ME equals zero.

Two additional ML algorithms were compared with RF: SVM and ANN. SVM is efficient with sparse data where the dimension of the input space is greater than the number of training samples (Kuter, 2021). While the training process of ANN is expensive and time-consuming, it can develop deep networks with growing training data which may increase predictive performance (Saikia et al., 2020). We used support vector regression to model CH₄ fluxes with the same predictor variables and dataset as used in the ensemble RF regressions. Multilayer perceptron regressor is an implementation of an ANN model that adjusts the weights of neurons using backpropagation to improve prediction accuracy. It uses the square error as the loss function and a stochastic gradient-based optimizer “Adam” for weight optimization. We used two hidden layers in the ANN model, each with 50 neurons. Data from all variables were normalized to achieve the best model performance for SVM and ANN.

2.3.2 CH₄ flux upscaling

We trained 500 ensemble RF models with all gridded predictors of the grid-level models from the general model design and with data from all sites for upscaling daily CH₄ fluxes. Each RF model was trained with the same optimized hyperparameters and different bootstrap samples. Ensemble models were then applied to 0.098° gridded predictors to produce the upscaling CH₄ flux intensities from the means of the 500 predictions and the prediction uncertainty from the standard deviations. Given that the CH₄ fluxes were modeled with data from the wetland EC sites, a wetland extent map was also needed to constrain the areas when scaling grid emissions (see Sect. 2.4). Final CH₄ emission and uncertainty maps associated with wetland extents were the result of multiplying the predicted means and standard deviations of flux intensities by wetland areas. All wetland maps were resampled to 0.098°×0.098° resolution with a conservative remapping method for producing the emission products.

2.4 Wetland extent maps and benchmark estimates of wetland CH₄ emissions

Wetland extent maps were applied to scale the modeled CH₄ flux intensities to the region. The Wetland Area and Dynamics for CH₄ Modeling (WAD2Mv2) product, representing spatiotemporal patterns of inundated vegetated wetlands at 0.25° resolution, was selected as the reference for dynamic wetland areas in this study (Zhang et al., 2021). Active- and passive-microwave-detected inundation data combined with static wetlands were used to delineate the monthly dynamics of wetland inundation between 2000 and 2020. Open waterbodies such as lakes, rivers, reservoirs, coastal wetlands, and rice paddies were excluded. We used monthly mean WAD2M fractions between 2010 and 2020 to represent seasonal wetland dynamics. Emission estimations are subject to differences in the wetland extent between maps (Saunois et al., 2020). We used monthly means of the Global Inundation Extent from Multi-Satellites (GIEMS2) product (Prigent et al., 2020) to represent temporal patterns in the restricted wetland extents at 0.25° resolution. The coarse-resolution maps were resampled to 0.098°×0.098° grids using the nearest-neighbor method. The static Global Lakes and Wetlands Database version 1 (GLWDv1) Level-3 1 km resolution map excluding the lakes, rivers, and reservoirs classes (Lehner and Döll, 2004) was included to quantify the upper limit of wetland cover. For all explicit GLWDv1 wetland classes, we assumed a 100 % wetland coverage in the original pixels, except for “intermittent wetland/lake” for which we assumed a 50 % coverage; for GLWDv1 classes represented as extent ranges, we used the average value of the range (i.e., 75 % for 50 %–100 % wetland, 37 % for 25 %–50 % wetland, and 12 % for 0 %–25 % wetland). To support domain emission comparisons, wetland cover was also extracted from the updated GLWD version 2 dataset (GLWD v2; Lehner et al., 2025) which provides the spatial extent of 33 waterbody and wetland classes at 500 m spatial resolution. All freshwater wetland classes that occur in our study area (classes 8–25) from GLWD v2 were included (i.e., excluding rivers, lakes, reservoirs, and other permanent open-water bodies as well as coastal saline/brackish wetlands). The original wetland areas per GLWD v2 pixel were summed across all included classes to derive a total wetland area per pixel. Furthermore, a regional freshwater wetland distribution dataset was calculated from a permafrost-region-specific land cover map (CALU – circum-Arctic land cover units) which classified 23 land covers including 3 wetland classes and 10 moist to wet tundra classes at 10 m resolution and aggregated to 1 km with the majority class (Bartsch et al., 2024). This regional wetland map was applied for CH₄ emission estimation in the North Slope region of Alaska to assess the impacts of different wetland maps on emission estimates in this area when compared against airborne measurements. Wetland areas from the finer-resolution maps were aggregated to 0.098°×0.098° grids for emission calculations.

We compared WetCH₄ emissions with benchmark domain or regional estimates from bottom-up process models, top-down atmospheric-observation-based inversions, and existing upscaling studies. We acquired data for the study domain from the ensemble mean of bottom-up process-based models from the Global Carbon Project (GCP) (Zhang et al., 2025) and the extended ensemble of wetland CH₄ estimates that were priors for the top-down GEOS-Chem atmospheric chemical and transport model (WetCHARTs) (Bloom et al., 2017; Friedlingstein et al., 2022). We also included the atmospheric inversions of northern high latitudes from an assimilation CarbonTracker-CH₄ system (Bruhwiler et al., 2014; update at https://gml.noaa.gov/ccgg/carbontracker-ch4/carbontracker-ch4-2023/, last access: 16 May 2023). We compared WetCH₄ with existing upscaled products of monthly CH₄ wetland fluxes based on Peltola et al. (2019) for the study domain. For regional wetland hotspots, CH₄ flux estimates were obtained from Carbon in Arctic Reservoirs Vulnerability Experiment (CARVE), which measured total atmospheric columns of CO₂, CH₄, and carbon monoxide over northern Alaska in spring, summer, and early fall between 2012 and 2014 (Chang et al., 2014; Miller et al., 2016). These were used to verify our seasonal emission estimates over the North Slope region (Zona et al., 2016).

3 Results

3.1 Model validation

3.1.1 Site-level modeling

Site-level modeling used all wetland sites to build an RF model and identified the 10 most important variables measured in situ that, if left out, decreased the valuation score of the model by more than 90 % based on the mean decrease in impurity (Fig. S3). With bootstrap sampling and using all candidate predictors (Fig. 1) in the model, the out-of-bag RMSE of the site-level model was 30.22 nmol m⁻² s⁻¹, and the out-of-bag R² between observed daily means of CH₄ fluxes and prediction was 0.73. Modeling with the 10 most important variables at the site level resulted in similar model performance, with an out-of-bag RMSE of 30.43 nmol m⁻² s⁻¹ and an out-of-bag R² of 0.73. Site-level model performance converged as the increment of predictor variables ordered by the importance rank (Fig. S4). We then tested building separate models according to wetland types because distinct CH₄ fluxes have been observed from wet tundra (Fig. S5; mean ± standard deviation: 13±14 nmol m⁻² s⁻¹), bogs (22±26 nmol m⁻² s⁻¹), and fens (56±88 nmol m⁻² s⁻¹). The out-of-bag R² (RMSE) was 0.85 (7.2 nmol m⁻² s⁻¹) for bogs, 0.84 (27.7 nmol m⁻² s⁻¹) for fens, and 0.57 (34.3 nmol m⁻² s⁻¹) for wet tundra. Modeling with the selected 10 predictors resulted in an out-of-bag R² (RMSE) of 0.84 (7.6 nmol m⁻² s⁻¹) for bogs, 0.84 (27.9 nmol m⁻² s⁻¹) for fens, and 0.53 (36.3 nmol m⁻² s⁻¹) for wet tundra. Next, we tested whether the inclusion of non-wetland sites (upland and rice sites) would affect model performance. This resulted in an out-of-bag R² decrease to 0.56 and an RMSE increase to 38.86 nmol m⁻² s⁻¹, suggesting that a generalized ML model over all land cover classes is not practical to reliably predict CH₄ fluxes with the current set of predictors and available data. This is most likely due to the distinctive features of CH₄ emissions between wetland and non-wetland classes (Fig. S5).

3.1.2 Grid-level modeling and remote-sensing constraints

Substituting in situ measurements of selected predictor variables with gridded MERRA-2 variables slightly reduced model accuracy. The out-of-bag R² decreased by 9.6 % to 0.65 and the RMSE increased by 15 % to 34.9 nmol m⁻² s⁻¹ compared to the site-level model. The coarse-resolution MERRA-2 reanalysis data capture less spatial variability in the selected physical variables and are less accurate at the grid level compared to in situ EC measurements.

Figure 3 Grid-level modeling. Panel (a) presents the distribution of the R², MAE, RMSE, and ME for all sites (size = 25) in a LOOCV scheme based on gridded data using five model settings: RF modeled using only MODIS NBAR bands, baseline features (MERRA-2 air pressure, latent heat flux, sensible heat flux, soil temperature, SMAP root zone soil wetness, topographic slope, spi, and cit), baseline features plus MODIS NBAR bands, baseline features plus correlated variables within the MERRA-2 and SMAP dataset, and all gridded input variables together. The model settings are ranked by the mean R², from lowest (left) to highest (right). Panel (b) shows the mean variable importance of baseline models (last column) in the LOOCV scheme and at each site (columns labeled with validation site ID). The values in each column are the means of the accumulation of the impurity decrease when a variable was taken out in the trees of an RF model, representing the importance of such variable to the model. The variable names and descriptions refer to Table 1.

Adding remote-sensing constraints to the gridded variables can improve the model predictive performance and reduce errors. Modeling on baseline features explained, on average, 46 % of daily CH₄ fluxes' variability at validation sites with the largest range of errors (Fig. 3a). The medians in the baseline model of the R², MAE, RMSE, and ME under the LOOCV scheme were 0.5, 16.4, 21.0, and 6.4 nmol m⁻² s⁻¹, respectively. Adding NBAR or covariates from MERRA-2 and SMAP input variables returned a higher mean R² or slightly lower mean errors than the baseline model, whereas modeling with all gridded input variables (the “all” model setting) achieved the highest mean R² of 0.51 with a comparable mean MAE (23.6 nmol m⁻² s⁻¹), RMSE (32.1 nmol m⁻² s⁻¹), and ME (0.9 nmol m⁻² s⁻¹) (Table S4). Although modeling with baseline features and covariates (the “base + CoVar” setting) received a comparable mean R² to modeling all variables, the latter had a higher median R² (0.53) and lower median errors (MA = 14.1 nmol m⁻² s⁻¹, RMSE = 19.8 nmol m⁻² s⁻¹, and ME = 4.0 nmol m⁻² s⁻¹). Our results suggest that including remote-sensing constraints or covariates improved the models' ability to predict spatial variability in wetland CH₄ fluxes and reduced prediction errors. These results confirm our selection of predictor variables for the upscaling model (Table 1).

Figure 4 Model predictive performance evaluation on RF-modeled CH₄ fluxes at the grid level under the LOOCV scheme: box plots of R², MAE, RMSE, and ME across validation sites by wetland types with mean values denoted using black squares at daily, weekly, and monthly (top, middle, and bottom panel, respectively) time steps.

The average importance of the baseline features shows their influence on the grid-level model predictive performance (Fig. 3b). The importance of independent predictors under the LOOCV scheme, although slightly varied between models, agreed with respect to selecting the MERRA-2 soil temperature (ts2) as the primary driver in predicting daily CH₄ fluxes in northern wetlands, followed by the SMAP root zone wetness (sm_r_wetness). The eight baseline features accounted for a 99 % reduction in the mean validation score of the baseline models. The average importance of all gridded variables used for upscaling (Fig. S14) was consistent with baseline models, emphasizing the importance of soil temperatures and root zone wetness. Additionally, air pressure and topography also contributed to explaining the daily variability in CH₄ fluxes. Nevertheless, all variables contributed to predicting variability in CH₄ fluxes, showing the complexity of environmental factors that would affect the rates of CH₄ production and the process of gas exchange.

Daily mean CH₄ fluxes exhibited great variability in wetlands across space and time (mean = 35 nmol m⁻² s⁻¹, σ=65 nmol m⁻² s⁻¹; Fig. S3). The model predictive performance (Fig. 4) was calculated for each site, and the average performance with respect to the daily variability in CH₄ fluxes was best at wet tundra sites with a mean R² of 0.56, followed by bog sites (0.51) and fen sites (0.45). Due to the large variability in fen daily fluxes, errors in daily predictions were highest for fen sites (mean RMSE = 54.2 nmol m⁻² s⁻¹ and mean MAE = 37.8 nmol m⁻² s⁻¹), followed by bog sites (mean RMSE = 27.6 nmol m⁻² s⁻¹ and mean MAE = 22.5 nmol m⁻² s⁻¹), whereas they were lowest for wet tundra sites (mean RMSE = 13.5 nmol m⁻² s⁻¹ and mean MAE = 10.3 nmol m⁻² s⁻¹). Our slight model overestimation of daily fluxes (mean ME = 0.9 nmol m⁻² s⁻¹) was driven by an underestimation of fen sites (mean ME = −12 nmol m⁻² s⁻¹) versus an overestimation of bog (mean ME = 14 nmol m⁻² s⁻¹) and wet tundra (mean ME = 3 nmol m⁻² s⁻¹) sites.

Figure 5 Example model predictive performance for seasonal cycles of daily FCH₄ at the validation sites of CA-SCB, CA-ARF, and US-NGB, representing bog, fen, and wet tundra, respectively.

Model predictive performance on aggregated monthly means of CH₄ fluxes increased by 37 % compared with daily means (mean R²=0.70; Fig. 4, Table S4). This improvement may be attributed to a better representation of the environmental conditions' average state over a month by the input variables compared to the daily variability. Performance was higher for wet tundra (mean R²=0.73) and bogs (mean R²=0.73), whereas it was lower for fen sites (mean R²=0.64, Fig. 4). Mean errors in monthly mean predictions were as follows: RMSE = 28.1 nmol m⁻² s⁻¹, MAE = 21.4 nmol m⁻² s⁻¹, and ME = 0.37 nmol m⁻² s⁻¹ (Table S4). Prediction residuals of daily and monthly CH₄ fluxes (Fig. S6) showed normal distributions for wet tundra sites, indicating that the spread of residuals resulted from random errors that increased with the flux magnitude. The residuals had a skewed normal distribution for bog sites, indicating likely overestimation. The long left-hand tails in the prediction residuals indicated that the intense emission fluxes from fens during summer peaks were underestimated (Fig. S6).

Site-by-site validation of daily flux predictions varied greatly between individual sites (Figs. 5 and S7). For example, US-UAF, an EC site in Interior Alaska with mature black spruce cover and full understory vegetation and mosses over permafrost (Ueyama et al., 2023a), which is the only one of the five forest bog sites in our dataset that had low CH₄ fluxes and weak seasonal cycles (less than 10 nmol m⁻² s⁻¹), was significantly overestimated by our model (RMSE = 58 nmol m⁻² s⁻¹ and MAE = 53 nmol m⁻² s⁻¹). Permafrost presence and groundwater below the soil surface may explain the low fluxes at this site (Iwata et al., 2015; Ueyama et al., 2023b).

Figure 6 Mean annual wetland CH₄ fluxes: panels (a)–(c) contain WetCH₄ upscaled fluxes between 2016 and 2022 weighted by wetland fractions for three wetland maps – WAD2Mv2, GIEMS2, and GLWDv1; panels (d)–(f) contain the bottom-up GCP ensemble mean, WetCHARTs, and top-down estimates of CarbonTracker-CH₄ natural microbial emissions.

3.2 Upscaled wetland CH₄ emissions

3.2.1 Wetland-area-weighted CH₄ emissions

Upscaled daily CH₄ fluxes were weighted by wetland fraction to estimate gridded daily CH₄ fluxes from northern wetlands based on WAD2Mv2, GIEMS2, and GLWDv1 between 2016 and 2022 (Fig. 6) and on GLWD v2 for comparison. The mean annual emissions and the RF-model-associated uncertainties are summarized with different wetland maps in Table S5. The estimate from WetCH₄ with WAD2Mv2 was 22.8±2.4 Tg CH₄ yr⁻¹, comparable to UpCH₄ (23.5±5.8 Tg CH₄ yr⁻¹). With GIEMS2, WetCH₄ estimated a minimum annual emission of 15.7±1.8 Tg CH₄ yr⁻¹. With GLWDv1 and GLWD v2, WetCH₄ estimated potential annual emissions of 46.0±5.1 Tg CH₄ yr⁻¹ and 51.6±2.2 Tg CH₄ yr⁻¹ for 2016–2022, respectively. The spatial patterns were similar to the post-2016 mean annual fluxes from the GCP process-based model ensemble means (28.6±21.6 Tg CH₄ yr⁻¹ for 2016–2020), WetCHARTs (29.5±30.0 Tg CH₄ yr⁻¹ for 2016–2019), and atmospheric inversions of CarbonTracker-CH₄ (40.9 Tg CH₄ yr⁻¹ for 2016–2022), highlighting the high-emission areas in the Hudson Bay Lowlands and Western Siberian Lowlands. The emissions from WetCH₄-GIEMS2 were lower in these two hotspots than other estimates. Differences in the distribution of CH₄ emissions between wetland products reflect the influence of wetland dynamics. Mean monthly wetland inundations are provided by WAD2Mv2 and GIEMS2, which set the dynamic limits for the wetland boundaries of the CH₄-emitting surface. While emissions resulting from inundation were captured, it appeared that saturated or wet subsoil conditions were not well represented by WAD2M and GIEMS2, resulting in low emissions in wet but non-inundated tundra (i.e., the Alaskan North Slope). To address this, we incorporated wetland fractions from the CALU high-resolution wetland map (Bartsch et al., 2024) that was specifically produced for the permafrost region in order to estimate Alaskan North Slope emissions. Wetland fractions from GLWD (both v1 and v2) represent a static maximum wetland distribution throughout time. Thus, estimates from GLWD may represent the upper bounds for all northern wetlands under contemporary conditions.

We compared spatial distributions of our upscaled fluxes (WetCH₄) with two alternative upscaled datasets. Using the same wetland weights, our product showed similar spatial patterns to UpCH₄ (McNicol et al., 2023) and the upscaled fluxes from Peltola et al. (2019) (Fig. S9). Spatially, the maximum mean flux for WetCH₄ with WAD2Mv2 for 2016–2022 was 69 mg CH₄ m⁻² d⁻¹, whereas UpCH₄ produced a maximum mean flux of 88 mg CH₄ m⁻² d⁻¹ between 2016 and 2018. While all three products predicted concentrated CH₄ exchange in the Hudson Bay Lowlands and Western Siberian Lowlands and low fluxes in the Western Canadian Arctic tundra, WetCH₄ predicted lower fluxes in the forested wetlands of Western Canada compared with UpCH₄ (Fig. S9a and b). With GLWDv1, WetCH₄ predicted similar fluxes to those of Peltola et al. (2019), with the exception of a number of potent emitting grids in the Western Siberian Lowlands (Fig. S9c and d) and a maximum mean flux of 132 mg CH₄ m⁻² d⁻¹ from WetCH₄.

Figure 7 Multiyear average seasonal cycles of wetland CH₄ emissions. Panel (a) presents a comparison of ML upscaled mean seasonal cycles in reference wetland areas (WAD2Mv2) with the cycles from process-based models in the northern mid-high latitudes (45–90° N). Panel (b) shows the same comparison for northern high latitudes (60–90° N) and the addition of atmospheric CarbonTracker-CH₄-attributed microbial emissions (2016–2022). Panel (c) presents a comparison of three ML upscaled mean seasonal cycles of CH₄ emissions with different wetland area maps (WAD2Mv2, WAD2Mv2 maximum extent, GIEMS2, and GLWDv1). Panel (d) shows a comparison of WetCH₄ mean seasonal cycles over land (black line), WetCH₄ weighted by wetland from the CALU map (light-green line), or WetCH₄ weighted by fractions from WAD2Mv2 (dark-green line), with estimates of CH₄ fluxes in growing seasons from CARVE retrievals in the North Slope area of Alaska (Zona et al., 2016).

3.2.2 Seasonal cycles of wetland CH₄ emissions

Mean seasonal cycles of wetland CH₄ emissions were consistent with bottom-up estimates in the domain and top-down inversions in high latitudes (Fig. 7). The amplitudes of two ML-based estimates agreed in the domain (WetCH₄ and UpCH₄ both within WAD2Mv2 wetland areas) and were lower than the ensemble means of GCP or WetCHARTs estimates during the growing season (Fig. 7a). In the northern high latitudes (60–90° N), the amplitudes of this study closely agree with WetCHARTs, and both amplitudes were lower than the ensemble means of GCP in the growing season (Fig. 7b). Our emissions in JJA were lower than the emissions attributed by the atmospheric inversion of CarbonTracker-CH₄, which does not discriminate between wetland and open-water sources. We did not use comparisons with CarbonTracker-CH₄ for 45–90° due to likely considerable contributions from aquatic systems and other non-wetland factors in the inversion estimates. Notably, uncertainties between ML-based approaches with the same wetland extents showed less variation than those between process-based models, especially during the growing season. The phase of our estimates (WetCH₄) agreed with bottom-up and top-down models, peaking in July followed by August (Fig. 7a and b), whereas UpCH₄ showed a month lag, probably due to the 2- or 3-week lag of the predictor variables selected in UpCH₄ (McNicol et al., 2023). Peak fluxes in July and August were commonly seen in tower measurements.

The seasonality in upscaled wetland CH₄ emissions corresponded to the intensities of fluxes and dynamics of wetland areas. We compared mean seasonal cycles of upscaled products with different dynamic or static wetland maps to constrain the impacts of wetland areas (Fig. 7c). As observed in spatial distributions (Fig. 7a and c), emissions from the potential emitting surface (WetCH₄_GLWDv1) were 95 % higher than those from reference inundated wetlands (WetCH₄_WAD2Mv2) during the growing season, and they doubled in winter. Within the GLWDv1 emitting surface, WetCH₄ predicted higher emissions than Peltola et al. (2019) in July (43 %), August (43 %), December (41 %), and January (61 %), but it predicted 15 % lower emissions in October. We decoupled the mean annual seasonal cycle for WAD2M from the emission seasonality using a fixed maximum WAD2M extent. The addition of maximum annual wetland extent further constrains the limitations of seasonal WAD2M extents in underestimating the CH₄-emitting surface for northern high-latitude wetlands, especially in the cold seasons. The resulting seasonal emissions, primarily driven by soil temperatures and moisture, manifested elevated emissions in all months and an intensified seasonal cycle. Reported emissions (Zona et al., 2016) and large bursts (Mastepanov et al., 2008) from the freezing active layer in permafrost areas in October (zero-curtain period) may not be well captured by our ML model. The differences in wetland areas between the two dynamic products (WAD2Mv2 and GIEMS2) mostly affected emissions in May and June in WetCH₄, but they significantly affected emission magnitudes in UpCH₄. Despite the differences in wetland areas, the phases of emissions cycles of WetCH₄ were consistent with those from Peltola et al. (2019), whereas UpCH₄ again lagged by a month.

We compared upscaled seasonal cycles with CH₄ fluxes estimated from regional airborne measurements taken during CARVE campaigns over the Alaskan North Slope (Fig. 7d). Given that the wetland area in this region is uncertain (Miller et al., 2016), we computed mean seasonal cycles over land assuming that all land in this area is water-saturated in the soil, over freshwater wetlands of CALU, and over WAD2M and HydroLAKES, representing three different scenarios. In the lowland area of the North Slope (74 295 km² spanning between 69.8 and 71.4° N and between 164.4 and 152.7° W), the wetland area was estimated at 10 611 km² from CALU, 4800 km² from GLWD v2, and 4049 km² from the maximum extent month in July of WAD2Mv2, respectively. The range of our upscaled estimates aligned with regional emissions derived from CARVE measurements. Chang et al. (2014) estimated 7±2 mg CH₄ m⁻² d⁻¹ of mean CH₄ fluxes during the growing season in the North Slope region from the column analysis of CARVE data. The mean fluxes (May–September) of WetCH₄ with CALU were estimated at 7.3±0.8 mg CH₄ m⁻² d⁻¹ (5.5±0.6 mgC CH₄ m⁻² d⁻¹), which is within the range of various CARVE estimations (Miller et al., 2016). The landscape is in the biome of the Arctic coastal tundra and is covered by sedges, grasses, mosses, and dwarf shrubs. A large number of lakes and freshwater ponds are scattered across the area. Studies carried out on the western Alaska lowland of the Yukon–Kuskokwim Delta found aquatic fluxes that were about 10 times higher than in wet tundra during September (Ludwig et al., 2023), suggesting that a major source of the airborne fluxes missing in WetCH₄ in the late growing season can be attributed to open-water fluxes. Remarkable increases could occur in summer and winter if we assume wetland over this region, as indicated by the range between the green and black lines in Fig. 8d. However, future emissions due to permafrost thaw still depend on the hydrological changes in the landscape.

Figure 8 Wetland CH₄ (a) annual emissions and associated uncertainties (colored shades) and (b) variations relative to multiyear means in the research domain (45–90° N). Wetland area data applied in WetCH₄ and UpCH₄ were WAD2Mv2. The time periods of the multiyear means are as follows: WetCH₄ (2016–2022), UpCH₄ (2016–2018), GCP bottom-up ensemble mean (2016–2020), and WetCHARTs (2016–2019).

3.2.3 Interannual variations in wetland CH₄ emissions

The mean annual emissions from ML-based estimates with WAD2M were lower than the GCP bottom-up ensemble mean and WetCHARTs over different years from 2016 forward in time (Fig. 8a). All products demonstrated similar emission patterns for the domain in the interannual trends and variations, which were highest in 2016 and lower for the 3 years from 2017 to 2019 (Fig. 8). The interannual variations in WetCH₄ were driven by the interannual variability in the upscaled fluxes, as only multiyear mean seasonal dynamics from WAD2Mv2 were used. All products identified intensified emissions in 2016, as indicated by the variations relative to period means (Fig. 8b). Higher-than-period-average emissions in 2020 were also modeled by WetCH₄ and the ensemble GCP.

Figure 9 Interannual variations and variability in subregions predicted by WetCH₄ with WAD2Mv2, GLWDv1, and GIEMS2, respectively: (a) interannual variations with respect to period means (2016–2022); (b) relative variability as the percentage of its period mean. Delta on the y axis denotes the annual emissions minus mean annual emissions in the 2016–2022 period. The box plots show the first quartile, the median, and the third quartile of the data, with the whiskers denoting 1.5× the interquartile range below or above the first or third quartile, respectively.

Subregional annual emissions and interannual variability (Fig. 9) for WetCH₄ were calculated for eight subregions in the northern high latitudes (Fig. S11): Siberian tundra, Eastern Siberia, Western Siberia, Fennoscandia, Canadian tundra, Eastern Canada, Western Canada, and Alaska. The main differences in WetCH₄-estimated emissions between WAD2Mv2 and GLWDv1 occurred in the Eastern Siberia, Eastern Canada, Western Canada, and Alaska subregions. However, interannual variabilities were similar. Interannual variation from Western Siberia accounted for 51 % of the variation in domain emissions (Fig. 9a). The positive change in Eastern Canada canceled the negative change in Western Siberia in 2021, resulting in low variability in the domain emissions for that year (Fig. 8). The relative interannual variability, which was calculated as the percentage of a subregional variation to its period mean emission, was attributed to variability from Western Siberia, Fennoscandia, Western Canada, and Alaska (Fig. 9b).

4 Discussion

This study provides new estimates of daily-scale 10 km wetland CH₄ fluxes for the northern terrestrial wetland region, upscaled from EC data. The upscaling framework was driven by MERRA-2 meteorological variables and soil temperatures and constrained by satellite products from SMAP soil moisture and MODIS NBAR, resulting in a good prediction accuracy (mean R²=0.70 and mean MAE = 27 nmol m⁻² s⁻¹) for monthly mean fluxes. Model agreement worsened at daily and weekly time steps due to higher variability in CH₄ fluxes at finer temporal resolutions. In our framework, we applied a rigorous criterion on the counts of half-hourly observations to control the selection quality of daily gap-filled data, which may filter out errors introduced by the gap-filling process or lack of observations for calculating daily means. The improvement in model performance can be partly attributed to the inclusion of soil temperature, satellite assimilation of soil moisture, and MODIS vegetation reflectance in the framework that represents controlling factors or proxies for CH₄ fluxes recognized in field experiments and synthesis studies (Fig. 3).

4.1 Important drivers to improve RF model predictive performance

Soil temperature plays an important role in microbial growth and dormancy (Chadburn et al., 2020) and exponentially affects microbial CH₄ emission rates, although the temperature sensitivity varies across space and time (van Hulzen et al., 1999; Knox et al., 2021). In northern wetlands, soil temperature is often more spatially variable relative to air temperature due to snow insulation and the active-layer depth (Smith et al., 2022; Wang et al., 2016; Yuan et al., 2022); thus, it should be considered in upscaling models. Compared to air temperature or land surface temperature, which have been used in previous upscaling studies (McNicol et al., 2023; Peltola et al., 2019), the inclusion of MERRA-2 soil temperatures in WetCH₄ likely contributed to a higher model predictive performance, although the impacts of the scale mismatch between the native MERRA-2 spatial resolution and the local footprints on the upscaled fluxes were not quantified. Independent validation studies found significant correlations in the temporal trend and seasonal cycles between MERRA-2 soil temperatures and in situ observations (Li et al., 2020; Ma et al., 2021) in the US and midlatitude Eurasia. However, lower correlations and overestimated monthly variability were found in the cold season throughout the Arctic (Herrington et al., 2024). This suggests that the impact of the uncertainty in MERRA-2 soil temperatures was concentrated in the cold season, when CH₄ fluxes were low. The agreement between the ensemble means of soil temperatures from eight reanalysis and land data assimilation system products and station measurements improved in the pan-Arctic region (Herrington et al., 2024), suggesting the potential to reduce upscaling uncertainty forced by the ensemble mean of reanalysis datasets.

Emergent vegetation with aerenchyma affects the recent substrate availability and the plant-mediated transport of CH₄ (Kyzivat et al., 2022; Melack and Hess, 2023). We used the full land bands of the MODIS NBAR product, rather than the derived vegetation indices employed in previous upscaling studies, as signals indicating that wetland vegetation functional characteristics may be lost when merging bands to derive simple vegetation indices (Chen et al., 2013). In our study, the near-infrared and shortwave-infrared bands (NBAR bands 2, 5, and 7) presented relatively high importance in the RF model due to their associations with vegetation characteristics and water table dynamics in northern peatlands (Baskaran et al., 2022; Burdun et al., 2023). Satellite inputs provide high-spatial-resolution constraints on the environmental variability and help reduce spatial model predictive errors (Fig. 3), indicating the requirement for high-spatial-resolution driving input to accurately model wetland CH₄ fluxes (Elder et al., 2021).

Surface and root zone soil moisture are important controls on ecosystem anaerobic metabolism. Low soil moisture implies aerobic conditions and allows methanotrophic bacteria to consume CH₄, whereas high soil moisture enables CH₄ production and suppresses consumption (Liebner et al., 2011; Olefeldt et al., 2013; Spahni et al., 2011). Soil wetness estimated in the root zone and the profile from SMAP measurements may be able to capture water table dynamics; hence, it ranked as important with respect to the WetCH₄ model performance. Validation of the SMAP Level-4 soil moisture data assimilation product has shown that it meets the performance requirement of an unbiased root-mean-square error < 0.04 m³ m⁻³ (Colliander et al., 2022). However, the validation sites are mostly located in North American grassland, cropland, and shrubland, requiring more in situ soil moisture observations in high-latitude tundra and peatland. Regional validation studies suggested that uncertainties in satellite-derived soil moisture including SMAP at high latitudes were high (Högström et al., 2018; Wrona et al., 2017) and remained to be addressed.

Underground processes of CH₄ production and oxidation are difficult to model (Ueyama et al., 2023b), especially for seasonal cycles in the northern high latitudes. A hysteresis effect that manifests intra-seasonal variability in the dependence of CH₄ fluxes on temperature has been observed at EC sites (Chang et al., 2021), but it was not reproduced in WetCH₄. Positive hysteresis and the difference in frozen status from topsoil to deep soil during the fall freeze results in zero-curtain periods that have been observed in high-latitude tundra (Bao et al., 2021; Zona et al., 2016), the occurrence of which was subsequently underestimated in our model.

The amount of additional substrate available for methanogenesis due to soil freezing/thawing, missing in our framework, could be a controlling factor on the occurrence of this phenomenon. Higher substrate availability elevates methanogen abundance and activities during the fall freeze (Bao et al., 2021). However, spatially explicit substrate data are not available. Using proxies such as net primary production or the EVI for substrate availability might be oversimplified (Larmola et al., 2010; Li et al., 2016; Peltola et al., 2019). In addition, the uncertainty in the deep-soil temperature of training inputs in late fall may hinder the model's ability to capture patterns of high emissions during zero-curtain periods observed in the Alaskan tundra (Fig. S10). More temporally accurate soil temperature data are needed to delineate the soil freezing progress and properly constrain predictions of CH₄ emission during the cold season (Arndt et al., 2019). The UpCH₄ results (McNicol et al., 2023) also suggest that simply imposing lags on temporal predictors in RF cannot capture complex intra-seasonal variability due to the complicated lag effects interacting with the water table depth (Turner et al., 2021). Without timestamps for predictors, RF treats time series' fluxes independently, which may limit its predictive performance. Deep learning models designed to account for temporal progress in data, such as long short-term memory (LSTM) neural networks, may improve the modeling accuracy of seasonal cycles (Reichstein et al., 2019; Yuan et al., 2022).

4.2 Data limitations in current EC CH₄ observations

Data deficiency in EC CH₄ flux observations in winter and in underrepresented areas limited the RF model's extrapolation ability. Data abundance and representativeness across space, time, and wetland types drives model performance and the ability to extrapolate for the data-driven approach. The 26 wetland EC sites included in this study are largely located in Fennoscandia, Eastern Canada, and Alaska (Fig. 2), leaving some regional emission hotspots underrepresented. For instance, the Western Siberian Lowlands, a large wetland complex and a major contributor of interannual variation in CH₄ in the region, comprise a region with little data. The nearest site (RU-VRK, not included in this study due to the fact that observations occurred before our study period) is situated on the western side of the Ural Mountains, within the Usa River depression. Cold-season emissions could contribute a substantial fraction of the annual Arctic tundra CH₄ budget (Mastepanov et al., 2008; Mavrovic et al., 2024; Zona et al., 2016). However, after filtering, 23 % of the EC data in high latitudes (>60° N) were recorded between November and March, which could be insufficient for accurately modeling and upscaling zero-curtain-period fluxes.

A total of 10 bog and fen sites used for modeling contain all-season daily flux records with more than 11 half-hourly observations per day, all from Fennoscandia and Canada. Although Alaska is represented by 11 wetland sites, sufficient winter observations with good quality are still needed. The Western Siberian Lowlands are underrepresented by EC CH₄ sites. Missing data in MODIS NBAR due to snow cover or gaps in SMAP reduced training data by 31 % and 48 % in the study domain, respectively. Filling data of MODIS NBAR to account for snow cover information and gap-filling SMAP soil moisture products can make full use of available EC observations and help improve model performance in the cold seasons. As gaps in winter SMAP data were filled with zero values, our approach has limitations with respect to the estimation of winter soil moisture gaps in areas where the zero-curtain period and talik were not represented by our interpolated soil temperatures (e.g., in coastal areas).

Many wetland sites in the study are located in areas with peatland presence, with 35 % of sites in peatland-rich areas with >50 % peatland cover (Hugelius et al., 2020). Mineral soil (soil containing less than 12 % organic carbon by weight) marshes, although covering only 5 % of the total wetland area in the northern high latitudes, need to be considered when deploying new EC sites due to their high CH₄ emissions (Kuhn et al., 2021; Olefeldt et al., 2021). This study identified regional CH₄ emission hotspots and areas undergoing strong interannual variations that are yet not part of the current FLUXNET network. However, the 10 km resolution of the RF estimates prohibits the identification of local hotspots that may occur at <1–10 m scales (Elder et al., 2021). The wall-to-wall flux maps also provide spatially continuous information to effectively further develop the CH₄ flux tower network.

4.3 Budget comparison

WetCH₄ estimated annual and seasonal mean emissions that were comparable to existing data-driven products in the study domain (Table S5). With the dynamic WAD2Mv2 map, our estimation was 0.7 Tg CH₄ yr⁻¹ smaller than UpCH₄ due to the mean seasonal cycles between 2010 and 2020 from WAD2M applied in our estimation. With the same static GLWDv1 map, our estimation was about 22 % larger than the estimate from Peltola et al. (2019) (37.5±12 Tg CH₄ yr⁻¹ for 2013–2014), despite the different periods. This is attributed to higher fluxes estimated by WetCH₄ in the DJF and JJA seasons. With two versions of the static GLWD maps, we estimated potential annual emissions of between 46.0 and 51.6 Tg CH₄ yr⁻¹. Compared to GLWDv1, version 2 of GLWD mapped smaller wetland fractions in the Hudson Bay Lowlands with intense CH₄ fluxes and more wetlands in the northwest of the Ural Mountains, Eastern Siberia, and the Sanjiang Plain, where CH₄ intensities were weaker, resulting in a larger estimate of the annual emission (Fig. S13). The wide range of data-driven estimates stemmed from the differences in wetland maps. While WAD2M provides crucial information on wetland inundation dynamics controlling interannual and inter-seasonal changes in CH₄-emitting areas, areas with saturated soil in the Arctic tundra are likely severely underestimated (Fig. 7d), requiring more accurate maps delineating wet tundra communities at higher spatial resolution (e.g., <1 km). Incorporating wetland fractions derived from high-resolution thematic maps (e.g., CALU) can improve the use of WAD2M in cold regions. Development or improvement of higher-resolution microwave remote-sensing products capable of tracking dynamic changes in local soil moisture conditions is also needed. Together, these two components likely currently yield the largest sources of uncertainty in high-latitude terrestrial CH₄ budgets.

Bottom-up estimates of wetland CH₄ emissions from data-driven GCP ensemble means and WetCHARTs are smaller than the top-down CarbonTracker-CH₄ estimate of natural microbial emissions, as the latter includes emissions from aquatic systems. Aquatic CH₄ emissions for this region have been estimated at 5.5 Tg CH₄ yr⁻¹ from rivers and streams (Rocher-Ros et al., 2023) and at 16.6 Tg CH₄ yr⁻¹ from lakes (Johnson et al., 2022). The total emission budget for wetlands and open water, based on this study and the aquatic estimates, is about 44.9 Tg CH₄ yr⁻¹, which is 4 Tg CH₄ yr⁻¹ more than the CarbonTracker-CH₄ estimate. The amplitudes of the WetCH₄ seasonal mean fluxes align with bottom-up and top-down estimates. Differences in the seasonal dynamics of wetland maps are the major source of upscaling uncertainty and result in various uncertainties between regional estimates. While atmospheric inversion models need bottom-up estimates as priors, data-driven, upscaled CH₄ products offer alternatives to process-based estimates to assist with inversion models in regions where data-driven models perform well (Bloom et al., 2017; Melton et al., 2013).

4.4 Future directions

Future development of EC networks in the northern high latitudes is urgently required to provide the additional observations needed to improve the model-based upscaling of CH₄ flux budgets and to address current gaps in the ecosystem- and regional-scale representation. Deploying new sites in underrepresented areas will not only benefit flux upscaling efforts but also our understanding of how ecosystem metabolism responds to the changing climate (Baldocchi, 2020; Pallandt et al., 2022; Villarreal and Vargas, 2021). With the availability of long-term predictor variable data, it is possible to expand upscaling frameworks over longer periods (e.g., 2000 to the present), when adequate flux observations in 2000–2010 from chambers are compiled, as 96 % of the data in FLUXNET-CH₄ were recorded after 2010 (McNicol et al., 2023).

Several data products exist for the meteorological predictor variables. Quantification of the measurement uncertainties between products of predictor variables and how the uncertainties propagate to upscaling products need to be addressed in future work. The mismatch of spatial scales between tower footprints and predictor variables may cause the underestimation of abruptly high fluxes measured in tower landscapes when environmental conditions are averaged over 0.5° grids (Chu et al., 2021; McNicol et al., 2023). Therefore, downscaling predictor variables to develop higher-resolution products is needed, especially for the Arctic region, where thermokarst development is shaping permafrost landscapes with fragments of wetlands, thermokarst ponds, and forests (Miner et al., 2022; Osterkamp et al., 2000; Wik et al., 2016). For example, Fang et al. (2022) downscaled global SMAP surface soil moisture to 1 km resolution, and optical/thermal and microwave fusion methods have been developed to downscale soil moisture (Peng et al., 2017). Nevertheless, downscaled products for root zone or profile soil moisture are needed for upscaling CH₄ fluxes, as are soil temperature products.

Beyond the ML-based upscaling framework, hybrid modeling of the data-driven approach and process-based models is a promising but challenging direction of future study (Reichstein et al., 2019). One practice constrained regional data-driven fluxes with top-down estimates via auto-learned weights on per-pixel fluxes in a region (Upton et al., 2024). Another practice pretrained a time-dependent ML algorithm with initialization from process-based synthetic data and then fine-tuned the model with observations (Liu et al., 2022). Finally, leveraging physical constraints to increase the interpretability of data-driven models and computational efficiency is still an important factor to consider in all hybrid modeling.

5 Code and data availability

The daily CH₄ flux intensities in the northern wetlands at a spatial resolution of 0.098°×0.098° and the associated uncertainties, along with daily emissions weighted by WAD2M, GIEMS2, and GLWDv1, can be accessed at https://doi.org/10.5281/zenodo.10802153 (Ying et al., 2024). The source code for the ML modeling and upscaling is publicly available at https://doi.org/10.5281/zenodo.10882613 (Ying, 2024).

6 Conclusions

We developed an ML framework (WetCH₄) to upscale daily wetland CH₄ fluxes of mid-high northern latitudes at 10 km spatial resolution that combines EC tower measurements with satellite observations and climate reanalysis. WetCH₄ is novel in that it is the first upscaling framework to introduce SMAP soil moisture and MODIS reflectance in modeling wetland CH₄ fluxes to improve accuracy (mean R²=0.70). The remote-sensing products provided high-spatial-resolution constraints associated with the abiotic controllers of CH₄ fluxes, indicating the importance of using high-spatial-resolution inputs in models for accurately simulating the spatiotemporally variable CH₄ emissions from heterogeneous northern wetland landscapes. The framework highlights the importance of soil temperature, vegetation, and soil moisture for modeling CH₄ fluxes with a data-driven approach. Using WetCH₄, an average annual CH₄ emissions of 22.8±2.4 Tg CH₄ yr⁻¹ with WAD2Mv2 was estimated and ranged between 15.7±1.8 Tg CH₄ yr⁻¹ with GIEMS2 and 51.6±2.2 Tg CH₄ yr⁻¹with GLWD v2 from vegetated wetlands (>45° N) for 2016–2022, approximately 14 %–32 % of the global wetland CH₄ budget (Saunois et al., 2020). Differences in estimates of wetland CH₄ emissions due to different wetland maps applied, highlighting the need for high-resolution wetland maps and the accurate delineation of wet soil dynamics. Emissions were relatively lower in 2017–2019 and intensified in 2016, 2020, and 2022, with the largest interannual variations coming from Western Siberia. Spatiotemporal distributions of CH₄ fluxes find emission hotspots and regions of intensified interannual variations that are not currently measured with EC. Based on a comparison with current EC sites, we suggest the need for tower observations in the wetlands of Western Siberia and Western Canada and diversified observations across wetland types. More site observations of soil-water-related variables are needed to improve the understanding of flux controls in northern wetland ecosystems. Future wetland CH₄ upscaling work could benefit from improved soil moisture products and hybrid modeling.