# Improved global wetland carbon isotopic signatures support post-2006 microbial methane emission increase

Youmi Oh, Qianlai Zhuang, Lisa R. Welp, Licheng Liu, Xin Lan, Sourish Basu, Edward J. Dlugokencky, Lori Bruhwiler, John B. Miller, Sylvia E. Michel, et al.

► To cite this version:

Youmi Oh, Qianlai Zhuang, Lisa R. Welp, Licheng Liu, Xin Lan, et al.. Improved global wetland carbon isotopic signatures support post-2006 microbial methane emission increase. Communications Earth & Environment, 2022, 3, 10.1038/s43247-022-00488-5. insu-03824347

HAL Id: insu-03824347

https://insu.hal.science/insu-03824347

Submitted on 21 Oct 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Improved global wetland carbon isotopic signatures support post-2006 microbial methane emission increase

Youmi Oh<fnref n="1" />, ², ³✉, Qianlai Zhuang<fnref n="1" />, ⁴✉, Lisa R. Welp<fnref n="1" />, ⁵, Licheng Liu<fnref n="1" />, ¹², Xin Lan<fnref n="1" />, ², ³, Sourish Basu<fnref n="1" />, ⁶, ⁷, Edward J. Dlugokencky³, Lori Bruhwiler³, John B. Miller<fnref n="1" />, ³, Sylvia E. Michel⁸, Stefan Schwietzke<fnref n="1" />, ⁹, Pieter Tans³, Philippe Ciais¹⁰ & Jeffrey P. Chanton¹¹

Atmospheric concentrations of methane, a powerful greenhouse gas, have strongly increased since 2007. Measurements of stable carbon isotopes of methane can constrain emissions if the isotopic compositions are known; however, isotopic compositions of methane emissions from wetlands are poorly constrained despite their importance. Here, we use a process-based biogeochemistry model to calculate the stable carbon isotopic composition of global wetland methane emissions. We estimate a mean global signature of -61.3 ± 0.7‰ and find that tropical wetland emissions are enriched by ~11‰ relative to boreal wetlands. Our model shows improved resolution of regional, latitudinal and global variations in isotopic composition of wetland emissions. Atmospheric simulation scenarios with the improved wetland isotopic composition suggest that increases in atmospheric methane since 2007 are attributable to rising microbial emissions. Our findings substantially reduce uncertainty in the stable carbon isotopic composition of methane emissions from wetlands and improve understanding of the global methane budget.

<footnote n="1">¹Department of Earth, Atmospheric, and Planetary Sciences, Purdue University, West Lafayette, IN, USA. ²Cooperative Institute for Research in Environmental Sciences, University of Colorado, Boulder, CO, USA. ³Global Monitoring Laboratory, NOAA, Boulder, CO, USA. ⁴Department of Agronomy, Purdue University, West Lafayette, IN, USA. ⁵Purdue Climate Change Research Center, West Lafayette, IN, USA. ⁶Global Modeling and Assimilation Office, NASA Goddard Space Flight Center, Greenbelt, MD, USA. ⁷Earth System Science Interdisciplinary Center, University of Maryland, College Park, MD, USA. ⁸Institute of Arctic and Alpine Research, University of Colorado, Boulder, CO, USA. ⁹Environmental Defense Fund, Berlin, Germany. ¹⁰Laboratoire des Sciences du Climat et de l'Environnement, Gif-sur-Yvette, France. ¹¹Department of Earth, Ocean, and Atmospheric Sciences, Florida State University, Tallahassee, FL, USA. ¹²Present address: Department of Bioproducts and Biosystems Engineering, University of Minnesota, St. Paul, MN, USA.</footnote>

<footnote>✉email: youmi.oh@noaa.gov; qzhuang@purdue.edu</footnote>
Methane (CH₄) is a powerful greenhouse gas, and its atmospheric abundance (in nmol mol⁻¹, abbreviated ppb) has increased by about 170% since the 1750s¹,². Unlike the steady increases of atmospheric CO₂ and N₂O, atmospheric CH₄ nearly stabilized from 1998 to 2006 and then rapidly increased with a growth rate averaging ~6 ppb yr⁻¹ between 2007 and 2013 and ~11 ppb yr⁻¹ between 2014 and 2021. Since 2007, CH₄ has increased while its stable carbon isotopic composition (δ¹³C-CH₄, Eq. 1) has trended to more negative values, after increasing for 200 years³⁻⁵. Diagnosing the mechanisms behind these changes continues to generate considerable attention and controversy⁶⁻¹².

Measurements of atmospheric CH₄ abundance and δ¹³C-CH₄, in combination with isotopic signatures of sources and sinks, allow partitioning of CH₄ budgets into different source categories. This is because isotopic signatures of source categories differ substantially, where the δ¹³C-CH₄ of microbial sources (mean of -61.7 with variability of 6.2‰) is isotopically more depleted than fossil (mean of -44.8 with variability of 10.7‰) and biomass burning (mean of -26.2 with variability of 4.8‰) sources⁹,¹³. The destruction of CH₄, primarily by reaction with hydroxyl radical (OH), isotopically enriches atmospheric CH₄ relative to the emission-weighted source signature⁷,¹⁴,¹⁵. Due to a wide range of δ¹³C-CH₄ in each source category¹³, spatial and temporal distributions must be known to reduce the uncertainty in source partitioning. Wetlands are the largest single natural CH₄ source and strongly influence atmospheric δ¹³C-CH₄ changes¹², but the spatial and temporal information of wetland δ¹³C-CH₄ is limited, and often a single uniform value is assumed¹⁵,¹⁶. Studies show that source partitioning in atmospheric modeling is highly sensitive to spatiotemporal understanding of wetland δ¹³C-CH₄.⁹

Observations of global wetland δ¹³C-CH₄ show that CH₄ emitted from boreal wetlands is isotopically more depleted than CH₄ emitted from the tropics¹⁷⁻¹⁹; proposed causes include the abundance of C₄ plants influencing the δ¹³C of precursor organic matter (POM) (δ¹³C-POM), differences in CH₄-producing archaea (methanogen) communities, and different CH₄ transport processes¹⁸,²⁰⁻²². Ganesan et al.²³ produced a spatially-resolved global wetland δ¹³C-CH₄ distribution, but their study did not simulate temporal variability and did not fully represent fractionation processes that change based on meteorology, soil and vegetation properties.

Here, we incorporate a carbon isotope module into a bio-
geochemistry model, the Terrestrial Ecosystem Model (TEM)²⁴,²⁵
to simulate and mechanistically understand the global wetland
δ¹³C-CH₄ distribution. The model is evaluated using site-level
and regional observations. We then use this model to understand
the mechanisms behind the spatial and temporal variability of
wetland δ¹³C-CH₄, and conduct uncertainty and sensitivity tests.
Finally, we investigate the effect of new wetland isotope maps on
atmospheric δ¹³C-CH₄ and global CH₄ emissions by using an
atmospheric model and atmospheric observations⁵,²⁶.

Results

Modeling wetland δ¹³C-CH₄ dynamics. TEM simulates CH₄ production, oxidation, and transport between soils and the atmosphere (Eqs. 3–10)²⁴,²⁵,²⁷,²⁸. A carbon isotope-enabled module is incorporated into TEM, referred to as isoTEM, which explicitly considers carbon isotopic fractionation processes in wetlands (Fig. 1). The isotopic fractionation factor (α) for each process is defined in Eq. ²⁰, where α is larger than 1 when the product is isotopically more depleted than the reactant.

δ¹³C-POM is determined by the global C₃ and C₄ plant
distribution (Supplementary Fig. 1)²⁹, where C₄ vegetation is
isotopically enriched due to its photosynthetic pathway³⁰. We

<figure><img src="image_3.png" /><figcaption>**Fig. 1** Schematic diagram of wetland CH₄ dynamics and fractionations for isoTEM. The model simulates δ¹³C of precursor organic matter (POM) (δ¹³C-POM), CH₄ production, oxidation, and transport to the surface. δ¹³C-POM is determined by global C₃/C₄ plant distribution and long-term trends of atmospheric δ¹³C-CO₂. CH₄ is produced by two pathways, one using H₂ and CO₂ and another using acetate, with fractionation factors (α) for HMs (α<sub>HM</sub>) ≈ 1.030–1.080 and for AMs (α<sub>AM</sub>) ≈ 1.000–1.040. Produced CH₄ is partly oxidized by methanotrophs with a fractionation factor α<sub>MO</sub> ≈ 1.015–1.035. Residual produced CH₄ is emitted to the surface via three processes, plant-mediated transport (TP), diffusion (TD), and ebullection (TE), with different fractionations, α<sub>TP</sub> ≈ 1.000–1.030, α<sub>TD</sub> ≈ 1.005, α<sub>TE</sub> ≈ 1.000, respectively (Supplementary Tables 2–4 and Method “Model development, Model optimization”). Bold and dashed lines in the figure refer to chemical and transport processes, respectively.</figcaption></figure>

incorporated observed long-term trends of atmospheric
δ¹³C-CO₂ into soil δ¹³C-POM (Supplementary Fig. 2)³¹⁻³³.
CH₄ is produced from POM in anaerobic soils by two distinct
methanogen communities: hydrogenotrophic methanogens
(HMs) which use H₂ and CO₂ and acetoclastic methanogens
(AMs) which use acetate³⁴. The fractional contribution of these
pathways is important because HMs produce isotopically more
depleted CH₄ compared to AMs (α<sub>HM</sub> and α<sub>AM</sub> in Eq. 12)¹⁹,³⁵. To
quantify the fractional contribution, we used in situ observations
from Holmes et al.¹⁹ and conducted a regression analysis between
the fractional contribution and main environmental factors,
including soil pH, carbon, and latitude (Eq. 11, Supplementary
Fig. 3, and Supplementary Table 1). Total produced δ¹³C-CH₄ is
then calculated using a mixing of CH₄ pools from the two
methanogen communities (Eqs. 13–14). The CH₄ produced is
partly oxidized by methanotrophs in aerobic soil layers³⁶ with
¹²CH₄ being oxidized preferentially relative to ¹³CH₄ (α<sub>MO</sub> in
Eq. 15). Then, the remaining CH₄ is emitted to the atmosphere
through three processes: plant-mediated transport, diffusion, and
ebullection, with fractionation factors of α<sub>TP</sub>, α<sub>TD</sub>, and α<sub>TE</sub>,
respectively (Eq. 16)²⁰. We calculated oxidized and emitted
δ¹³C-CH₄ using the ratio of oxidation and transport processes
and their fractionation factors (Eqs. 17–22) (Method “Model
development”).

We optimized four fractionation factors related to CH₄ production, oxidation, and plant-mediated transport (α<sub>HM</sub>, α<sub>AM</sub>, α<sub>MO</sub>, α<sub>TP</sub>) using field observations in boreal (50–90°N), temperate (30–50°N/S), and tropical (<30°N/S) wetlands<sup>35,37,38</sup> (Eqs. 12, 15, 16, Supplementary Table 2–4 and Supplementary Figs. 4, 5). We
<figure><img src="image_4.png" /><figcaption>**Fig. 2** Global distribution of wetland δ¹³C-CH₄ and its latitudinal and long-term gradients simulated by isoTEM. **a** Modeled global wetland δ¹³C-CH₄ for wetland grid cells with static inundation data⁴⁹. **b** Mean latitudinal distribution of δ¹³C of POM (yellow), produced CH₄ (red), and CH₄ emitted to the atmosphere for all grid cells (blue) and flux-weighted grid cells (purple). **c** Long-term trends of global mean wetland δ¹³C-CH₄ with and without incorporating long-term trend in δ¹³C-POM (blue and purple, respectively). The shaded area in (b, c) represents one standard deviation determined from 20 ensembles of simulations where the optimized parameters were varied.</figcaption></figure>

set α<sub>TE</sub> to 1.000 and α<sub>TD</sub> to 1.005 based on previous studies²⁰ since ebullition and diffusion are governed by physical processes. To quantify uncertainties in model simulations, we used 20 ensemble members of optimization. We simulated global wetland CH₄ fluxes and their isotopic signatures during 1984–2016 at a spatial resolution of 0.5° with a 50-year spin-up to let δ¹³C-CH₄ of carbon pools come to a steady state (Methods “Model optimization, Simulation setup”).

**Simulated wetland δ¹³C-CH₄ and its comparison with observations.** We estimated the mean global wetland source signature to be -61.3 ± 0.7‰ during 1984–2016 (Fig. 2a). This value is more enriched than the mean wetland signature of -62.3‰ in Ganesan et al.²³ but similar to the mean value of -61.5‰ reported in Sherwood et al.¹³ (Supplementary Figs. 8, 9). The latitudinal distribution of δ¹³C-CH₄ ranges from a mean of -57 ± 3‰ in the tropics to -68 ± 4‰ in boreal regions (Fig. 2b). Our model simulates isotopically depleted global δ¹³C-CH₄ during the summer due to larger emissions from boreal regions (Supplementary Fig. 10) and a long-term trend of -0.7 ± 0.1‰ during 1984–2016 (blue line in Fig. 2c) when incorporating the long-term trend in δ¹³C-POM (Supplementary Fig. 2).

We compared the magnitude and spatial variability of the simulated wetland δ¹³C-CH₄ with site-level observations (Method “Model-data comparison”). We used 70 in situ measurements of global wetland δ¹³C-CH₄ from previous studies after excluding the measurements applied for optimization (Supplementary Data 1, Supplementary Fig. 11)¹³,¹⁹. We showed that isoTEM reduced the root mean square error (RMSE) by 40% compared to Ganesan et al.²³ (2.2 vs. 3.6) (Fig. 3a, b). Compared to a static isoTEM map in July, 2016, temporally-varying isoTEM reduced the RMSE slightly (2.2 vs. 2.4) (Supplementary Fig. 12). Ganesan et al.²³ prescribed maximum and minimum values as boundary conditions, resulting in unrealistic clusters of wetland δ¹³C-CH₄

near -65‰ for boreal and -60‰ for tropical sites (Fig. 3a and Supplementary Fig. 9).

Furthermore, we compared the spatial variability of simulated wetland δ¹³C-CH₄ with estimated signatures from airborne measurements for three regions in Alaska during 2012-2013 and 2015 using Miller-Tans plots (Fig. 3c-e) (Method “Model data comparison”)³⁹. In situ flux observations collected across Alaskan wetlands show an average of -65‰, but with a large 9‰ variance⁴⁰, which could be due to changes in wetland habitat including soil nutrients, pH, and vegetation distribution. The estimated signatures from observation also show that compared with δ¹³C-CH₄ from the North Slope of Alaska (-65 ± 1‰), δ¹³C-CH₄ from interior Alaska is more depleted (-69 ± 6) and δ¹³C-CH₄ from southwest Alaska is more enriched (-59 ± 4‰) (Supplementary Fig. 13 and Supplementary Table 5). IsoTEM reproduces the spatial variability (-67 ± 1, -68 ± 1, and -61 ± 2‰) for North Slope, interior, and southwest Alaska, respectively), whereas Ganesan et al.²³ simulated no spatial variability around -65‰ (Fig. 3e). IsoTEM simulates the spatial variability as the model optimized parameters for vegetated and non-vegetated sites separately and incorporated meteorology and soil inputs that vary spatially and temporally.

**Mechanistic understanding of spatial and temporal variability of wetland δ¹³C-CH₄.** We investigated the relative importance of the isotopic fractionation processes that affect the latitudinal gradient of wetland δ¹³C-CH₄ (Fig. 2b and Supplementary Fig. 14). First, compared to the boreal zone, δ¹³C-POM is enriched in the tropics by 5 ± 2‰ as C₄ plants are more prevalent (yellow line in Fig. 2b, Supplementary Figs. 1, 14a). Second, due to a larger fraction of AM in the tropics (Supplementary Fig. 3), the δ¹³C-CH₄ produced by methanogens is enriched by 12 ± 3‰ (red line in Fig. 2b, Supplementary Fig. 14b). Third, δ¹³C-CH₄ emitted from wetlands is 6 ± 4‰ more depleted in the tropics due to a larger proportion of plant-mediated transport causing higher
<figure><img src="image_5.png" /><figcaption>**Fig. 3** Site-level and regional model-data comparison of wetland $\delta^{13}C-CH_4$. **a, b** Site-level model-data comparison of observations with (a) Ganesan et al.23 and (b) temporally-varying isoTEM. **c-e** Regional model-data comparison of simulated wetland $\delta^{13}C-CH_4$ in Alaska by (c) Ganesan et al.23 and (d) isoTEM, and (e) their comparison with observation-based source signatures from NOAA aircraft measurements. The source signature is derived using Miller-Tans plots39. All observation data used for site-level comparison are listed in Supplementary Data 1. Error bars for observations in (a, b, e) represent one standard deviation of measured/inferred wetland $\delta^{13}C-CH_4$. Error bars for isoTEM in panel e represent one standard deviation determined from 20 ensemble simulations where the optimized parameters were varied.</figcaption></figure>

effective transport fractionation ($\alpha_T$) (blue line in Fig. 2b, Eq. 19, Supplementary Figs. 14d, 15, 16). Thus, in our simulation, $\delta^{13}C-CH_4$ emitted from tropical wetlands is enriched by ~11‰ compared to boreal wetlands. This difference is strengthened due to the distribution of $C_4$ plants (+5 ± 2‰) and the fractional contribution of differing methanogen communities (+12 ± 3‰) but weakened due to plant-mediated transport (-6 ± 4‰).

The long-term decrease in wetland $\delta^{13}C-CH_4$ simulated by isoTEM is mostly due to the decrease in atmospheric $\delta^{13}C-CO_2$32. The decreasing trend is incorporated into $\delta^{13}C$-POM (Supplementary Fig. 2) and causes the long-term decrease in wetland $\delta^{13}C-CH_4$ of ~0.7‰ from 1984 to 2016 (blue line in Fig. 2c)31. We conducted a simulation without the decreasing trend in $\delta^{13}C$-POM, which showed that increased temperature caused plant productivity and plant-mediated transport to increase and $\delta^{13}C-CH_4$ to decrease by ~0.1‰ during 1984–2016 (purple line in Fig. 2c and Supplementary Fig. 15). This implies that wetland $\delta^{13}C-CH_4$ could further change in the future due to decreases in $\delta^{13}C$-POM and increases in plant-mediated transport.

There is no continuous long-term measurements of wetland $\delta^{13}C-CH_4$ to verify our simulated long-term trend. Instead, we ran a regression analysis using observations collected from various wetland locations since the early 1980s (Supplementary Data 1) (Method "Uncertainty and sensitivity tests"). The results show that the representation of data increases when adding year as a parameter for the regression analysis (Supplementary Table 6), and the observed data show a long-term decreasing trend with year (~-0.1‰ year⁻¹) (Supplementary Fig. 17). More

continuous long-term observations of wetland $\delta^{13}C-CH_4$ are necessary to further verify the simulated long-term trends in wetland $\delta^{13}C-CH_4$.

**Uncertainty and sensitivity tests.** The version of TEM that we use for this study explicitly simulates soil $CO_2$ and $CH_4$ but not soil $H_2$ and acetate pools27, because the spatial and temporal soil $H_2$ and acetate pools are highly uncertain, and it is hard to verify the simulated pool changes with limited observations. On the contrary, the $CH_4$ production, oxidation, and transport processes in TEM have been thoroughly validated for global regions from previous studies24,25,27,41–44. Therefore, instead of introducing additional uncertainty from explicitly simulating $H_2$ and acetate pools that cannot be validated, we applied the observed fraction of different methanogen communities ($f_{HM}$) based on regression to the total $CH_4$ production rates simulated by TEM (Supplementary Fig. 3 and Supplementary Table 1). Thus, in our simulation, the fraction of HM and AM ($f_{HM}$) changes spatially but not temporally.

To quantify the uncertainty of our regression analysis of $f_{HM}$, we ran additional sensitivity tests by varying the $f_{HM}$ based on the uncertainty from Markov Chain Monte Carlo approach (Method "Uncertainty and sensitivity tests" and Supplementary Table 1)45. The results show that varying the parameters do not change the wetland $\delta^{13}C-CH_4$ substantially (<1%) (Supplementary Table 7). We acknowledge that this simplification would cause uncertainty in our model results, and future studies should explicitly measure
changes in H₂ and acetate concentrations in soils to incorporate the detailed processes into the model.

The simplification of CH₄ production processes may also cause uncertainty in the fractionation as we do not explicitly simulate fractionation processes from POM to CO₂/acetate and from CO₂/ acetate to CH₄. However, studies show that fractionation factors of the fermentation (POM to CO₂) and syntrophy (POM to acetate) processes are minor (α ≈ 1.00)¹⁹,⁴⁶,⁴⁷. There may be additional CO₂ produced by acetoclastic methanogenesis that have large fractionation (α ≈ 1.05), but the fraction is negligible from observations¹⁹. Thus, we believe our fractionation factors for HMs and AMs (αHM and αAM, respectively) reasonably represent the major fractionation processes of CH₄ production.

Furthermore, to quantify the influence of the uncertainty of our model inputs on simulation results, we varied temperature, precipitation, net primary productivity (NPP), atmospheric CH₄, and applied transient inundation maps⁴⁸ (Method “Uncertainty and sensitivity tests”). The results show that meteorology and substrate inputs alter mean wetland δ¹³C-CH₄ by ±1‰ (Supplementary Table 7). Our TEM simulations showed that CH₄ fluxes are sensitive to these inputs²⁷. However, δ¹³C-CH₄ shows small changes because the fractionation is determined by the fraction of CH₄ oxidation and transport processes (Eqs. 21, 22), that are calculated as a function of soil CH₄ production and the resultant CH₄ concentration changes (Cₘ in Eqs. 6–10). When CH₄ production increases due to input changes, CH₄ oxidation and transport increase simultaneously, causing minor variation in the fraction of oxidation and transport (Supplementary Fig. 16). Inundation datasets also alter wetland δ¹³C-CH₄ by changing the areas where wetland emissions occur (±2‰) (Supplementary Table 7 and Supplementary Figs. 6, 7).

**Implication for atmospheric modeling and global CH₄ budget.**

We constructed four scenarios with different wetland emissions and isotopic signature maps as inputs for TM5 atmospheric modeling during 1984–2016 to understand the impacts of spatially- and temporally-resolved wetland δ¹³C-CH₄ (Table 1). Scenario A uses a globally uniform value of wetland δ¹³C-CH₄; Scenario B uses a temporally static but spatially variable wetland isotope map from Ganesan et al.²³; and Scenario C uses spatially- and temporally-resolved maps from isoTEM. We used the same wetland fluxes²⁷ with a static inundation map⁴⁹ for Scenarios A–C that applied a step increase in fluxes in 2007 and 2014 by hypothesizing that microbial wetland emissions are the dominant driver of the post-2006 atmospheric CH₄ increase⁹,²⁶,⁵⁰ (46 Tgyr⁻¹ increase in total 2016 emissions across the global wetlands compared to the

averaged total emissions in 1999–2006) (Supplementary Fig. 19). However, since other studies have suggested an increase in fossil emission as a dominant driver for post-2006 CH₄ increases¹², we created scenario D that uses isoTEM wetland isotope maps with increases in both microbial and fossil emissions since 2007 (Table 1).

For Scenarios A–D, we adjusted global mean fossil and ruminant fluxes simultaneously to satisfy the long-term average mass balance of atmospheric CH₄ (Fig. 4a) and δ¹³C-CH₄ (Method “Forward modeling using TM5 atmospheric model”), as done by Lan et al.²⁶. These adjustments bring the long-term global average δ¹³C-CH₄ from simulation to the observed atmospheric levels without changing the post-2006 trends in simulated δ¹³C-CH₄⁹,²⁶. After adjustments, global mean fossil fluxes in scenarios A–D are between 170 and 190 Tgyr⁻¹ (Supplementary Fig. 19), within the uncertainty range in Schwietzke et al.⁹. For all other fluxes, their isotopic signatures, and CH₄ sinks that include OH, Cl, and O(¹D)¹⁴,⁵¹,⁵², we used the same setup in our model as in Lan et al.²⁶ (Supplementary Table 8). We compared simulated CH₄ and δ¹³C-CH₄ with observations from NOAA/INSTAAR global flask-air measurements (Supplementary Table 10)²,⁵.

The atmospheric simulation showed that Scenarios A–C follow the observed δ¹³C-CH₄ trend reasonably closely (Fig. 4b). However, Scenario D, which hypothesizes a post-2006 increase in microbial and fossil fluxes, does not follow the decreasing trend in global mean δ¹³C-CH₄. As pointed out earlier⁸,⁹,²⁶,⁵⁰, the magnitude of the δ¹³C-CH₄ decrease suggests that the increase in microbial emissions dominates fossil emissions in the post-2006 global CH₄ increase. We also confirmed a dominant increase in post-2006 microbial emissions, even though the long-term decrease in wetland δ¹³C-CH₄ of ~0.7‰ allow for a larger fossil emission increase. An additional simulation of Scenario C without including the long-term decrease in wetland δ¹³C-CH₄ shows differences of ~0.1‰ in simulated atmospheric δ¹³C-CH₄ in 2016 compared with model results with long-term wetland δ¹³C-CH₄ trend (Supplementary Fig. 23). This difference can accommodate more post-2006 emission increases from isotopically enriched fossil sources for Scenario C.

We differentiated Scenarios A–C by comparing their simulated latitudinal gradients of atmospheric δ¹³C-CH₄ with observations (Fig. 4c and Supplementary Fig. 20). The observed mean latitudinal gradient during 1998–2016 shows more negative δ¹³C-CH₄ at northern high latitudes compared to the Southern Hemisphere by 0.45 ± 0.05‰ (Supplementary Table 9), resulting from the dominance of northern emissions combined with the subsequent fractionation by reaction with OH during transport to the Southern Hemisphere¹⁷. Scenario C, which uses IsoTEM maps,

<table><caption>**Table 1 Setup of TM5 atmospheric modeling for Scenarios A-D.**</caption><thead><tr><th>Scenario</th><th>Wetland isotope map</th><th>Assumption of post-2006 CH<sub>4</sub> increase</th><th>Global mass balance of CH<sub>4</sub> and δ<sup>13</sup>C-CH<sub>4</sub><sup>a</sup></th></tr></thead><tbody><tr><td>A: Uniform w/Microbial Increase</td><td>One uniform value (–62.3‰, a mean signature of Ganesan et al.<sup>23</sup>)</td><td>Wetland emission increase (46 TgCH<sub>4</sub>yr<sup>–1</sup> increase from 1999–2006 to 2016)</td><td>Yes</td></tr><tr><td>B: Ganesan w/Microbial Increase</td><td>One spatial map from Ganesan et al.<sup>23</sup> (mean of –62.3‰)</td><td></td><td></td></tr><tr><td>C: isoTEM w/Microbial Increase</td><td>Spatio-temporally-resolved maps from isoTEM (mean of –61.3‰)</td><td></td><td></td></tr><tr><td>D: isoTEM w/ Microbial + Fossil Increase</td><td>(this study)</td><td>Wetland (60%) + fossil (40%) emission increase<sup>12</sup> (28 TgCH<sub>4</sub>yr<sup>–1</sup> increase from wetland, 18 TgCH<sub>4</sub>yr<sup>–1</sup> increase from fossil, from 1999–2006 to 2016)</td><td></td></tr></tbody></table>

<footnote>ªUsing a global mass balance model from previous studies⁹,²⁶, the long-term mean fossil and ruminant fluxes were adjusted from EDGAR 4.3.2 inventory to match the observed atmospheric growth rate of CH₄ during 1984–2016 and the 1998–2016 mean of δ¹³C-CH₄. By conducting the mass balance for all scenarios, we intended to reduce the spin-up time for atmospheric δ¹³C-CH₄ to be stabilized and compare all scenarios fairly (Method “Forward modeling using TM5 atmospheric model”).</footnote>
<figure><img src="image_6.png" /><figcaption>Fig. 4 Observed and simulated atmospheric CH₄ and δ¹³C-CH₄ from TM5 atmospheric modeling. a, b Model-data comparison of long-term trend of (a) atmospheric CH₄ from 1985 to 2016 (in ppb) and (b) δ¹³C-CH₄ from 1999 to 2016 (in ‰) by observation (gray) and simulations from Scenario A (yellow), B (red), C (blue), and D (skyblue). c Model-data comparison of normalized north-south gradient of atmospheric δ¹³C-CH₄ for Scenario A (yellow), B (red), and C (blue) in 2012. The north-south δ¹³C-CH₄ was calculated by zonally-averaging the surface δ¹³C-CH₄ and normalized based on the mean δ¹³C-CH₄ at 60-90 °S. The normalized north-south δ¹³C-CH₄ for other years is in Supplementary Fig. 20 and Supplementary Table 9. d Histogram of the difference between simulated and observed δ¹³C-CH₄ for Scenario A (yellow), B (red), and C (blue) for 6 measurement sites located in the northern hemisphere. The histogram plots for all measurement sites are in Supplementary Fig. 22. Information about Scenarios A-D is in Table 1.</figcaption></figure>

best reproduces the observed north-south gradient (0.48‰); Scenarios A and B under- and over-estimate the gradient by ~0.1‰ (0.37‰, and 0.59‰, respectively). The difference is also clear when comparing simulated atmospheric δ¹³C-CH₄ of Scenarios A–C at 10 measurement sites (Supplementary Figs. 21, 22 and Supplementary Table 10). The simulated and observed atmospheric δ¹³C-CH₄ differ the most at Northern Hemispheric sites, where Scenario C best reproduces the atmospheric δ¹³C-CH₄ data, but Scenario A and Scenario B simulate more negative and positive δ¹³C-CH₄, respectively (Fig. 4d).

The difference in north-south gradient of atmospheric δ¹³C-CH₄ between scenarios in Fig. 4c has an implication on regional partitioning of sources. Our sensitivity test of atmospheric modeling showed that all scenarios with transient inundation data⁴⁸ (Scenarios E–G) underestimated the north-south δ¹³C-CH₄ gradient (0.27 ± 0.06‰) compared with observations (0.45 ± 0.05‰) (Method “Forward modeling using TM5 atmospheric model”, Supplementary Table 11, Supplementary Figs. 26–30). Thus, we ran an additional Scenario H that increased emissions from boreal wetlands by 2.5 times over the original transient data (Supplementary Fig. 26 and Supplementary Table 11), which increased the north-south gradient by ~0.1‰ and improved the match with the observed north-south δ¹³C-CH₄ gradient (0.39‰) (Supplementary Figs. 29, 30).

Discussion

The atmospheric CH₄ burden has grown rapidly since 2007, and the largest annual increase since NOAA began measurements in 1983 was observed in 2021⁵³,⁵⁴. Since 2019, δ¹³C-CH₄ decreased more steeply⁵⁵, suggesting a further increase in microbial emissions

as this and other studies suggest⁸,⁹,²⁶,⁵⁰. The microbial sources include anthropogenic emissions from ruminants, agriculture, and waste, and natural emissions from wetlands and other aquatic ecosystems¹². Our simulation with increase in wetland emissions can reproduce the observed post-2006 δ¹³C-CH₄ decrease (Fig. 4), and our additional sensitivity test with increase in anthropogenic microbial emissions also tracks the post-2006 δ¹³C-CH₄ decrease (Supplementary Figs. 24, 25). However, the scenario with emission increase from both microbial and fossil sources did not reproduce the decreasing trend in atmospheric δ¹³C-CH₄ (Scenario D in Fig. 4). Other atmospheric studies that use atmospheric δ¹³C-CH₄ observations also showed that fossil emission increase is not a dominant reason of recent CH₄ increase²⁶,⁵⁶.

Atmospheric δ¹³C-CH₄ measurements have not been widely used to inform global methane budget because of uncertainty and spatiotemporal variation in source signatures, specifically citing limitation in wetland source signatures¹¹. In this study, we mechanistically explain the spatiotemporal variations of wetland δ¹³C-CH₄ and validate the simulation using regional, latitudinal, and global measurements, which substantially reduce the uncertainty in δ¹³C-CH₄ source signatures (Fig. 3). The small decreasing trend in wetland δ¹³C-CH₄ allow for more fossil emission increase in our estimate, but cannot change the conclusion that fossil emission increases are not the dominant driver for post-2006 global CH₄ increases.

This study considers wetland δ¹³C-CH₄ during the historical period only, but the future changes in wetland δ¹³C-CH₄ will depend on multiple factors. First, our simulation shows that changes in δ¹³C-POM affect wetland δ¹³C-CH₄ as SOC is mostly derived from new carbon from vegetation. The simulated active layer depth from a previous study⁵⁷ shows that the active layer
depth had a minor change during our simulation period (mean of old stored carbon in Arctic permafrost may play an important role as a substrate for methanogens in the future⁵⁸. Also, studies found the importance of microbial fossil CH₄ emissions from Arctic regions in the future⁵⁹,⁶⁰. The emissions are partially included as geologic seep emissions in our atmospheric modeling simulation (Supplementary Fig. 19 and Supplementary Table 8), and we also considered microbial fossil emissions with depleted δ¹³C-CH₄ in our total fossil emission estimates²⁶. Lastly, our simulation shows that the increase in NPP cause more plant-mediated transport. This effect will be more important in the future as plant functional types and plant growth change due to temperature increase.

There are several aspects of the model that could be improved. First, our optimization of fractionation factors was based on limited observations; additional long-term measurements of wetland δ¹³C-CH₄ would reduce the uncertainty. Second, the fractional contribution of two methanogen communities (HMs and AMs) changes spatially but not temporally in the model. We need a better understanding of temporal changes in methanogen communities especially following permafrost thaw and disturbance³⁵, and explicitly measure changes in H₂ and acetate concentrations in soils to incorporate detailed CH₄ production processes into the model. Third, various vertical methanogenic and non-methanogenic processes change δ¹³C of CH₄ and CO₂, the vertical CO₂/CH₄ ratios, and thus δ¹³C-CH₄ emitted from wetlands, since CO₂ is a substrate for HM⁶¹,⁶². We need to identify detailed vertical subsurface processes by conducting manipulation experiments using isotopic labeling analysis and inhibitor techniques to include those fractionation processes in future modeling studies⁶³. Fourth, current wetland models do not simulate large CH₄ emissions and δ¹³C-CH₄ from tropical tree stems and aquatic sources properly⁶⁴⁻⁶⁶. More measurements from these sources are crucial to improve the estimate of natural CH₄ emission and δ¹³C-CH₄ changes.

## Conclusion

To the best of our knowledge, this study is the first to use a biogeochemistry model to mechanistically explain and reduce the uncertainty in global wetland δ¹³C-CH₄. IsoTEM explains the latitudinal gradient of wetland δ¹³C-CH₄ that is increased by the distribution of C₃/C₄ plants and methanogen community type but decreased by plant-mediated transport. The long-term trends of the simulated wetland δ¹³C-CH₄ is controlled by δ¹³C-POM and plant-mediated transport. Our results suggest that rising microbial emissions is the dominant driver for the post-2006 global CH₄ increase and the concurrent decrease in atmospheric δ¹³C-CH₄, and the isoTEM spatial distribution of wetland δ¹³C-CH₄ better reproduces the observed atmospheric δ¹³C-CH₄ latitudinal gradient.

## Methods

**Model development.** We incorporated a carbon isotope module of methane (CH₄) into an existing process-based biogeochemistry model, the TEM (Fig. 1). The stable carbon isotope in delta notation (δ) describes the ratio of the heavy isotope to the light isotope in the sample (R<sub>sam</sub> = (¹³C/¹²C)<sub>sam</sub>) relative to a known standard ratio, R<sub>std</sub>, which is Vienna Pee Dee Belemnite (VPDB) for carbon²⁰ (Eq. 1). The deviation of this ratio-of-ratios from one is multiplied by 1000 to express isotope variations in parts per thousand (‰, permil). To express isotopic fractionation for the reaction A → B, we used a fractionation factor (α) defined in Eq. 2²⁰, where reactant A is in the numerator and product B is in the denominator. If α is larger than 1, the δ¹³C of product is isotopically more depleted in the heavy isotope than the δ¹³C of reactant, and if α is smaller than 1, the δ¹³C of product is more enriched in ¹³C than the δ¹³C of reactant.

$$ \delta^{13}C = (R_{sam}/R_{std}) - 1 \quad (1) $$

$$ \alpha = \frac{R_A}{R_B} = \left( \frac{\delta^{13}C_A}{1000} + 1 \right) / \left( \frac{\delta^{13}C_B}{1000} + 1 \right) \quad (2) $$

**Terrestrial Ecosystem Model (TEM).** TEM is a commonly used biogeochemistry model and its CH₄, soil, thermal, and hydrological dynamics have been evaluated in previous studies²⁴,²⁸,⁴¹⁻⁴⁴. The CH₄ dynamics module of TEM simulates CH₄ production, oxidation, and three transport processes—diffusion, ebullition, and plant-mediated transport—between soil and atmosphere. Please refer to the details of TEM in Oh et al.²⁵ and Liu et al.²⁷.

In TEM wetland model, changes in CH₄ concentrations (Cₘ) at depth z and time t (∂Cₘ(z,t)/∂t) are governed by Eq. 3, where Mₚ(z,t), M₀(z,t), Rₚ(z,t), and Rₐ(z,t) are CH₄ production, oxidation, plant-mediated transport, and ebullition rates, respectively, and ∂Fₐ(z,t)/∂z represents flux divergence from gaseous and aqueous diffusion. CH₄ is produced by methanogens in anaerobic soils (Mₚ) and is calculated by multiplying maximum potential production rate (Mₐₗ₀) and limiting functions of substrate, soil temperature, pH, and redox potentials (Sₘₐ, Mₘₗ, pH, and Rx, respectively) (Eq. 4). For this study, we assume that substrates for methanogens are mainly from soil organic carbon (SOC) derived from vegetation (Net Primary Productivity, NPP), where NPP(mon) is monthly NPP (gC m⁻² month⁻¹), NPPₘₐₓ is ecosystem-specific maximum monthly NPP, and f(Cₐₗₐₗ(z)) describes the relative availability of organic carbon substrate at depth z (Eq. 5). The substrate availability changes depending on atmospheric CO₂, meteorology, and soil properties⁶⁷.

$$ \frac{\partial C_M(z, t)}{\partial t} = M_P(z, t) - M_O(z, t) - \frac{\partial F_D(z, t)}{\partial z} - R_P(z, t) - R_E(z, t) \quad (3) $$

$$ M_{P,TEM}(z,t) = M_{GO}f(S_{OM}(z,t))f(M_{ST}(z,t))f(pH(z,t))f(R_x(z,t)) \quad (4) $$

$$ f(S_{OM}(z,t)) = \left( 1 + \frac{NPP(mon)}{NPP_{max}} \right) f(C_{DIS}(z)) \quad (5) $$

The produced CH₄ is partly oxidized by methanotrophs and is calculated by the multiplying the maximum potential oxidation rate (Oₘₐₓ) and limiting functions of CH₄ concentration, soil temperature, soil moisture, redox potential, nitrogen deposition, diffusion limited by high soil moisture, and oxygen concentration (Cₘₗ, Tₘₐₗₐₗ, Eₘₗ, Rₐₓ, Nₐₗ, Dₘₗ, and Cₐ₂ respectively) (Eq. 6). We use Michaelis-Menten kinetics with kₚₚₚₐₗₐₗ of 5 µM for the CH₄ limitation (Eq. 7).

$$ M_{O,TEM}(z,t) = O_{MAX}f(C_M(z,t))f(T_{SOIL}(z,t))f(E_{SM}(z,t))f(R_{OX}(z,t))f(N_{ap}(z,t))f(D_{MO}(z,t))f(C_{O₂}(z)) \quad (6) $$

$$ f(C_M(z,t)) = \frac{C_M(z,t)}{k_{CH4,LAM} + C_M(z,t)} \quad (7) $$

The remaining CH₄ is emitted to the surface with three different transport processes. First, gaseous and aqueous diffusion (Fₐ) occur due to concentration gradients of CH₄ (∂Cₘ(z,t)/∂t) (Eq. 8). The molecular diffusion coefficient (D) in different soil layers depends on soil texture and soil moisture. Ebullition (Rₐ) occurs when CH₄ bubble forms with Cₘ greater than µmol L⁻¹, and is calculated with a constant rate of Kₐ (1.0 h⁻¹) (Eq. 9). Plant-mediated transport (Rₚ) occurs for plants that function as a direct conduit for CH₄ to the atmosphere, and is functions of rate constant of 0.01 h⁻¹, vegetation type, root density, vegetation growth, and soil CH₄ concentrations (Kₚ, TRₚ, Vₚ, Rₚ, Fₚ, Rₐ, and Cₘ, respectively) (Eq. 10)⁶⁸. Rₚ depends on ecosystem-specific plant functional types and increases in a warmer soil due to the increase in vegetation growth. In TEM model, the soil profile was divided into 1-cm layers, and soil temperature, moisture, and CH₄ dynamics of TEM were simulated at an hourly time step²⁴,²⁷.

$$ F_D(z,t) = -D(z) \frac{\partial C_M(z,t)}{\partial t} \quad (8) $$

$$ R_E(z,t) = K_E f(C_M(z,t)) \quad (9) $$

$$ R_P(z,t) = K_P TR_{v,gg} f_{ROOT}(z) f_{GROW}(t) C_M(z,t) \quad (10) $$

**Methane stable carbon isotope module in TEM (isoTEM).** IsoTEM explicitly considers carbon isotopic fractionation processes for precursor organic matter (POM) and CH₄ during production, oxidation, and transport process. The δ¹³C of POM (δ¹³C-POM) is determined by the global C₃ and C₄ vegetation distribution²⁹ and is set to -27‰ and -13‰ for C₃- and C₄-only vegetation areas, respectively. The δ¹³C-POM for areas with mixed C₃ and C₄ vegetation is determined by the proportion of each type of photosynthetic pathway (Supplementary Fig. 1). We also incorporated long-term trends of atmospheric δ¹³C-CO₂ into soil δ¹³C-POM changes. Atmospheric δ¹³C-CO₂ became depleted in ¹³C by ≈2‰ during 1951-2016³,³³, and this signal is transferred to photosynthesates and POM for CH₄ emissions in wetlands⁶⁹. We incorporated this trend with a 6-year carbon residence time between photosynthesis and CH₄ emission in wetlands (Supplementary Fig. 2)³¹.

The CH₄ is then produced in anaerobic soils by two distinct methanogen communities: hydrogenotrophic methanogens (HMs) use H₂ and CO₂ and
acetoclastic methanogens (AMs) use acetate (CH₃COO⁻) for CH₄ production³⁴. Both mechanisms produce equimolar amounts of CO₂ and CH₄ from cellulose-like substrates. Using in situ observations from Holmes et al.¹⁹ the fractional contribution of the two methanogen communities is calculated based on a multiple regression analysis with the main environmental factors (Eq. 11). From the principal component analysis, Holmes et al.¹⁹ found a combination of environmental parameters including pH, vegetation type, soil organic carbon (SOC), and latitude are correlated with the dominant methanogenic pathway. The regression results show that fractional contribution of HMs (f<sub>HM</sub>) is positively correlated with latitude with a steep increase at 60°N (slope of 0.11 and 5.19 for latitudes below and above 60°N, respectively), and negatively correlated with pH (slope of -9.23) and SOC (slope of -0.7) (R² of 0.41, p < 0.001) (Eq. 11, Supplementary Table 1, and Supplementary Fig. 3).

$$f_{HM} = \begin{cases} a_1 \times \text{lat} + b \times \text{pH} + c \times \text{SOC} + d & \dots \text{for latitude} < \text{latitude}_{\text{step}} \\ a_1 \times \text{lat} + a_2 \times (\text{latitude} - \text{latitude}_{\text{step}}) + b \times \text{pH} + c \times \text{SOC} + d & \dots \text{for latitude} > \text{latitude}_{\text{step}} \end{cases} \qquad (11)$$

The δ¹³C-CH₄ produced by HMs and AMs more negative than the δ¹³C-POM, with the fractionation factors for HMs (α<sub>HM</sub>) ≈ 1.030–1.080 and for AMs (α<sub>AM</sub>) ≈ 1.000–1.040 (Eq. 12). The produced δ¹³C-CH₄ is calculated using a binary mixing of CH₄ pools from the two methanogen communities (Eqs. 13, 14).

$$\alpha_{HM} = \frac{1000 + \delta^{13}C_{POM}}{1000 + \delta^{13}CH_{4,prod,HM}}, \quad \alpha_{AM} = \frac{1000 + \delta^{13}C_{POM}}{1000 + \delta^{13}CH_{4,prod,AM}} \qquad (12)$$

$$\delta^{13}CH_{4,prod,HM} = \delta^{13}C_{POM} - 1000 \times \ln(\alpha_{HM}), \quad \delta^{13}CH_{4,prod,AM} = \delta^{13}C_{POM} - 1000 \times \ln(\alpha_{AM}) \qquad (13)$$

$$\delta^{13}CH_{4,prod} = f_{HM} \times \delta^{13}CH_{4,prod,HM} + (1 - f_{HM}) \times \delta^{13}CH_{4,prod,AM} \qquad (14)$$

The produced CH₄ is partly oxidized by methanotrophs in aerobic soils, which prefer ¹²CH₄, thus α for CH₄ oxidation (α<sub>MO</sub>) ≈ 1.015–1.035 (Eq. 15). Then, the produced CH₄ is transported to the atmosphere through three processes, plant-mediated transport, diffusion, and ebullition, with different fractionation factors α<sub>TP</sub> ≈ 1.000–1.030, α<sub>TD</sub> ≈ 1.000–1.010, α<sub>TE</sub> ≈ 1.000–1.005, respectively²⁰ (Eq. 16).

$$\alpha_{MO} = \frac{1000 + \delta^{13}CH_{4,prod}}{1000 + \delta^{13}CH_{4,oxid}} \qquad (15)$$

$$\alpha_{TP} = \frac{1000 + \delta^{13}CH_{4,prod}}{1000 + \delta^{13}CH_{4,TP}}, \quad \alpha_{TE} = \frac{1000 + \delta^{13}CH_{4,prod}}{1000 + \delta^{13}CH_{4,TE}}, \quad \alpha_{TD} = \frac{1000 + \delta^{13}CH_{4,prod}}{1000 + \delta^{13}CH_{4,TD}} \qquad (16)$$

We calculated the oxidized and transported δ¹³C-CH₄ based on “open system equations” at steady state to consider residual enriched CH₄ after oxidation and transport processes⁷⁰⁻⁷³. We approximated that CH₄ produced in the entire vertical soil column is either oxidized or transported in each hourly time step (Eq. 17). In Eqs. 17, 18, M<sub>P</sub>(z,t), M<sub>o</sub>(z,t), R<sub>P</sub>(z,t), and R<sub>E</sub>(z,t) represent CH₄ production, oxidation, plant-mediated transport, and ebullition rates, respectively, and ∂F<sub>D</sub>(z,t)/∂z represents flux divergence due to gaseous and aqueous diffusion for each soil layer z and time t. For simplicity, we defined effective transport fractionation, α<sub>T</sub>, by flux-weighting the proportions of fractionation factors of three transport processes in Eq. 19. The isotopic difference between oxidation and transport processes can be described by a fractionation factor, α<sub>T/MO</sub>, in Eq. 20. Given these conditions, isotopic signatures for oxidation and transport to the atmosphere (emission) can be written in Eqs. 21, 22. For more details, refer to Hayes⁷⁴.

$$\sum_{z} M_{P}(z, t) = \sum_{z} M_{o}(z, t) + \sum_{z} \frac{\partial F_{D}(z, t)}{\partial z} + \sum_{z} R_{P}(z, t) + \sum_{z} R_{E}(z, t) \qquad (17)$$

$$f_{ox} = \frac{\sum_{z} M_{o}(z, t)}{\sum_{z} M_{P}(z, t)}, f_{TP} = \frac{\sum_{z} R_{P}(z, t)}{\sum_{z} M_{P}(z, t)}, f_{TE} = \frac{\sum_{z} R_{E}(z, t)}{\sum_{z} M_{P}(z, t)}, f_{TD} = \frac{\sum_{z} \frac{\partial F_{D}(z, t)}{\partial z}}{\sum_{z} M_{P}(z, t)} \qquad (18)$$

$$\alpha_{T} = \frac{(f_{TP}\alpha_{TP} + f_{TE}\alpha_{TE} + f_{TD}\alpha_{TD})}{f_{TP} + f_{TE} + f_{TD}} \qquad (19)$$

$$\alpha_{T/MO} = \frac{\alpha_{MO}}{\alpha_{T}} = \epsilon_{T/MO} + 1 \qquad (20)$$

$$\delta^{13}CH_{4,oxid} = \frac{\delta^{13}CH_{4,prod} - (1 - f_{ox})\epsilon_{T/MO}}{\alpha_{T/MO}(1 - f_{ox}) + f_{ox}} \qquad (21)$$

$$\delta^{13}CH_{4,emitted} = \frac{\alpha_{T/MO}\delta^{13}CH_{4,prod} + f_{ox}\epsilon_{T/MO}}{\alpha_{T/MO}(1 - f_{ox}) + f_{ox}} \qquad (22)$$

**Model optimization.** We optimized 4 fractionation factors, α<sub>HM</sub>, α<sub>AM</sub>, α<sub>MO</sub>, and α<sub>TP</sub>, using in situ observations for six wetland ecosystem types (Eqs. 12, 15, 16). Since the fractionation factors for ebullition and diffusion are governed by physical processes, we set them as constants based on literature (α<sub>TE</sub> = 1.000, α<sub>TD</sub> = 1.005)²⁰. The wetland ecosystems are divided into forested and non-forested wetlands for boreal (50–90°N), temperate (30–50°N/S), and tropical (<30°N/S) regions. To optimize parameters, we collected observation data from six sites representing each ecosystem (Supplementary Tables 2–4)³⁵,³⁷,³⁸. For tropical wetlands, we used observation data from Burke et al.³⁸,⁷⁵. For forested wetlands, we used data from “Willow Marsh Trail” station, a swamp wetland dominated by hardwoods and *Lemnaceae*. For non-forested wetlands, we used data from “St. Petersburg” site where Sawgrass is the dominant vegetation. For temperate wetlands, we used data from Kelly et al.³⁷. For forested wetlands, we used data from “S2 Bog” where is entirely forested with *Picea mariana*. For non-forested wetlands, we used data from “Junction Fen” where is treeless and dominated by *Carex oligosperma*. For Arctic wetlands, we used data from McCalley et al.³⁵. For forested wetlands, we could not find δ¹³C-CH₄ data from the well-drained “Palsa” occupied by woody plants, mosses, and ericaceous. Thus, we used δ¹³C-CH₄ data from “Sphagnum” site that is in the transition between the Palsa and Eriophorum sites, and showed similar CH₄ fluxes as the “Palsa” site. For non-forested wetlands, we used data from the “Eriophorum” site.

Besides the observed meteorology from field sites, we also used CRU time-series version 4.01 to fill missing meteorological inputs⁷⁶. We then used the Shuffled Complex Evolution Approach in R language (SCE-UA-R) to minimize the difference between simulated and observed δ¹³C-CH₄⁷⁷. For each site, 20 ensembles were run using SCE-UA-R with 10,000 maximum loops per parameter ensemble, and all of them reached steady state before the end of the loops. Our optimization results show that isoTEM captures the magnitude and seasonality of observed soil CH₄ fluxes and δ¹³C-CH₄ (Supplementary Fig. 4).

**Simulation setup.** To estimate spatially- and temporally-varying δ¹³C-CH₄ from global wetlands, we used spatially explicit data of land cover, soil pH and textures, meteorology and leaf area index (LAI)²⁴,²⁷. Land cover, soil pH and textures were used to assign vegetation-specific and texture-specific parameters to a grid cell⁷⁸⁻⁸⁰. Meteorological inputs were derived from historical air temperature, precipitation, vapor pressure, and cloudiness from gridded CRU time-series version 4.01⁷⁶. We used monthly LAI derived from satellite imagery⁸¹ to prescribe LAI for each 0.5° × 0.5° grid cell. All other parameters except fractionation factors were set the same as in Liu et al.²⁷. We simulated global wetland CH₄ fluxes and their isotopic ratios between 1984 and 2016 at a spatial resolution of 0.5° × 0.5° with a 50-year spin-up to let the carbon isotopic composition of carbon pools come to a steady state.

Because various wetland inundation data exist⁸², we first assumed that every global land grid cell can potentially be saturated, thus this product can be used with any wetland inundation data in future studies. To fill the grid cells without wetland types, we set forested and non-forested wetlands based on global vegetation types (Supplementary Fig. 5). In our analyses, simulated ecosystem-specific δ¹³C-CH₄ from wetlands was flux weighted for each grid cell, based on CH₄ emissions simulated by TEM defined over the static inundation data from Matthews and Fung (Supplementary Fig. 6a)⁴⁹.

## Model-data comparison

**Site level.** We compared our model results with previously published data from 58 in situ measurements compiled by Holmes et al.¹⁹ and 66 in situ measurements by Sherwood et al.¹³. Holmes et al.¹⁹ compiled latitude, fraction of HM and AM, pH, vegetation, and δ¹³C-CH₄ to understand factors affecting the methanogenic pathway in global wetlands. The wetland database of Sherwood et al.¹³ includes literature reference, latitude, wetland types, and measurement methods. After combining overlapped data of Holmes et al.¹⁹ and Sherwood et al.¹³ and excluding data that we used for our model optimization³⁵,³⁷,³⁸, 70 sites remained for site-level validation (Supplementary Fig. 11 and Supplementary Data 1). Due to a possible mismatch of soil and vegetation properties, and wetland distribution of grid cells between model and observation, we compared observed δ¹³C-CH₄ with simulated δ¹³C-CH₄ of the sampling year within two adjacent grid cells (1° × 1°) of the observation.

**Regional level.** We used aircraft air samples from 3 regions in Alaska from the Carbon in Arctic Reservoirs Vulnerability Experiment (CARVE)⁸³,⁸⁴. From 2012 to 2015, CARVE collected airborne measurements of atmospheric chemical components and relevant land surface parameters in the Alaskan Arctic to provide insights into Arctic carbon cycling. During the flights, flask-air samples were collected then sent to NOAA GML for measurements of 50 trace gases including CO₂, CH₄, CO, OCS, NMHCs, and then sent to INSTAAR for and the isotopic composition of CO₂ and CH₄. After excluding airborne data with flags, there are 1476 measurements during the sampling period.
In situ flux observations collected across Alaskan wetlands show an average of -65‰ but a large 9‰ variation, due to the complex vegetation and soil properties⁴⁰. To compare the spatial variability of wetland δ¹³C-CH₄, we divided the Alaskan continent into three regions: North Slope, interior, and southwest Alaska based on latitude (62–68 °N, 52–62 °N and 140–155 °W, and 52–62 °N and 155–170 °W for North Slope, interior, and southwest Alaska, respectively). We used Miller-Tans plots to identify the source signatures of δ¹³C-CH₄ from wetlands using the airborne measurements³⁹. To identify wetland isotopic signatures, we removed measurements that may have effects from fossil fuel emission (C₃H₈ < 300 ppt), biomass burning (CO < 300 ppb), and transport influence (Altitude < 1500 m), and we set the background altitude to >5000 m. After plotting the data, 2014 was excluded due to limited data and small R² (Supplementary Table 5).

### Uncertainty and sensitivity tests

Long-term trends in wetland δ¹³C-CH₄ from observations. We considered latitude, pH, and soil carbon as key parameters that determine variability of wetland δ¹³C-CH₄ to run a linear regression using the site-level observations collected from global wetlands since the early 1980s (Supplementary Data 1). We added year as additional parameter for the linear regression and see if it improves the fit with data. The regression results show that wetland δ¹³C-CH₄ is negatively correlated with year, latitude, and SOC (slope of -0.11, -0.10, and -0.20, respectively), and positively correlated with pH (slope of 2.21) (R² of 0.30, p < 0.001) (Eq. 23, Supplementary Fig. 17, and Supplementary Table 6). The regression without year as a parameter showed smaller coefficient (R² of 0.25, p < 0.001).

$$ \delta^{13}C - \text{CH}_4 = a \times \text{lat} + b \times \text{pH} + c \times \text{SOC} + d \times \text{year} + e \quad (23) $$

Markov Chain Monte Carlo for the fraction of HM ($f_{HM}$). We used a Markov Chain Monte Carlo (MCMC) approach for parameter uncertainty estimation for $f_{HM}$. MCMC is a method for estimating the posterior probability density function for asset of parameters, given priors those parameters and a set of observations⁴⁵. We used independent, uniform prior probability density functions for each parameter in Supplementary Table 1. Thirty-nine data points from Holmes et al.¹⁹ were used to constrain the model. Gaussian errors were assumed. We generated a Markov chain with 100,000 elements to estimate the joint posterior probability density functions. The chain converged after about 10,000 elements. We used the posterior probability density function to estimate the uncertainty of parameter (Supplementary Table 1).

Sensitivity test with meteorological and substrate inputs, $f_{HM}$, and inundation. We conducted 8 sensitivity tests of meteorology and substrate inputs. Specifically, we altered air temperature by ±3 °C, precipitation by ±30%, and atmospheric CH₄ abundance by ±30%, and NPP by ±30%, uniformly for each grid cell, while maintaining all other variables at their default isoTEM values. We also varied parameters for $f_{HM}$ based on the uncertainty range from MCMC (Supplementary Table 1). We further varied a wetland inundation using satellite-driven Surface Water Microwave Product Series-Global Lakes and Wetlands Database (SWAMPS-GLWD)⁴⁸.

### Forward modeling using TM5 atmospheric model

Global mass balance for bottom-up inventory. We adjusted global long-term mean fossil fluxes to match the simulated growth rate of CH₄ during 1984–2016 and the 1998–2016 mean of δ¹³C-CH₄ with observation (Table 1 and Supplementary Table 11). Lan et al.²⁶ showed that there is an offset of simulated global mean δ¹³C-CH₄ when using EDGAR 4.3.2 inventory as the inventory underestimates fossil fluxes. To remove the offset and compare our scenarios fairly, we adjusted fossil fluxes between 170 and 190 TgCH₄yr⁻¹ (Supplementary Fig. 19), within the uncertainty range in Schwietzke et al.⁹. To satisfy the global mass balance, we ran one box model that included CH₄ sources of biogenic, fossil and biomass/biofuel emissions, with corresponding isotopic signatures, and CH₄ sinks due to reaction with OH, Cl, and O(¹D) and soil bacteria, all with different fractionation factor. When we increased or decreased fossil fluxes, we accordingly decreased or increased ruminant flux, respectively, so the total annual CH₄ fluxes followed the observed atmospheric CH₄ growth rate, and the long-term mean total emission was set to 536–538 TgCH₄yr⁻¹ during 1984–2016. For more details on the setup and equations for global mass balance, refer to Lan et al.²⁶.

Data sources for CH₄ emissions and its isotopic source signatures. We used the bottom-up inventory constructed by Lan et al.²⁶ (Supplementary Table 8). In specific, for CH₄ emissions, we used GFED 4.1 s for biomass burning for 1997–2016⁸⁵ and annual emissions from the Reanalysis of Trophospheric chemical composition project before 1997, and the EDGAR 4.3.2 inventory for other anthropogenic emissions for 1984–2016⁸⁶. For emissions from geological seeps, we used gridded emission from Etiope et al.⁸⁷. Emission estimates from wild animals and termites were adopted from Bergamaschi et al.⁸⁸. For δ¹³C-CH₄ source signature, fossil fuel source signature data were based on the global δ¹³C-CH₄ source signature inventory 2020⁹⁰, where the data were categorized by coal gas, conventional gas, and shale gas. Biomass burning, biofuel burning, ruminant, and wild animal δ¹³C-CH₄ data were based on the global maps of C₃/C₄ distribution²⁹. The geological seeps δ¹³C-CH₄ data were from Etiope et al.⁸⁷.

TM5 atmospheric modeling of CH₄ and δ¹³C-CH₄. Atmospheric CH₄ mole fractions and δ¹³C-CH₄ were simulated from 1984 to 2016 by coupling the surface fluxes and isotope source signatures from the bottom-up inventory with the TM5 tracer transport model driven by ECMWF ERA Interim meteorology with the 4DVAR branch of the TM5 model⁹⁰,⁹¹. TM5 was run globally at 6° × 4° over 25 vertical sigma-pressure hybrid levels, for total CH₄ and ¹³C-CH₄. For each source type, ¹³C-CH₄ fluxes were derived from total CH₄ fluxes and source-specific isotope source signatures. We spun up our model during 1984–1999 and selected 2000–2016 to compare with atmospheric observations to ensure our spin-up period was sufficient for equilibration of atmospheric δ¹³C-CH₄ inter-hemispheric gradient²⁶,⁹². As per Lan et al.²⁶ we applied tropospheric Cl sink of Hossaini et al.⁵¹ and the OH field from Spivakovsky et al.¹⁴ with a fractionation factor of -3.9‰. The CH₄ sinks varied spatially and seasonally but did not change interannually. For more details on setup for TM5 modeling, refer to Lan et al.²⁶.

Atmospheric CH₄ and δ¹³C-CH₄ measurement. Observational data of atmospheric CH₄ and δ¹³C-CH₄ used to evaluate model results are from flask-air measurements from NOAA's Global Greenhouse Gas Reference Network²⁶,⁵⁴. The flask-air samples was analyzed for δ¹³C-CH₄ at the Institute of Arctic and Alpine Research (INSTAAR), University of Colorado, Boulder. Gas chromatography-Isotope-ratio mass spectrometry (GC-IRMS) is used for δ¹³C-CH₄ analysis⁵⁵. The δ¹³C-CH₄ in air measurements are referenced against the Vienna Pee Dee Belemnite (VPDB) standard (Eq. 1). A subset of the observation sites predominantly influenced by well-mixed background air is used to construct a Marine Boundary Layer (MBL) zonally averaged surface using methods developed by Masarie and Tans (1995)⁹³, to represent the observational-based global long-term trend and north-south gradient. This includes 31 sites with CH₄ measurements during study period of 1984–2016 and 10 of which with δ¹³C-CH₄ measurements starting in 1998 (Supplementary Fig. 21 and Supplementary Table 10). More details on the MBL data products and uncertainties can be found at https://www.esrl.noaa.gov/gmd/cgg/mbl/mbl.html. For model-observation comparisons, model results from the same set of MBL sites are sampled, and the same calculation methods are applied to model results and observations for global long-term and north-south gradient. The north-south gradient was calculated as the difference of atmospheric δ¹³C-CH₄ between 60–90 °S and 60–90 °N.

Atmospheric modeling with transient inundation data for Scenarios E-H. Since we used static wetland inundation data⁴⁹ for our default Scenarios A–D, we used transient wetland inundation data from Poulter et al.⁴⁸ and ran TM5 atmospheric model (Supplementary Figs. 26–30 and Supplementary Table 11). Same as Scenarios A–C, we constructed Scenarios E-G with different wetland isotopic signature maps as inputs for TM5 atmospheric modeling in 1984–2016. In specific, the first uses a globally uniform wetland δ¹³C-CH₄ of -62.3‰, the mean wetland signature from Ganesan et al.²³ (referred to as Scenario E), the other uses a static wetland isotope spatial map from Ganesan et al.²³ (referred to as Scenario F), and the last used spatially- and temporally-resolved maps from isoTEM (referred to as Scenario G).

The wetland fluxes for Scenarios E-G are based on Liu et al.²⁷ and transient inundation⁴⁸ but applied an increase in fluxes after 2006 by hypothesizing that the microbial wetland emission is a dominant driver of post-2006 atmospheric CH₄ increase (Supplementary Fig. 26), same as Scenarios A–C. We also conducted the global mass balance by adjusting global long-term mean fossil fluxes between 160 and 180 TgCH₄yr⁻¹ for Scenarios E-G to match the simulated growth rate of CH₄ during 1984–2016 and the 1998–2016 mean of annual δ¹³C-CH₄ with observations.

Scenarios E-G reproduced the observed global CH₄ growth rate during 1984–2016 and the global long-term mean δ¹³C-CH₄ with observation during 1998–2016 (Supplementary Fig. 28), as we set the fluxes based on the mass balance. However, Scenarios E-G with transient inundation data underestimated the north-south δ¹³C-CH₄ gradient (0.27 ± 0.06‰) compared with observations (0.45 ± 0.05‰) (Supplementary Fig. 29). Thus, we ran an additional scenario H that increased emissions from boreal wetlands by 2.5 times over the original transient data (Supplementary Fig. 26 and Supplementary Table 11), which improved the match with the observed north-south δ¹³C-CH₄ gradient (0.39‰) (Supplementary Fig. 29). The site-level comparison with atmospheric δ¹³C-CH₄ from 10 observation sites also confirmed that Scenario H more closely reproduced the observation (Supplementary Fig. 30). This implies that the transient inundation data from Poulter et al.⁴⁸ may need more wetland emissions from boreal regions as found in static inundation data⁴⁹ (Supplementary Fig. 6) and other satellite-derived inundation data⁹⁴.

### Data availability

Supplementary Data 1 is available at: https://figshare.com/articles/dataset/Supplementary_Data_1_of_Oh et al_2022./19929965. The stable carbon isotopic composition of wetland emissions is available at: https://doi.org/10.25925/9s6n-g811.

### Code availability

The code is also archived and available at: https://doi.org/10.15138/cem6-ka15.
Received: 9 November 2021; Accepted: 22 June 2022

Published online: 12 July 2022

References

1. Etheridge, D. M., Steele, L., Francey, R. J. & Langenfelds, R. L. Atmospheric methane between 1000 AD and present: Evidence of anthropogenic emissions and climatic variability. *J. Geophys. Res. Atmos.* **103**, 15979–15993 (1998).

2. Dlugokencky, E.J., Crotwell, A.M., Mund, J.W., Crotwell, M.J. & Thoning, K.W. NOAA global greenhouse gas reference network flask-air sample measurements of CO₂, CH₄, CO, N₂O, H₂, SF₆ and isotopic ratios at global and regional background sites, 1967—Present. https://doi.org/10.15138/VNCZ-M766 (2022).

3. Dlugokencky, E. J., Nisbet, E. G., Fisher, R. & Lowry, D. Global atmospheric methane: budget, changes and dangers. *Philos. Trans. R. Soc. A Math. Phys. Eng. Sci.* **369**, 2058–2072 (2011).

4. Ferretti, D. F. et al. Atmospheric science: unexpected changes to the global methane budget over the past 2000 years. *Science* **309**, 1714–1717 (2005).

5. White, J.W.C., Vaughn, B. H. & Michel, S. E. University of Colorado, Institute of Arctic and Alpine Research (INSTAAR), Stable Isotopic Composition of Atmospheric Methane (¹³C) from the NOAA ESRL Carbon Cycle Cooperative Global Air Sampling Network, 1998–2017, Version: 2018-09-24. ftp://aftp.noaa.gov/data/trace_gases/ch4c13/flask/ (2018).

6. Hausmann, P., Sussmann, R. & Smale, D. Contribution of oil and natural gas production to renewed increase in atmospheric methane (2007–2014): top–down estimate from ethane and methane column observations. *Atmos. Chem. Phys.* **16**, 3227–3244 (2016).

7. Naus, S. et al. Constraints and biases in a tropospheric two-box model of OH. *Atmos. Chem. Phys.* **19**, 407–424 (2019).

8. Schaefer, H. et al. A 21st-century shift from fossil-fuel to biogenic methane emissions indicated by 13CH₄. *Science* **352**, 80–84 (2016).

9. Schwietzke, S. et al. Upward revision of global fossil fuel methane emissions based on isotope database. *Nature* **538**, 88–91 (2016).

10. Worden, J. R. et al. Reduced biomass burning emissions reconcile conflicting estimates of the post-2006 atmospheric methane budget. *Nat. Commun.* **8**, 1–11 (2017).

11. Turner, A. J., Frankenberg, C. & Kort, E. A. Interpreting contemporary trends in atmospheric methane. *Proc. Natl. Acad. Sci.* **116**, 2805–2813 (2019).

12. Saunois, M. et al. The global methane budget 2000–2017. *Earth system science data*, **12**, 1561–1623 (2020).

13. Sherwood, O. A., Schwietzke, S., Arling, V. A. & Etiope, G. Global inventory of gas geochemistry data from fossil fuel, microbial and burning sources, version 2017. *Earth Syst. Sci. Data* **9**, 639–656 (2017).

14. Spivakovsky, C. M. et al. Three-dimensional climatological distribution of tropospheric OH: Update and evaluation. *J. Geophys. Res. Atmos.* **105**, 8931–8980 (2000).

15. Strode, S. A. et al. Strong sensitivity of the isotopic composition of methane to the plausible range of tropospheric chlorine. *Atmos. Chem. Phys.* **20**, 8405–8419 (2020).

16. Rice, A. L. et al. Atmospheric methane isotopic record favors fossil sources flat in 1980s and 1990s with recent increase. *Proc. Natl. Acad. Sci.* **113**, 10791–10796 (2016).

17. Feinberg, A. I., Coulon, A., Stenke, A., Schwietzke, S. & Peter, T. Isotopic source signatures: Impact of regional variability on the δ¹³CH₄ trend and spatial distribution. *Atmos. Environ.* **174**, 99–111 (2018).

18. Brownlow, R. et al. Isotopic ratios of tropical methane emissions by atmospheric measurement. *Global Biogeochem. Cycles* **31**, 1408–1419 (2017).

19. Holmes, M. E., Chanton, J. P., Tfaily, M. M. & Ogram, A. CO₂ and CH₄ isotope compositions and production pathways in a tropical peatland. *Global Biogeochem. Cycles* **29**, 1–18 (2015).

20. Chanton, J. P. The effect of gas transport on the isotope signature of methane in wetlands. *Org. Geochem.* **36**, 753–768 (2005).

21. Fisher, R. E. et al. Measurement of the ¹³C isotopic signature of methane emissions from northern European wetlands. *Global Biogeochem. Cycles* **31**, 605–623 (2017).

22. Nakagawa, F., Yoshida, N., Nojiri, Y. & Makarov, V. Production of methane from alases in eastern Siberia: Implications from its ¹⁴C and stable isotopic compositions. *Global Biogeochem. Cycles* **16**, 14–1–14–15 (2002).

23. Ganesan, A. L. et al. Spatially resolved isotopic source signatures of wetland methane emissions. *Geophys. Res. Lett.* **45**, 3737–3745 (2018).

24. Zhuang, Q. et al. Methane fluxes between terrestrial ecosystems and the atmosphere at northern high latitudes during the past century: a retrospective analysis with a process-based biogeochemistry model. *Global Biogeochem. Cycles* **18**, GB3010 (2004).

25. Oh, Y. et al. Reduced net methane emissions due to microbial methane oxidation in a warmer Arctic. *Nat. Clim. Chang.* **10**, 317–321 (2020).

26. Lan, X. et al. Improved constraints on global methane emissions and sinks using δ¹³C-CH₄. *Global Biogeochem. Cycles* https://doi.org/10.1029/2021gb007000 (2021).

27. Liu, L. et al. Uncertainty quantification of global net methane emissions from terrestrial ecosystems using a mechanistically-based biogeochemistry model. *J. Geophys. Res. Biogeosci.* **125**, e2019JG005428 (2020).

28. Zhuang, Q. et al. Response of global soil consumption of atmospheric methane to changes in atmospheric climate and nitrogen deposition. *Global Biogeochem. Cycles* **27**, 650–663 (2013).

29. Still, C. J., Berry, J. A., Collatz, G. J. & DeFries, R. S. Global distribution of C₃ and C₄ vegetation: Carbon cycle implications. *Global Biogeochem. Cycles* **17**, 6–1–61–4 (2003).

30. Holmes, M. E., Chanton, J. P., Bae, H. S. & Ogram, A. Effect of nutrient enrichment on δ¹³CH₄ and the methane production pathway in the Florida Everglades. *J. Geophys. Res. Biogeosci.* **119**, 1267–1280 (2014).

31. Lassey, K. R., Etheridge, D. M., Lowe, D. C., Smith, A. M. & Ferretti, D. F. Centennial evolution of the atmospheric methane budget: What do the carbon isotopes tell us? *Atmos. Chem. Phys.* **7**, 2119–2139 (2007).

32. Keeling, R. F. et al. Atmospheric evidence for a global secular increase in carbon isotopic discrimination of land photosynthesis. *Proc. Natl. Acad. Sci. USA* **114**, 10361–10366 (2017).

33. Graven, H. et al. Compiled records of carbon isotopes in atmospheric CO₂ for historical simulations in CMIP6. *Geosci. Model Dev.* **10**, 4405–4417 (2017).

34. Horn, M. A., Matthies, C., Küsel, K., Schramm, A. & Drake, H. L. Hydrogenotrophic methanogenesis by moderately acid-tolerant methanogens of a methane-emitting acidic peat. *Appl. Environ. Microbiol.* **69**, 74–83 (2003).

35. McCalley, C. K. et al. Methane dynamics regulated by microbial community response to permafrost thaw. *Nature* **514**, 478–481 (2014).

36. Le Mer, J. & Roger, P. Production, oxidation, emission and consumption of methane by soils: a review. *Eur. J. Soil Biol.* **37**, 25–50 (2001).

37. Kelly, C. A., Carolina, N., Hill, C., Dise, B. & Martens, C. S. Temporal variations in the stable carbon isotopic composition of methane emitted from Minnesota peatlands. *Global Biogeochem. Cycles* **6**, 263–269 (1992).

38. Burke, R. A., Barber, T. R. & Sackett, W. M. Methane flux and stable hydrogen and carbon isotope composition of sedimentary methane from the Florida Everglades. *Global Biogeochem. Cycles* **2**, 329–340 (1988).

39. Miller, J. B. & Tans, P. P. Calculating isotopic fractionation from atmospheric measurements at various scales. *Tellus Ser. B Chem. Phys. Meteorol.* **55**, 207–214 (2003).

40. Chanton, J. P., Fields, D. & Hines, M. E. Controls on the hydrogen isotopic composition of biogenic methane from high-latitude terrestrial wetlands. *J. Geophys. Res. Biogeosci.* **111**, 1–9 (2006).

41. Tang, J., Zhuang, Q., Shannon, R. D. & White, J. R. Quantifying wetland methane emissions with process-based models of different complexities. *Biogeosciences* **7**, 3817–3837 (2010).

42. Zhu, X. et al. Rising methane emissions in response to climate change in Northern Eurasia during the 21st century. *Environ. Res. Lett.* **6**, 45211 (2011).

43. Lu, X. & Zhuang, Q. Modeling methane emissions from the Alaskan Yukon River basin, 1986–2005, by coupling a large-scale hydrological model and a process-based methane model. *J. Geophys. Res. Biogeosci.* **117**, G2 (2012).

44. Jin, Z., Zhuang, Q., He, J.-S., Zhu, X. & Song, W. Net exchanges of methane and carbon dioxide on the Qinghai-Tibetan Plateau from 1979 to 2100. *Environ. Res. Lett.* **10**, 85007 (2015).

45. Gilks, W. R., Richardson, S. & Spiegelhalter, D. *Markov chain Monte Carlo in practice*. (CRC press, 1995).

46. Conrad, R. Quantification of methanogenic pathways using stable carbon isotopic signatures: A review and a proposal. *Org. Geochem.* **36**, 739–752 (2005).

47. Blaser, M. & Conrad, R. Stable carbon isotope fractionation as tracer of carbon cycling in anoxic soil ecosystems. *Curr. Opin. Biotechnol.* **41**, 122–129 (2016).

48. Poulter, B. et al. Global wetland contribution to 2000–2012 atmospheric methane growth rate dynamics. *Environ. Res. Lett.* **12**, 094013 (2017).

49. Matthews, E. & Fung, I. Methane emission from natural wetlands: Global distribution, area, and environmental characteristics of sources. *Global Biogeochem. Cycles* **1**, 61–86 (1987).

50. Nisbet, E. G., Dlugokencky, E. J. & Bousquet, P. Methane on the rise—again. *Science* **343**, 493–495 (2014).

51. Hossaini, R. et al. A global model of tropospheric chlorine chemistry: Organic versus inorganic sources and impact on methane oxidation. *J. Geophys. Res. Atmos.* **121**, 14–271 (2016).

52. Saueressig, G. et al. Carbon 13 and D kinetic isotope effects in the reactions of CH₄ with O (¹D) and OH: new laboratory measurements and their implications for the isotopic composition of stratospheric methane. *J. Geophys. Res. Atmos.* **106**, 23127–23138 (2001).

53. Stein, T. Increase in atmospheric methane set another record during 2021. *NOAA Res. News* (2022).

54. Dlugokencky, E. Trends in atmospheric methane. *NOAA/ESRL https://gml.noaa.gov/ccgg/trends_ch4* (2022).
55. Stein, T. New analysis shows microbial sources fueling rise of atmospheric methane. NOAA Research News. https://research.noaa.gov/article/ArtMID/587/ArticleID/2769/New-analysis-shows-microbial-sources-fueling-rise-of-atmospheric-methane (2021).

56. Nisbet, E. G. et al. Very strong atmospheric methane growth in the 4 years 2014–2017: implications for the paris agreement. *Global Biogeochem. Cycles* **33**, 318–342 (2019).

57. Qiu, C. et al. A strong mitigation scenario maintains climate neutrality of northern peatlands. *One Earth* **5**, 86–97 (2022).

58. Schuur, E. A. G. et al. The effect of permafrost thaw on old carbon release and net carbon exchange from tundra. *Nature* **459**, 556–559 (2009).

59. Sullivan, T. D. et al. Influence of permafrost thaw on an extreme geologic methane seep. *Permafr. Periglac. Process* **32**, 484–502 (2021).

60. Walter Anthony, K. M., Anthony, P., Grosse, G. & Chanton, J. Geologic methane seeps along boundaries of Arctic permafrost thaw and melting glaciers. *Nat. Geosci.* **5**, 419–426 (2012).

61. Corbett, J. E. et al. Partitioning pathways of CO₂ production in peatlands with stable carbon isotopes. *Biogeochemistry* **114**, 327–340 (2013).

62. Deng, J. et al. Adding stable carbon isotopes improves model representation of the role of microbial communities in peatland methane cycling. *J. Adv. Model. Earth Syst.* **9**, 1412–1430 (2017).

63. Hodgkins, S. B. et al. Changes in peat chemistry associated with permafrost thaw increase greenhouse gas production. *Proc. Natl. Acad. Sci. USA* **111**, 5819–5824 (2014).

64. Pangala, S. R. et al. Large emissions from floodplain trees close the Amazon methane budget. *Nature* **552**, 230–234 (2017).

65. Barba, J. et al. Methane emissions from tree stems: a new frontier in the global carbon cycle. *New Phytol.* **222**, 18–28 (2019).

66. Rosentreter, J. A. et al. Half of global methane emissions come from highly variable aquatic ecosystem sources. *Nat. Geosci.* **14**, 225–230 (2021).

67. Zhuang, Q. et al. Modeling soil thermal and carbon dynamics of a fire chronosequence in interior Alaska. *J. Geophys. Res. D Atmos.* **108**, (2003).

68. Walter, B. P. & Heimann, M. A process-based, climate-sensitive model to derive methane emissions from natural wetlands: Application to five wetland sites, sensitivity to model parameters, and climate. *Global Biogeochem. Cycles* **14**, 745–765 (2000).

69. Wingate, L. et al. Photosynthetic carbon isotope discrimination and its relationship to the carbon isotope signals of stem, soil and ecosystem respiration. *New Phytol.* **188**, 576–589 (2010).

70. Monson, K. D. & Hayes, J. M. Biosynthetic control of the natural abundance of carbon 13 at specific positions within fatty acids in Escherichia coli. Evidence regarding the coupling of fatty acid and phospholipid synthesis. *J. Biol. Chem.* **255**, 11435–11441 (1980).

71. Mahieu, K., De Visscher, A., Vanrolleghem, P. A. & Van Cleemput, O. Modelling of stable isotope fractionation by methane oxidation and diffusion in landfill cover soils. *Waste Manag.* **28**, 1535–1542 (2008).

72. Preuss, I., Knoblauch, C., Gebert, J. & Pfeiffer, E. M. Improved quantification of microbial CH₄ oxidation efficiency in arctic wetland soils using carbon isotope fractionation. *Biogeosciences* **10**, 2539–2552 (2013).

73. Throckmorton, H. M. et al. Pathways and transformations of dissolved methane and dissolved inorganic carbon in Arctic tundra watersheds: Evidence from analysis of stable isotopes. *Global Biogeochem. Cycles* **29**, 1893–1910 (2015).

74. Hayes, J. M. An introduction to isotopic calculations. *Woods Hole Oceanogr. Inst.* 1–10 (2004).

75. Burke, R. A. Jr, Barber, T. R. & Sackett, W. M. Seasonal variations of stable hydrogen and carbon isotope ratios of methane in subtropical freshwater sediments. *Global Biogeochem. Cycles* **6**, 125–138 (1992).

76. Harris, I., Jones, P. D., Osborn, T. J. & Lister, D. H. Updated high-resolution grids of monthly climatic observations - the CRU TS3.10 Dataset. *Int. J. Climatol.* **34**, 623–642 (2014).

77. Duan, Q. Y., Gupta, V. K. & Sorooshian, S. Shuffled complex evolution approach for effective and efficient global minimization. *J. Optim. Theory Appl.* **76**, 501–521 (1993).

78. Melillo, J. M. et al. Global climate change and terrestrial net primary production. *Nature* **363**, 234 (1993).

79. Global Soil Data Task. 2000. Global Soil Data Products CD-ROM (IGBP-DIS). CD-ROM. International Geosphere-Biosphere Programme, Data and Information System, Potsdam, Germany. Available from Oak Ridge National Laboratory Distributed Active Archive Center, Oak Ridge, Tennessee, U.S.A. [http://www.daac.ornl.gov].

80. Zhuang, Q. et al. Carbon cycling in extratropical terrestrial ecosystems of the Northern Hemisphere during the 20th century: A modeling analysis of the

influences of soil thermal dynamics. *Tellus Ser. B Chem. Phys. Meteorol.* **55**, 751–776 (2003).

81. Myneni, R. B. et al. Global products of vegetation leaf area and fraction absorbed PAR from year one of MODIS data. *Remote Sens. Environ.* **83**, 214–231 (2002).

82. Melton, J. R. et al. Present state of global wetland extent and wetland methane modelling: conclusions from a model inter-comparison project (WETCHIMP). *Biogeosciences* **10**, 753–788 (2013).

83. Miller, S. M. et al. A multiyear estimate of methane fluxes in Alaska from CARVE atmospheric observations. *Global Biogeochem. Cycles* **30**, 1441–1453 (2016).

84. Chang, R. Y. W. et al. Methane emissions from Alaska in 2012 from CARVE airborne observations. *Proc. Natl. Acad. Sci. USA* **111**, 16694–16699 (2014).

85. Werf, G. R. et al. Global fire emissions estimates during 1997–2016. *Earth Syst. Sci. Data* **9**, 697–720 (2017).

86. Janssens-Maenhout, G. et al. EDGAR v4.3.2 Global Atlas of the three major Greenhouse Gas Emissions for the period 1970–2012. *Earth Syst. Sci. Data* **11**, 959–1002 (2019).

87. Etiope, G. & Schwietzke, S. Global geological methane emissions: an update of top-down and bottom-up estimates. *Elem. Sci. Anth.* **7**, 47 (2019).

88. Bergamaschi, P., et al. Satellite chartography of atmospheric methane from SCIAMACHY on board ENVISAT: 2. Evaluation based on inverse model simulations. *J. Geophys. Res. Atmos.* **112**, https://doi.org/10.1029/2006JD007268 (2007).

89. Sherwood, O. A. et al. Global δ¹³C-CH₄ Source Signature Inventory. https://doi.org/10.15138/qn55-e011 (2020).

90. Basu, S. et al. Global CO₂ fluxes estimated from GOSAT retrievals of total column CO₂. *Atmos. Chem. Phys. Discuss.* **13**, 8695–8717 (2013).

91. Meirink, J. F., Bergamaschi, P. & Krol, M. C. Four-dimensional variational data assimilation for inverse modelling of atmospheric methane emissions: method and comparison with synthesis inversion. *Atmos. Chem. Phys.* **8**, 6341–6353 (2008).

92. Tans, P. P. A note on isotopic ratios and the global atmospheric methane budget. *Global Biogeochem. Cycles* **11**, 77–81 (1997).

93. Masarie, K. A. & Tans, P. P. Extension and integration of atmospheric carbon dioxide data into a globally consistent measurement record. *J. Geophys. Res. Atmos.* **100**, 11593–11610 (1995).

94. Prigent, C., Jimenez, C. & Bousquet, P. Satellite-derived global surface water extent and dynamics over the last 25 years (GIEMS-2). *J. Geophys. Res. Atmos.* **125**, e2019JD030711 (2020).

## Acknowledgements

This work was supported by NASA Earth and Space Science Fellowship Program (#80NSSC17K0368 P00001) and Interdisciplinary Research in Earth Science (#NNX17AK20G). We thank Carmody K. McCalley for providing data and John Mund for technical support.

## Author contributions

Y.O., Q.Z., and X.L. conceived the study. Y.O., Q.Z., L.L., and L.R.W. built the model. E.J.D., S.E.M., J.B.M., S.S., and P.C. provided unpublished or raw data. Y.O. conducted model runs. S.B., L.B., P.T., and J.P.C., and all other authors contributed to data interpretation and preparation of paper text.

## Competing interests

The authors declare no competing interests.

## Additional information

Supplementary information The online version contains supplementary material available at https://doi.org/10.1038/s43247-022-00488-5.

Correspondence and requests for materials should be addressed to Youmi Oh or Qianlai Zhuang.

Peer review information *Communications Earth & Environment* thanks the anonymous reviewers for their contribution to the peer review of this work. Primary handling editors: Leiyi Chen and Clare Davis.

Reprints and permission information is available at http://www.nature.com/reprints

Publisher's note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
<figure><img src="image_7.png" /></figure>

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2022