Ensemble estimates of global wetland methane emissions over 2000–2020

Abstract

Due to ongoing climate change, methane (CH₄) emissions from vegetated wetlands are projected to increase during the 21st century, challenging climate mitigation efforts aimed at limiting global warming. However, despite reports of rising emission trends, a comprehensive evaluation and attribution of recent changes remains limited. Here we assessed global wetland CH₄ emissions from 2000–2020 based on an ensemble of 16 process-based wetland models. Our results estimated global average wetland CH₄ emissions at 158 $\pm$ 24 (mean $\pm$ 1$\mathit{\sigma }$) Tg CH₄ yr⁻¹ over a total annual average wetland area of 8.0 $\pm$ $\mathrm{2.0}\times {\mathrm{10}}^{\mathrm{6}}$ km² for the period 2010–2020, with an average increase of 6–7 Tg CH₄ yr⁻¹ in 2010–2019 compared to the average for 2000–2009. The increases in the four latitudinal bands of 90–30° S, 30° S–30° N, 30–60° N, and 60–90° N were 0.1–0.2, 3.6–3.7, 1.8–2.4, and 0.6–0.8 Tg CH₄ yr⁻¹, respectively, over the 2 decades. The modeled CH₄ sensitivities to temperature show reasonable consistency with eddy-covariance-based measurements from 34 sites. Rising temperature was the primary driver of the increase, while precipitation and rising atmospheric CO₂ concentrations played secondary roles with high levels of uncertainty. These modeled results suggest that climate change is driving increased wetland CH₄ emissions and that direct and sustained measurements are needed to monitor developments.

1 Introduction

Wetlands are the largest single source in the global methane (CH₄) budget, representing $\sim$ 25–35 % of the total combined natural and anthropogenic sources (Kirschke et al., 2013; Saunois et al., 2016, 2020), with an uncertainty range of 100–230 Tg CH₄ yr⁻¹ (Cao et al., 1996; Gedney et al., 2004; Bousquet et al., 2006; Petrescu et al., 2010; Spahni et al., 2011; Melton et al., 2013; Bridgham et al., 2013; Bloom et al., 2017; Poulter et al., 2017). Covering 8 %–10 % of the global land surface (Zhang et al., 2021a), wetland area is sensitive to climate variations (Zhang et al., 2018; Zhu et al., 2017). Over the last deglaciation, wetlands played an important role in driving the rise of atmospheric CH₄ concentrations (Hopcroft et al., 2017; Nisbet et al., 2023; Kleinen et al., 2023). In recent decades, wetlands have experienced unprecedented and ongoing changes, including continuous thawing of permafrost (Natali et al., 2019; Treat et al., 2018), land-use change (Fluet-Chouinard et al., 2023), a lengthening of the growing season in the Arctic (Arndt et al., 2019), and expansion in tropical areas due to enhanced precipitation (Fleischmann, 2023). Recent evidence from in situ measurements (Rößger et al., 2022), data-driven estimates (Yuan et al., 2024; Ying et al., 2024), and satellite observations (Feng et al., 2022) suggests that these ongoing changes could enhance wetland CH₄ emissions and thus affect the trajectory of the atmospheric CH₄ concentration. Furthermore, atmospheric ${\mathit{\delta }}^{\mathrm{13}}$C-CH₄ records also show a trend toward increased depletion since the late 2000s (Lan et al., 2021; Nisbet et al., 2019), indicating that isotopically light biogenic sources, such as wetlands (Basu et al., 2022; Feng et al., 2022), agriculture, and waste (Schaefer et al., 2016; Zhang et al., 2021b), have become dominant contributors to the rise in atmospheric CH₄. Current estimates of wetland CH₄ emissions (hereafter denoted as eCH₄) in response to climate change are projected to increase by up to 15 %–30 % by 2050 (Koffi et al., 2020; Zhang et al., 2017), accounting for 25 %–40 % of the pledged reduction in anthropogenic emissions (Shindell et al., 2019). These trends and projections suggest that the emerging wetland CH₄ climate feedback that influences atmospheric CH₄ concentration requires a better understanding of long-term changes in eCH₄.

Directly diagnosing the variations and trends of eCH₄ at large scales is challenging. Site-level measurements, such as those from chamber and eddy covariance techniques, are useful for identifying underlying mechanisms and monitoring CH₄ fluxes at the landscape scale but are difficult to upscale due to large uncertainties in extrapolation and the high spatial heterogeneity of wetland CH₄ fluxes (Chu et al., 2021; Kuhn et al., 2021). Interpreting eCH₄ using satellite observations and inversions of atmospheric concentration data is also subject to uncertainties in anthropogenic sources, other natural sources, atmospheric chemistry, and model errors associated with atmospheric transport (Gatti et al., 2021; Gloor et al., 2021; Palmer et al., 2022; Patra et al., 2011; Zhang et al., 2021c). Global wetland models, integrated within land biosphere models, can serve to bridge our understanding of wetland CH₄ processes and diagnosing wetland CH₄ dynamics at large scales (Melton et al., 2013; Wania et al., 2013). These models provide mechanistic explanations for the causes of changes in eCH₄ dynamics. Furthermore, recent advances in wetland models (Arora et al., 2018; Kaiser et al., 2017; Shu et al., 2020; Grant, 2017; Chang et al., 2020) show significant potential for improving our understanding of eCH₄ through the incorporation of complex biogeochemical processes.

Current studies have reached various conclusions on the change in eCH₄ over the last decades. Studies based on single biogeochemical models (Zhang et al., 2018; Zhu et al., 2017) suggest a significant increase in eCH₄ from 2000–2006 to 2007–2017, while atmospheric inversions (Zhang et al., 2021c; Yin et al., 2021; Basu et al., 2022; Feng et al., 2022) suggested even higher rate increases from 2 to 3 Tg CH₄ yr⁻¹ yr⁻¹ during the post-2010 period. Poulter et al. (2017) reported no significant change between the 2000–2006 and 2007–2012 periods based on an ensemble of wetland models, while Saunois et al. (2020) show a slight increase ($\sim$ 2 Tg CH₄ yr⁻¹) on average for 2007–2017 compared to the 2000–2006 level using a large set of wetland CH₄ models. However, these models demonstrate considerable differences in estimated eCH₄, both spatially and temporally (Ma et al., 2021; Parker et al., 2020; Chang et al., 2023), primarily due to the sensitivity of their estimations to the wetland areal extent, the implemented biogeochemical structures, and parameterizations. The multi-model ensemble approach is applied to increase the skill, reliability, and consistency of model forecasts, potentially offsetting individual model errors (Schaefer et al., 2012). However, a recent study (Chang et al., 2023) found that downsampling atmospheric inversion and wetland model CH₄ predictions based on a comparison to eddy covariance data did not reduce uncertainty in global eCH₄ estimates. Therefore, it has become necessary to thoroughly evaluate the performance of these models using the most recent generation of wetland models against an increasingly dense network of observations (Delwiche et al., 2021; Knox et al., 2019) from eddy covariance sites.

Here we conducted ensemble simulations of 16 wetland biogeochemical models following a common modeling protocol to provide monthly integrated global eCH₄ for the period of 2000–2020, as part of the Global Carbon Project's methane budget activity. The inundation dynamics of each model were simulated using a model-specific prognostic hydrological modeling approach as well as a set of diagnostic satellite-driven simulations. A set of factorial simulations were carried out to isolate the effects of temperature, precipitation, and rising atmospheric CO₂ concentration. The modeled temperature sensitivity was evaluated against the global eddy covariance database FLUXNET-CH₄ (Delwiche et al., 2021; Knox et al., 2019) and a data-driven global wetland CH₄ upscaling dataset UpCH₄ (McNicol et al., 2023) based on FLUXNET-CH₄. In addition, we examined the changes in eCH₄ for the year 2020, which was characterized as an extremely warm and wet year with the highest growth rate of atmospheric CH₄ observed over the study period.

2 Methods

2.1 Wetland model ensemble

Sixteen wetland models participated in the ensemble simulations (Table S1). Wetland CH₄ models can be generally described as functions describing the biogeochemical processes that control CH₄ production and oxidation through methanogenesis and methanotrophy, as well as the biophysical processes that regulate CH₄ transport from the soil to the atmosphere (Table S1). Methanogenesis in the models is linked to different proxies (e.g., carbon substrate, heterotrophic respiration, net primary production) with a wide range of model complexity – more sophisticated models include wetland plant functional types (PFTs) and explicitly simulate the processes of CH₄ production, consumption, and transport, while the simplified models use generalized empirical equations to simulate net fluxes without explicitly calculating individual components of the CH₄ flux.

Wetlands were defined as naturally vegetated forested and non-forested ecosystems with saturated/inundated areas, excluding coastal wetlands; cultivated wetlands such as rice paddies; and open water systems such as rivers, lakes, ponds, and reservoirs. A prognostic wetland inundation scheme and a diagnostic wetland dataset (Wetland Area and Dynamics for Methane Modeling (WAD2M v2); Zhang et al., 2021a) were applied to identify the wetland areal dynamics. The prognostic wetland areal dynamics were independently determined by each model's hydrological modules, which use water table depth or soil moisture, combined with subgrid topographic conditions, to determine saturated areas within a land surface grid cell (Zhang et al., 2016; Xi et al., 2022). Among the participating models, there was a large variation in complexity and in the level of comprehensiveness with which wetland extents were characterized. The modules for simulating inundation ranged from simplified TOPMODEL approaches to more sophisticated representations of water table variation, with the estimated magnitude being influenced by the hydrologic schemes utilized and the sensitivities to precipitation. The prognostic modeled wetland extent showed large variability in estimated magnitude but was consistent with satellite-based inundation products in predicting different phases of inundation (Xi et al., 2022; Zhang, et al., 2021a). The ensemble mean of the modeled wetland extent is close to $\mathrm{7.5}\times {\mathrm{10}}^{\mathrm{6}}$ km² as estimated by WAD2M but higher than the $\mathrm{4.6}\times {\mathrm{10}}^{\mathrm{6}}$ km² by the satellite-based product Global Surface Water Extent and Dynamics version 2 (GIEMS-2; Prigent et al., 2020). The modeled temporal variations in wetland areas show high correlations with satellite-based products for temperate regions and high latitudes (Fig. S1), except in the tropics. The limited agreement in the tropics may be due to the influence of aerosols and clouds on satellite-based measurements, as well as the process-based model's performance limitations in representing wetland areas. The diagnostic runs are exclusively used for temperature dependence calculations due to a discontinuity issue in WAD2Mv2 over a few tropical hotspots, which affect a subset of wetland models that are particularly sensitive to inundation in the hotspots.

2.2 Modeling protocol and simulation setups

The modeling protocol aimed to provide wetland CH₄ fluxes and quantify the associated uncertainties arising from model differences, meteorological forcing, and wetland extent dynamics. To quantify meteorological forcing uncertainty, we used two climate inputs: a ground-based monthly climate dataset from the Climatic Research Unit (CRU) (Harris et al., 2014) up to 2020 and a harmonized daily dataset from the Global Soil Wetness Project Phase 3 (GSWP3-W5E5) through the year 2019, which is a multiple-source-based daily dataset (Cucchi et al., 2020; Dirmeyer et al., 2006) used in the Inter-Sectoral Impact Model Intercomparison Project 3a (ISIMIP3a). For models that require 6-hourly meteorological forcings, a temporal-interpolation dataset (CRU-JRA) was applied based on the Japanese Reanalysis Agency (JRA55) aligned with CRU. The atmospheric CO₂ concentration values for 1861–2020 were obtained from the CMIP6 experimental protocol (Meinshausen et al., 2017). Ancillary data, such as soil texture and CH₄-related parameter sets, used model-specific inputs. All the models were run in “natural vegetation” mode without the transient effects of land-use or land-cover change. Methane oxidation in wetland soils was implicitly included in the estimate, but the upland oxidative sink was not included as it was not part of the net wetland emission calculations. Models included the spin-up period to pre-industrial conditions, assuming net ecosystem exchange equilibrium before 1860, by recycling fixed CO₂ concentrations (1860 level of 286.42 ppm) and meteorology (1901–1920).

2.3 FLUXNET-CH₄ and the machine learning-based upscaling product UpCH₄

FLUXNET-CH₄ is the first global dataset of CH₄ eddy covariance measurements that includes $\sim$ 80 sites globally, including different wetland types such as peatlands (e.g., bog, fen), mineral wetlands (e.g., marsh, swamp), and rice paddies. For this study, a subset of natural freshwater wetland sites was selected for the analysis. All the eddy covariance measurements used in this study were gap-filled daily and total fluxes were filled using an artificial neural network (ANN) approach (Knox et al., 2019). In addition, a data-driven gridded dataset (UpCH₄; McNicol et al., 2023) for 2001–2018, which is based on 119 site-years of CH₄ fluxes from the FLUXNET-CH₄ dataset, was applied in the comparison. This dataset used a random forest model to upscale ground-based eddy covariance CH₄ flux data and was then forced with globally gridded predictor data and two wetland extent products to predict wetland CH₄ emissions. The predictors included data sources from climate, biometeorological, and soil properties.

2.4 Time series decomposition and statistical analyses

To attribute the time series of global wetland CH₄ emissions to what we consider the dominant drivers of change (i.e., temperature, precipitation, and CO₂ concentration), we applied a multiple regression approach (Piao et al., 2013) to estimate the parameters of global wetland CH₄ sensitivity to climate drivers using the following equation:

1$$y=\mathit{\beta }{\mathrm{CO}}_{\mathrm{2}}+\mathit{\gamma }\mathrm{Tmp}+\mathit{\delta }\mathrm{Pre}+c+\mathit{\epsilon },$$ where $y$ is the global annual total wetland CH₄ emission of each model from the transient run or from the observation-based upscaling dataset UpCH₄, and Tmp, Pre, and CO₂ are the mean annual temperature, total annual precipitation, and mean atmospheric CO₂ concentration for that year, respectively. $\mathit{\gamma }$, $\mathit{\delta }$, $\mathit{\beta }$, and $c$ are regression coefficients, and $\mathit{\epsilon }$ is the residual error term. The regression coefficients were calculated using a maximum likelihood estimate. Changes in other meteorological forcings may also influence the estimation of $y$. These confounding drivers, such as solar radiation and wind speed, although they are considered to have minor impacts on the variations of eCH₄, were implicitly accounted for in the regression coefficients.

2.5 Model factorial experiment

To further separate the contribution of different controls on the change in methane emissions ($\mathrm{\Delta }$eCH₄) by climate variations and rising CO₂, we used a subset of four models that conducted factorial experimental simulations by holding each factor constant during part of the transient runs. This subset of the wetland models (i.e., four wetland models: ELM-ECA, LPJ-wsl, SDGVM, and VISIT) performed a set of factorial simulations to specifically attribute the effect of temperature, precipitation, and rising CO₂ concentration on wetland CH₄ fluxes with the climatology of 2000–2006 for 2007–2020. The simulations were performed by running the model while keeping one factor constant at a time to estimate the contribution of each component to the total range of variations (Table S2). For these factorial simulations, we evaluated the annual amplitude of wetland eCH₄ as a relative percentage change to minimize the impacts of different modeling implementation choices, such as different input variables among models. The effect of the total changes on the relative change in amplitude was represented by the difference between the transient (one factor is time-varying) and baseline (static at 2000–2006 levels) runs. For simplicity, the relative contribution of a single driver to eCH₄ variations was quantified as the transient run minus the specific control run. To calculate the contribution of each driver using the subset of the models, we calculated weighting factors per year across the models, with lower bias resulting in higher weight relative to the full ensemble mean using an inverse function.

2.6 Temperature dependence calculation

To further evaluate the response of eCH₄ to rising temperatures, we calculated the modeled seasonal eCH₄ temperature dependence, referred to as the apparent $Q_{\mathrm{10}}$ metric at the locations of 34 FLUXNET-CH₄ sites. This seasonal $Q_{\mathrm{10}}$ differs from the intrinsic $Q_{\mathrm{10}}$ prescribed in the parameterization of respiratory processes in each model. Here it represents the overall response of eCH₄ along geographic temperature gradients. The apparent $Q_{\mathrm{10}}$ is defined as eCH₄ sensitivity to temperature change. We calculated apparent $Q_{\mathrm{10}}$ based on CH₄ emitting strength over a standard wetland area, which was calculated as the CH₄ fluxes divided by inundated area on a per-pixel basis to exclude the effect of inundation dynamics. To derive the temperature dependence of eCH₄ at the soil or ecosystem level, we applied the following equation:

2$$R\left(i\right)=R_{\mathrm{b}}\left(i\right)Q_{\mathrm{10}}^{\frac{T\left(i\right)-T_{\mathrm{ref}}}{\mathrm{\Gamma }}},$$ where $R\left(i\right)$ is the net wetland flux at the location of site $i$, $R_{\mathrm{b}}\left(i\right)$ is the basal net CH₄ flux at the reference temperature $T_{\mathrm{ref}}$, and $T\left(i\right)$ is ambient temperature. The parameters $Q_{\mathrm{10}}$, $\mathrm{\Gamma }=\mathrm{10}$°, and $T_{\mathrm{ref}}=\mathrm{15}\mathit{°}$ are all time-independent constants. The $Q_{\mathrm{10}}$ acting on specific timescales can be obtained from eCH₄ at corresponding specific timescales (i.e., seasonal total and annual total) by fitting an exponential regression with modeled eCH₄ and air temperature from CRU or GSWP3-W5E5. To quantify the uncertainty in observed apparent $Q_{\mathrm{10}}$, we employed 1000 sets of resampled FLUXNET-CH₄ observations generated based on a Gaussian distribution. The uncertainty range in measured seasonal mean CH₄ fluxes was determined by aggregating the uncertainty of daily total fluxes obtained through ANN gap filling.

3 Results and discussion

3.1 Changes in eCH₄ during the period 2000–2020

The multi-model ensemble based on the prognostic inundation schemes shows that the average annual global eCH₄ over the period 2000–2020 was 156 $\pm$ 24 Tg CH₄ yr⁻¹ (mean $\pm$ 1$\mathit{\sigma }$). The average annual eCH₄ increased from 153 $\pm$ 23 Tg CH₄ yr⁻¹ during 2000–2009 to 158 $\pm$ 24 Tg CH₄ yr⁻¹ during 2010–2020. And 15 out of 22 model simulations show significant positive linear trends ($p\mathit{<}\mathrm{0.01}$), with an ensemble mean increase rate of 0.6 $\pm$ 0.3 Tg CH₄ yr⁻¹ yr⁻¹ over 2000–2020 (Fig. 1a; Table 1; Fig. S2). Differences in total annual emissions between the two sets of simulations driven by two different climate datasets (CRU and GSWP3-W5E5) agree well in the magnitude of the annual anomalies. Notable eCH₄ variations to climate events were observed, such as the rise during the 2010 La Niña ($+$5.2 Tg CH₄ yr⁻¹) and the decline during the 2015 El Niño ($-$4.6 Tg CH₄ yr⁻¹) after removing the positive linear trends. The multi-model ensemble wetland eCH₄ response to climate events is consistent with those reported by earlier studies (Zhang et al., 2018; Zhu et al., 2017) using single wetland models, indicating a modulation of the phase of eCH₄ anomaly ($\mathrm{\Delta }$eCH₄) by the El Niño–Southern Oscillation. The model ensemble demonstrates a consistent increase in interannual variability (IAV) in $\mathrm{\Delta }$eCH₄ from 3.6 $\pm$ 1.6 Tg CH₄ yr⁻¹ during 2000–2009 to 4.7 $\pm$ 1.5 Tg CH₄ yr⁻¹ during 2010–2020, suggesting a potential increase in eCH₄ variability under climate change.

The models consistently show that 2020 is the strongest positive anomaly year during 2000–2020, with a net increase of 2 [$-$2, 7] Tg CH₄ yr⁻¹ (mean [min, max]) in 2020 compared to 2019. This positive anomaly in 2020 (Table 1) is broadly consistent with a recent study (Peng et al., 2022) that reported 6 $\pm$ 2.3 Tg CH₄ yr⁻¹ based on simulations of two bottom-up models with different climate datasets. The discrepancy in estimated magnitude between the Peng et al. (2022) and our results is partly due to the parameterizations of the CH₄ module that causes lower annual magnitude in this study ($\sim$ 162 $\pm$ 23 Tg CH₄ yr⁻¹ in 2020) compared to the Peng et al. (2022) study (177 $\pm$ 31 Tg CH₄ yr⁻¹ in 2020). Additionally, the precipitation inputs in the climate forcing used in this study show a lower positive anomaly ($\sim$ of 20 mm yr⁻¹ in CRU over global wetland) in precipitation in 2020 compared to the reanalysis-based estimates ($\sim$ 40–117 mm yr⁻¹) over global wetlands used in the study by Peng et al. (2022), which leads to lower estimates of wetland area and consequently lower emissions in this study. Moreover, our model ensemble does not indicate a strong increase ($-$0.2 [$-$1.5–0.7] Tg CH₄ yr⁻¹) in eCH₄ in Africa in 2020. This contrasts with recent atmospheric inversions (Feng et al., 2023; Qu et al., 2022), which suggest a large increase of 11–17 Tg CH₄ yr⁻¹ above 2019 levels in African CH₄ emissions for 2020. The estimated increase from these inversions is equivalent to 55 %–85 % of total eCH₄ in Africa during 2010–2019 in our study (Fig. 2). These discrepancies highlight the need for further studies to investigate the differences between these two approaches, including uncertainty in climate inputs in process-based bottom-up models and partitioning different sources in atmospheric inversions.

Figure 1

Simulated global wetland CH₄ emissions from the model ensemble for 2000–2020. (a) Time series of annual total emissions during 2000–2020, with the shaded area representing the range between minimum and maximum modeled emissions. The horizontal lines represent the ensemble means of 2000–2009 (152 Tg CH₄ yr⁻¹) and 2010–2019 (158 Tg CH₄ yr⁻¹), respectively. (b) Latitudinal gradient of eCH₄ difference ($\mathrm{\Delta }$eCH₄), with the mean annual total $\mathrm{\Delta }$eCH₄ for each of the 30° latitude bins from the two sets of simulations shown. The change is calculated relative to the mean of the 2000–2009 level from the two sets of simulations with prognostic wetland emission models grouped by different climate datasets, CRU and GSWP3-W5E5. (c) Boxplots of mean seasonal $\mathrm{\Delta }$eCH₄ for the three regions. The central mark and the bottom and top edges of the box indicate the median and the 25th and 75th percentiles of the ensemble, respectively. The colored lines represent the average seasonal cycle of 2000–2009 from the simulations grouped by two climate datasets, CRU and GSWP3-W5E5.

[Figure omitted. See PDF]

There were widespread net increases in eCH₄ across all latitudinal bands during 2010–2020, compared to the average of 2000–2009, with the largest magnitudes occurring in the 90° S–30° N bands (there are relatively few wetlands in the southern extratropics of 90–30° S, contributing 0.1–0.2 Tg CH₄ yr⁻¹) and temperate regions (30–60° N) (Fig. 1b). The annual magnitude of eCH₄ increased by 3.7–3.8, 1.8–2.4, and 0.6–0.8 Tg CH₄ yr⁻¹ in the tropical, temperate, and Arctic wetlands, respectively. The tropics have experienced the largest increases in annual total emissions with an increase of 3 % relative to 2000–2009 (Table 1). This finding is aligned with the results of several recent atmospheric inversions (Basu et al., 2022; Feng et al., 2022; Lan et al., 2021) using satellite observations and/or isotopic measurements that suggest a large increase in microbial emissions for the post-2007 period in the tropics. While the increase in annual total emissions from temperate wetlands is lower than that from the tropics, they nevertheless show a larger relative increase of 5 %–8 % compared to 2000–2009. Arctic wetlands also show an increased rate of 5 %–7 % relative to the same period.

The increase in eCH₄ occurs in parallel with differing patterns of enhanced seasonal cycles between tropical and extratropical wetlands (30–90° N) (Fig. 1c). In temperate and Arctic wetlands, the majority of the increase in emissions (60 %–92 %) occurred primarily during the growing season (May–October). Specifically, increases in Arctic wetlands occurred during the early growing season (May–July), aligning with findings from a data-driven estimate (Yuan et al., 2024) and a long-term eddy-covariance-based study (Rößger et al., 2022) that observed early-growing-season increases in eCH₄ due to continuous warming in a Siberian wetland. In contrast, the increase in emissions within the 90° S–30° N band exhibited relatively minor seasonal variations throughout the year, with the May–October period accounting for a 24 % greater increase in $\mathrm{\Delta }$eCH₄ compared to the November–April period (Fig. S3).

Table 1

Summary of wetland CH₄ emissions (Tg CH₄ yr⁻¹) over different time periods by latitudinal bands for the prognostic wetland simulations. The ensemble mean with minimum and maximum (numbers within square brackets) are listed, respectively.

Time period	Forcing	90–30° S	30° S–30° N	30–60° N	60–90° N	Global
2000–2009	CRU	3[1–5]	107[63–141]	31[16–60]	11[4–29]	152[119–187]
	GSWP3-W5E5	3[1–5]	106[60–142]	33[18–57]	11[4–29]	153[116–188]
2010–2019	CRU	3[1–6]	110[67–144]	34[17–64]	12[4–30]	158[126–193]
	GSWP3-W5E5	3[1–6]	110[64–146]	35[18–60]	12[4–29]	158[118–203]

3.2 Spatial distribution of eCH₄

A few key regions contribute significantly to global emissions (Fig. 2a, c). These regions are mainly floodplains located along major river basins such as the Amazon, Ganges, Mississippi, and Yangtze; tropical peatlands in the Congo and Southeastern Asia; and high-latitude peatlands in the Hudson Bay Lowland (HBL) and West Siberian Lowland (WSL). However, inter-model variabilities in eCH₄ reveal varying levels of spatial agreement between models, with the largest discrepancies coming from South America and Africa. South America is one of the largest contributors to the global total eCH₄. Still, the net change in that region shows only a moderate increase, with diverging trends within the Amazon Basin during the 2010s (Fig. 2b, d). The uncertain temporal trends are consistent with a long-term, large-scale atmospheric inversion based on airborne campaigns (Basso et al., 2021). South Asia and Africa are among the regions with the largest increases in the tropics, next to North America, but have high uncertainty with a lower level of agreement among the models (Fig. S4). The model ensemble shows that northwestern South Asia has a significant percentage increase in eCH₄ during 2010–2019 relative to its average levels from 2000–2009, suggesting a possible high sensitivity of eCH₄ to climate change in this region.

The comparison with previous estimates from bottom-up approaches and top-down atmospheric inversions (Table S3) suggests that the model ensemble mean generally captures well the spatial distribution of annual eCH₄, with a potential underestimation for a few methane hotspots (Fig. S5). The model ensemble means for the Amazon Basin, HBL, and WSL show good agreement with atmospheric inversions (Bergamaschi et al., 2013; Pickett-Heaps et al., 2011; Ringeval et al., 2014; Tunnicliffe et al., 2020; Wilson et al., 2016, 2021) and bottom-up modeling estimates (Bansal et al., 2023; Bloom et al., 2017; Bohn et al., 2015), with relatively low uncertainty. The model ensemble highlights WSL and HBL as CH₄ hotspots in the high latitudes, with good agreements of annual magnitudes with atmospheric inversions and in situ observations (Bohn et al., 2015; Glagolev et al., 2011; Pickett-Heaps et al., 2011), while the models have lower estimates for Alaska compared to the inversions (Chang et al., 2014; Miller et al., 2016). However, for the two hotspots of the Pantanal and Sudd wetlands, the models tended to underestimate the annual eCH₄ compared to a few recent satellite-based estimates (Gerlein-Safdi et al., 2021; Gloor et al., 2021; Lunt et al., 2021; Pandey et al., 2021), with a large uncertainty range of up to 2 orders of magnitude across the model ensemble (Fig. S5). In addition to the regions where eCH₄ are being underestimated, recent studies (France et al., 2022; Shaw et al., 2022) based on aircraft measurements suggest that the bottom-up models likely underestimate high eCH₄ fluxes in some little-studied wetlands, such as those in Zambia and Bolivia. The underestimations by process-based wetland models can be attributed to (1) the challenge in accurately capturing the areal dynamics of wetlands under varying hydrological conditions, such as in flat terrains that receive lateral transport of water from upper streams (Li et al., 2024; Lunt et al., 2021; Gerlein-Safdi et al., 2021); (2) existing knowledge gaps in mapping wetlands in remote areas, which affect the parameterization of inundation modeling; and (3) the limited representation of water table regulation (Chen et al., 2021) and wetland PFTs (Bastviken et al., 2023) on eCH₄ in biogeochemical models.

Figure 2

Spatial distribution of eCH₄ and the average change between the 2010s and 2000s. (a) Map of mean eCH₄ (unit: g CH₄ m⁻² yr⁻¹ per 0.5° grid cell) for 2000–2020. The regions defined in panels (c) and (d) and regional CH₄ hotspots in Table S3 are outlined in black and red, respectively. (b) Map of change in mean annual wetland emissions ($\mathrm{\Delta }$eCH₄) between the 2010s and 2000s. (c) Boxplot of mean annual eCH₄ and (d) $\mathrm{\Delta }$eCH₄ by regions for 2000–2020 in ascending order for median estimates. Afr: Africa; CAs: Central Asia; EAs: East Asia; Eur: Europe; NAm: North America; NAs: North Asia; Oz: Oceania; SAm: South America; SAs: South Asia; SEAs: Southeast Asia.

[Figure omitted. See PDF]

Figure 3

Attributions of $\mathrm{\Delta }$eCH₄ during 2000–2020. (a) Histogram showing the sensitivity coefficients derived from a multiple regression approach (see the Methods section) for temperature ($\mathit{\gamma }$), precipitation ($\mathit{\delta }$), and atmospheric CO₂ concentration ($\mathit{\beta }$). The curves represent probability distributions of sensitivity coefficients across the models, assuming a Gaussian distribution. Vertical lines represent estimates from the machine learning-based dataset UpCH₄, with different colors corresponding to different climate datasets. (b) Time series of anomalies for annual mean temperature ($\mathrm{\Delta }T$), annual total precipitation ($\mathrm{\Delta }P$), and annual mean wetland extent ($\mathrm{\Delta }$Fw) for 2000–2020 for CRU and 2000–2019 for GSWP3. The shaded areas in $\mathrm{\Delta }$Fw represent the minimum and maximum ranges from the prognostic model simulations. Dashed lines are linear fitted trends for the corresponding variables.

[Figure omitted. See PDF]

3.3 Attribution of wetland CH₄ changes

To evaluate the relative contribution of different factors on global eCH₄, we calculated the sensitivity of eCH₄ to mean annual temperature (denoted as $\mathit{\gamma }$), annual total precipitation (denoted as $\mathit{\delta }$), and CO₂ concentration (denoted as $\mathit{\beta }$) using a multiple regression approach for each model run over the period of 2000–2020. The same approach was applied to the upscaled gridded machine learning dataset UpCH₄, which uses eddy covariance measurements from FLUXNET-CH₄ as training inputs. The model ensemble suggests that temperature is the primary driver of the increase in eCH₄ (Fig. 3a). The regression coefficient for $\mathit{\gamma }$ is 4.6 Tg CH₄ yr⁻¹ °C⁻¹, with a range of $-$0.4 and 9.0 Tg CH₄ yr⁻¹ °C⁻¹ between the 10th and 90th percentiles among all models. This mean temperature sensitivity is slightly higher than the $\mathit{\gamma }$ coefficient of 3.2–4.1 Tg CH₄ yr⁻¹ °C⁻¹ estimated for UpCH₄. In contrast, precipitation contributed little to the increase from the prognostic simulations, with a coefficient $\mathit{\delta }$ of 0–0.3 Tg CH₄ yr⁻¹ mm⁻¹. The coefficient $\mathit{\delta }$ was lower at $-$0.05–0 Tg CH₄ yr⁻¹ mm⁻¹ for UpCH₄, as precipitation was not chosen as a model training predictor through its feature selection, based on site-level eddy covariance measurements (McNicol et al., 2023). However, precipitation is a more dominant factor at large scales, especially for tropical floodplains, which contribute the largest proportion of emissions but are poorly represented by eddy covariance measurements. The model-ensemble-estimated $\mathit{\beta }$ remains small, ranging from 0 to 0.3 Tg CH₄ yr⁻¹ ppm⁻¹, while UpCH₄ suggests a $\mathit{\beta }$ at $-$0.01 Tg CH₄ yr⁻¹ ppm⁻¹. However, other confounding drivers might influence eCH₄ as well, such as solar radiation, wind speed, and nitrogen deposition. Thus, the inferred sensitivities are implicitly accounted for in the regression coefficients despite their relatively small impacts compared to the major drivers.

Generally, the factorial simulations of the four-model subset indicated a consistently positive contribution (three out of four) from rising temperature to $\mathrm{\Delta }$eCH₄, with a large variability (SD $=$ 4.3 Tg CH₄ yr⁻¹) of contributions from precipitation (Fig. S6). The strength of the CO₂ fertilization effect varied among models and was moderate but positive in all models. Two models (ELM-ECA and SDGVM) were among the models with higher sensitivity to climate variations, while LPJ-wsl and VISIT were close to the full ensemble mean. ELM-ECA produced a negative temperature effect on eCH₄, likely due to its modeled nutrient constraints and higher temperature sensitivity for methanotrophic compared to methanogenic processes. Considering the deviation of each model from the full ensemble mean, the weighted mean (Fig. S7) contributions for temperature, precipitation, and CO₂ concentration from the subset models were 3.2, 1.8, and 1.4 Tg CH₄ yr⁻¹, respectively. The results from the subset of the models consistently demonstrate that temperature is the primary factor influencing eCH₄.

Overall, the interannual variations of modeled eCH₄ were primarily associated with rising temperature, altered precipitation patterns, and rising atmospheric CO₂ concentrations that stimulated ecosystem productivity through the CO₂ fertilization effect (Yvon-Durocher et al., 2014). We note that a recent study found strong hysteresis in the seasonal temperature dependence of observed eCH₄ using the FLUXNET-CH₄ dataset (Chang et al., 2021). Those hysteretic features likely result in uncertainty in annual temperature sensitivity estimates but would not bias the conclusion of temperature as a dominant controller of eCH₄ at the decadal timescale. The links between rising temperature and enhanced net CH₄ fluxes are evident (as described below), as the annual global average temperature over wetland areas has significantly ($p\mathit{<}\mathrm{0.01}$) increased by 0.5–0.7 °C from 2000–2020 (Fig. 3b). The modeled interannual variations of wetland extent dynamics reproduced the response to strong climate events (e.g., positive anomaly during the La Niña phase in 2010/2011 (Boening et al., 2012) and 2020). Both climate-forcing datasets suggest no significant trend in the anomaly of annual mean wetland area globally over the same period based on the prognostic hydrological simulations (Fig. 3b). Similarly, no significant regional trends in wetland area were found for most of the subregions, with the exception of South America, which shows a decrease, and East Asia, which shows a slight increase (Fig. S8). Considering that the extent of modeled wetland areas is primarily driven by precipitation, we do not detect a substantial contribution of changes in wetland extent to the long-term increase in eCH₄ over 2000–2020 based on the climate datasets. However, considerable differences in annual and seasonal precipitation estimates between the climate datasets used in this study and those derived from reanalysis or satellite-based products (Zhang et al., 2023a) result in large uncertainties in the estimated trends in wetland extent.

Figure 4

Temperature dependence of simulated seasonal eCH₄ across locations of FLUXNET-CH₄ sites. (a) Model ensemble mean (“Model Ensmean”) of simulated eCH₄ against seasonal mean temperature for the JJA season along the temperature gradient at the locations of FLUXNET-CH₄ sites in comparison to the estimates from eddy covariance measurements (“Obs”; Fig. S10; Table S4) and UpCH₄. Each dot represents the value at one site for an individual year when observations are available. The unit of the simulated CH₄ emissions is gCH₄ m⁻¹ month⁻¹ per standard wetland area to exclude the effect of inundation on eCH₄. The exponential fitted curves are shown. (b) Histogram of the seasonal $Q_{\mathrm{10}}$ for the 16 individual models for the months of DJF, MAM, JJA, and SON. Sample sizes are shown in the plot. The $Q_{\mathrm{10}}$ values derived from FLUXNET-CH₄, UpCH₄, and the model ensemble mean are shown as vertical solid lines, with a width of the bar for “Obs” indicating the uncertainty range of $Q_{\mathrm{10}}$ based on measurement uncertainty.

[Figure omitted. See PDF]

3.4 Temperature dependence of wetland CH₄ models

The modeled CH₄ emissions show an exponential relationship between eCH₄ and air temperature, with higher temperatures corresponding to higher mean eCH₄ during the peak growing season (JJA, June–July–August) in the Northern Hemisphere (Fig. 4a). The model ensemble mean of eCH₄ response to temperature shows good agreement within the range of the spread when compared to the site-level measurements from FLUXNET-CH₄ and the gridded product UpCH₄. The model ensemble mean has a higher CH₄ emitting strength (i.e., CH₄ emission per standard wetland area) for the high latitudes, leading to lower apparent $Q_{\mathrm{10}}$. This implies that the model-ensemble-estimated temperature dependence for the high latitudes could be potentially overestimated during the JJA season. The apparent $Q_{\mathrm{10}}$ values for individual models show a large spread (Fig. S9), with 11 out of the 16 models having statistically significant ($p\mathit{<}\mathrm{0.01}$) exponential relationships. The good agreement between the ensemble mean and observations suggests that the ensemble approach provides a better constraint compared to single models alone. Furthermore, it is important to acknowledge that the sparse spatial coverage of FLUXNET-CH₄ over low latitudes, especially for underrepresented areas such as Africa, Southeast Asia, and South America, limits our ability to evaluate temperature dependencies over high-temperature regions (Fig. S10).

The modeled apparent $Q_{\mathrm{10}}$ exhibits an average temperature dependence similar to that of ecosystem respiration, as reported by previous studies (Bloom et al., 2017; Mahecha et al., 2010; Yvon-Durocher et al., 2014), indicating that the underlying factors controlling the response of eCH₄ and ecosystem respiration to temperature covary. The modeled temperature dependences are more constrained with less spread for JJA and SON (September–October–November) than DJF (December–January–February) and MAM (March–April–May) when most site-level measurements have limited availability. The seasonal variations of modeled apparent $Q_{\mathrm{10}}$ differ from site-level observations or UpCH₄, reflecting discrepancies in the involved processes between eddy covariance and land surface models. Given that underrepresented processes such as substrate supply tend to have higher sensitivity of ecosystem metabolic processes to temperature, it is likely that the models do not entirely capture the fine-scale processes that affect the overall temperature response (Chang et al., 2021). In addition, the absence or underrepresentation of certain biophysical processes could lead to lower modeled apparent $Q_{\mathrm{10}}$. For instance, the ensemble mean of modeled apparent $Q_{\mathrm{10}}$ for SON months is underestimated, likely linked to the limited representation of processes during the freeze–thaw cycle (e.g., zero-curtain period), as suggested by previous observational studies (Mastepanov et al., 2008; Zona et al., 2016).

4 Conclusions

Our results estimated global average wetland CH₄ emissions at 158 $\pm$ 24 (mean $\pm$ 1$\mathit{\sigma }$) Tg CH₄ yr⁻¹ for the period 2010–2020, with an average decadal increase of 6–7 Tg CH₄ yr⁻¹ compared to the decade of 2000–2009. The increases in the four latitudinal bands of 90–30° S, 30° S–30° N, 30–60° N, and 60–90° N were 0.1–0.2, 3.6–3.7, 1.8–2.4, and 0.6–0.8 Tg CH₄ yr⁻¹, respectively, during the 2 decades. Our analysis reveals how global wetlands respond to variations in the primary climatic controls of temperature, precipitation, and rising CO₂ concentrations. The model average shows good agreement with eddy covariance measurements on temperature dependence, confirming the primary role of temperature in the rising trajectory of eCH₄ at decadal timescales. Furthermore, the modeled ensembles of prognostic wetland extents offer a complementary approach to satellite-based estimates (Prigent et al., 2020; Zhang, et al., 2021a) and enable further investigation into the uncertainties in wetland area estimation. These differences can motivate improvements to inundation schemes through an improved water table position (Chen et al., 2021) and lateral flow representation. Note that a large portion of tropical wetlands comprise inundated floodplains connecting rivers, where the leaching of methane production from wetlands to river networks is not accounted for in the wetland models. The prognostic models estimate an annual mean maximum wetland area of 8.0 $\pm$ $\mathrm{2.0}\times {\mathrm{10}}^{\mathrm{6}}$ km², with a seasonal cycle (annual maximum minus annual minimum) of 4.7 $\pm$ $\mathrm{2.0}\times {\mathrm{10}}^{\mathrm{6}}$ km². Resolving the large uncertainty in wetland areas and seasonal variation remains a high priority to refine bottom-up estimates of eCH₄. Lastly, our results highlight the important but highly uncertain CO₂ fertilization effect on eCH₄. The mean sensitivity coefficient $\mathit{\beta }$ and results from the factorial experiment suggest a net increase in eCH₄ of 0.1 %–2.3 % relative to the annual total under an average $\sim$ 20 ppm increase in atmospheric CO₂ concentration. In comparison, a synthesis study based on field experiments (van Groenigen et al., 2011) shows a narrower range of 0.3 %–0.6 % average increase for every 20 ppm increase, assuming a linear fertilization effect between CO₂ concentration and eCH₄.

Our results show that an ensemble of process-based wetland methane models provides quantification for uncertainty in eCH₄, as well as better constraints than a single model on the predicted trend and magnitude of eCH₄. However, nominally distinct models might have similar biases because of similarities in the way they represent a subset of processes (see Table S1 for the model summary). Future evaluation of modeled processes, such as oxidation, production, and transport pathways, along with model errors across different timescales using statistical tools, could help identify similarities in model behaviors to guide model development (Zhang, 2023b). Furthermore, the eCH₄ estimates are subject to forcing uncertainty, given that the two climate datasets applied in the simulation protocol do not cover the full magnitude and variability of climatic variables. Specifically, precipitation has a significant impact on wetland extent and anaerobic soil conditions but has large uncertainty in spatiotemporal patterns (Sun et al., 2018). Thus, we recommend that future ensemble simulations consider the uncertainty in climate variables among different datasets. In addition, the sensitivity parameters derived from the multiple regression are not independent of climate datasets. Thus, they are affected by the choice of meteorological drivers. Overall, quantitatively accounting for model performance and dependence and thoroughly evaluating the effectiveness (Chang et al., 2023) could improve the wetland model ensemble estimation in future studies.

Code and data availability

The code for the wetland models is available upon request from the respective model groups. The wetland ensemble results are publicly available in the Zenodo repository at 10.5281/zenodo.11309188 (Zhang, 2024). The wetland estimates from individual models are available upon request from the respective model groups. The FLUXNET-CH₄ dataset is publicly available at the following link: https://fluxnet.org/data/fluxnet-ch4-community-product/ (last access: 28 December 2024, FLUXNET-CH4 Community Product, 2024). The UpCH₄ dataset can be found at the link in McNicol et al. (2023).

The supplement related to this article is available online at: https://doi.org/10.5194/bg-22-305-2025-supplement.

Author contributions

BP and ZZ designed the simulation experiment with contributions from JM and WR. ZZ conducted data collection and data analysis. JM, WR, GB, PC, NG, PH, AI, AJ, FJ, TK, TL, XL, PM, JM, CP, SP, ZQ, QS, HT, XX, YY, XY, WZ, QinZ, QiuZ, QiaZ, and ZZ performed the simulations. ZZ prepared the manuscript with contributions from all co-authors.

Competing interests

At least one of the (co-)authors is a member of the editorial board of Biogeosciences. The peer-review process was guided by an independent editor, and the authors also have no other competing interests to declare.

Disclaimer

Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.

Acknowledgements

This paper is the result of a collaborative international effort under the umbrella of the Global Carbon Project, a project of Future Earth, and a research partner of the World Climate Research Programme. Joe R. Melton thanks Jade Skye for her assistance in running and processing the CLASSIC model simulations.

Financial support

Zhen Zhang acknowledges support from the National Natural Science Foundation of China Basic Science Center for Tibetan Plateau Earth System project. Yi Xi and Shushi Peng were funded by NSFC (41830643, 41722101). Thomas Kleinen acknowledges support from the German Federal Ministry of Education and Research ((BMBF) grant no. 01LP1921A). Akihiko Ito and Prabir Patra were partly supported by the Arctic Challenge for Sustainability II (ArCS-II) project (grant no. JPMXD1420318865), funded by Japan's Ministry of Education, Culture, Sports, Science and Technology (MEXT). Xiangyu Liu and Qianlai Zhuang are supported by a NASA project (NNX17AK20G). Qiuan Zhu and Changhui Peng are supported by the Second Tibetan Plateau Scientific Expedition (2019QZKK0304). Jurek Müller, Qing Sun, and Fortunat Joos were supported by the Swiss National Science Foundation (#200020_200511). Qing Zhu and William Riley were supported by the Reducing Uncertainties in Biogeochemical Interactions through Synthesis and Computation (RUBISCO) Scientific Focus Area and Energy Exascale Earth System Modeling Project, which are sponsored by the Earth and Environmental Systems Modeling (EESM) program under the Office of Biological and Environmental Research of the U.S. Department of Energy Office of Science. Yuanzhi Yao and Hanqin Tian are funded in part by an NSF program (award number: #1903722) and the NASA CMS Program (award number: NX14AO73G). Tingting Li was supported by the National Key Scientific and Technological Infrastructure project “Earth System Science Numerical Simulator Facility” (EarthLab) and the Open Research Program of the International Research Center of Big Data for Sustainable Development Goals (grant no. CBAS2023ORP02). Peter Hopcroft was supported by a Birmingham Fellowship and the University of Birmingham's BlueBEAR high-performance computing (HPC) service. Wenxin Zhang acknowledges support from the LUNARC computation project LU 2021/2-114 and the Swedish Research Council (Vetenskapsrådet), starting grant 2020-05338. Wenxin Zhang and Paul A. Miller acknowledge this study as a contribution to the strategic research areas Modeling the Regional and Global Earth System (MERGE) and Biodiversity and Ecosystem Services in a Changing Climate (BECC) at Lund University. Robert B. Jackson acknowledges support from the United Nations Environment Programme (UNEP) to Stanford University DTIE21-EN3143. Nicola Gedney was supported by the Newton Fund through the Met Office Climate Science for Service Partnership Brazil (CSSP Brazil). Atul K. Jain and Xiaoming Xu were supported by the US National Science Foundation (NSF-831361857) and would like to acknowledge the high-performance computing support from Cheyenne (https://www.cisl.ucar.edu/ncar-supercomputing-history/cheyenne, last access: 15 January 2025, Cheyenne HPE, 2016), provided by NCAR’s Computational and Information Systems Laboratory, sponsored by the National Science Foundation. Philippe Ciais acknowledges support from the space Agency Climate Change Initiative (ESA CCI) RECCAP2 project (grant no. ESRIN/4000123002/18/I-NB). Josep G. Canadell acknowledges the support of the Australian National Environmental Climate Science Program – Climate Systems hub. Gavin McNicol acknowledges support from the NASA CMS program (award number: NNH20ZDA001N).

Review statement

This paper was edited by Tyler Cyronak and reviewed by two anonymous referees.