Introduction

In the years 2020 and 2021, the methane growth rate (MGR) in the atmosphere reached 15.2 ± 0.5 and 17.8 ± 0.5 parts per billion per year (ppb yr−1) respectively, hitting record highs since systematic measurements started in the early 1980s by NOAA’s Global Monitoring Laboratory (GML)1. The unprecedented methane growth during 2020 and 2021 coincided with the reduced human activities and pollutant emissions during COVID-19 lockdowns and the gradual recovery afterwards2,3,4,5, together with the occurrence of a moderate and prolonged La Niña event6, which offers a unique opportunity to examine the drivers of methane variabilities on a year-to-year basis.

Both process-based studies of sources (“bottom-up” estimates) and atmospheric-based inverse analyses (“top-down” estimates) pointed to pronounced emission growth in 2020 compared to 2019, arising from tropical and northern sources, likely driven by enhanced wetland emissions7,8,9,10, the main source component of natural methane emissions. From 2019 to 2020, the increase of MGR was earlier and larger in the Tropics and Northern high-latitudes than in the Southern extra-tropics from NOAA GML marine boundary layer sites (Figs. 1b and S1; see Methods)11, consistent with an activation of wetland emissions. The exceptionally high methane growth persisted in 2021, and prevailed over most latitude bands (Fig. 1a, b). Observations of total column methane concentrations (XCH4) by Greenhouse Gases Observing SATellite (GOSAT), whether obtained from the National Institute for Environmental Studies (Japan) full physics retrievals (hereafter “GSNIES”)12 or from the University of Leicester proxy retrievals (hereafter “GSUoL”)13, confirmed the unexpected methane surges in 2020 and 2021 (Fig. 1c, d). Note that both GOSAT retrievals showed a larger MGR increase in 2020 over the Southern extra-tropics than that observed from surface network, with GSUoL exhibiting higher MGR globally and in the Tropics as well (Figs. 1 and S1).

Fig. 1: Variations of atmospheric methane growth rates between 2010 and 2021.
figure 1

a–d Methane (CH4) growth rates were derived from zonally averaged observations of NOAA Global Monitor Laboratory (NOAA/GML) marine boundary layer (MBL) sites11, GOSAT NIES retrievals12 or GOSAT UoL retrievals13 of total column methane concentrations (XCH4), following the curve-fitting routines of ref. 71. For c, d, results are not shown north of 50°N or south of 50°S due to data gaps of GOSAT retrievals over these regions during winter. Source data are provided as a Source Data file.

The zonally-averaged changes in MGR across latitude bands reveal the integrated variations of regional sources and sinks, atmospheric transport, and removal by the hydroxyl radical (OH). In this study, we use a global atmospheric inversion system (PYVAR-LMDZ-SACS; see Methods)14,15 that assimilated either surface or satellite-based CH4 observations to infer the spatiotemporal patterns in flux changes from 2019 to 2021. An ensemble of six inversions is performed, with different assimilated observation datasets (surface network, GSNIES or GSUoL) and transport model physical parameterizations (the “classic” and “advanced” versions) (Table S1; see Methods). This allows us to test the consistency of flux change patterns inferred from different types of measurements while accounting for the uncertainties due to imperfect representation of atmospheric mixing. Surface networks offer a good coverage of northern mid-to-high latitudes (especially over Europe and North America), whereas satellite data have improved data densities from 60°S to 60°N, particularly a better coverage of the tropics than surface networks (Figs. S2 and S3). We prescribe changes in OH concentration fields from a full chemistry transport model (LMDZ-INCA)16,17, driven by interannually varying meteorology18 with natural and anthropogenic emissions of NOx, CO and hydrocarbons updated to 2021 (Fig. S4; Table S2; see Methods)19,20,21. The ensemble of atmospheric inversions reveals anomalous and persistent emission surges from inundated areas in tropical Africa and Asia during 2020–2021, linked to water table rises and La Niña conditions. Meanwhile, we also quantify variations of methane emissions using bottom-up inventories for anthropogenic and fire sources, and process-based biogeochemical models for wetland emissions (see Methods). We find that the strong emission enhancements over tropical inundated areas are not captured by current process-based wetland models, likely pointing to models’ limitations in accurately representing dynamics of tropical wetland extents and processes driving tropical wetland methane emissions.

Results and Discussions

Variations of methane emissions from atmospheric inversions

All members of the ensemble of six inversions showed global increases in surface CH4 emissions by an average of 20.3±9.9 Tg CH4 yr-1 and 24.8±3.1 Tg CH4 yr-1 in 2020 and 2021 respectively, compared to 2019 (uncertainty being one standard deviation of the ensemble). A large portion of this surge in global emissions was accounted for by the northern tropics (0°–30°N), which contributed about 80% (16.2±8.3 Tg CH4 yr-1) and 95% (23.2±4.0 Tg CH4 yr-1) of the global increases in 2020 and 2021, respectively (Figs. 2a, b and S5; Table S3). Higher emission increases were found over tropical Africa and Southeast Asia in both years, according to all six inversions (Figs. 2a, b and S5–S7). Overall, most of the regions with strong and consistent emission changes overlapped with major wetland complexes and inundated areas. These emission anomalies aligned well with changes in liquid water equivalent (LWE) from the Gravity Recovery and Climate Experiment Follow-On (GRACE-FO) satellites (Figs. 2a–d and S5c–f). LWE is a proxy for water-table depth and water stored in wetland systems22.

Fig. 2: Methane emission anomalies in 2020 and 2021 relative to 2019 inferred from six inversions in relation to changes in satellite observed liquid water equivalent (LWE).
figure 2

a, b Spatial patterns of the mean CH4 emission anomalies (ΔCH4 emissions) averaged over the ensemble of six top-down (TD) inversions. The shaded areas indicate that posterior fluxes from all six inversions have the same changing direction. The inset bar plots summarize the net emission changes at the global scale and for four latitude bands. Error bars represent one standard deviation of the six inversions. c, d Spatial patterns of changes in LWE (ΔLWE) from the GRACE-FO satellites. e, g Scatterplot of CH4 emission anomalies versus changes in GRACE-FO LWE across 20 major wetland regions defined based on the regularly flooded wetland map23 in f. Each circle represents a major wetland region, with the size scaling with the magnitude of the region’s posterior methane emission in 2019 (unit: Tg CH4 yr-1). The color of the circle corresponds to the color code for each wetland region in f. Source data are provided as a Source Data file.

Focusing on the 20 major wetland regions of the globe that represent about 60% of all wetland areas based on the map of ref. 23, a significant correlation was noted between top-down estimates of emission anomalies and changes in LWE from GRACE-FO (Fig. 2e–g). The correlation with temperature or precipitation changes was not statistically significant across regions (Fig. S8), although these two factors could potentially impact emission variations between years for individual wetland regions (Table S4). Among the 20 wetland regions, the Niger River basin, the Congo basin, the Sudd swamp, the Ganges floodplains and the Southeast Asian deltas in the tropics, and the Hudson Bay lowlands in the boreal region, exhibited persistent emission enhancements in both 2020 and 2021. Emission increases over each of these regions exceeded by at least 1.5–2 times the interannual variability (1.5–2σ) of methane emissions during 2010–2019 (Table S5), based on the results of inversions with settings comparable to those in our study but covering a longer period24. These six wetland regions together contributed around 70% (14.1±4.2 Tg CH4 yr-1) and 60% (14.9±2.6 Tg CH4 yr-1) of the global emission increases for 2020 and 2021, respectively (Figs. 23e, f and S5; Table S5). The synchronous emission increases over these wetland regions were paralleled by overall wetter conditions (increased LWE or precipitation), and additionally by warming conditions for specific regions (e.g., Hudson Bay lowlands; Table S4).

Fig. 3: Top-down versus bottom-up estimates of methane emission anomalies in 2020 and 2021 relative to 2019.
figure 3

a, b Mean methane emission anomalies (ΔCH4 emissions) of four latitude bands derived from the ensemble of six top-down (TD) inversions and bottom-up (BU) estimates. The black dots represent the net global emission changes relative to 2019. c, d Spatial patterns of bottom-up CH4 emission anomalies summed up from process-based wetland models, and inventories of anthropogenic and fire emissions. The color scale is the same as for the top-down CH4 emission anomalies in Fig. 2a, b. The inset bar plots summarize the net emission changes at the global scale and for four latitude bands. e, f Top-down versus bottom-up estimates of methane emission anomalies for eight tropical inundated areas. The delineation of each inundated area is shown in Fig. 2f and Fig. S5i. The open circle indicates two times the interannual variability (+2σ or –2σ) of methane emissions during 2010–2019 derived from a previous study using the inversion system PYVAR-LMDZ-SACS and GSUoL as constraints for CH4 (ref. 24). Error bars in all panels denote one standard deviation of the methane emission anomalies from the ensemble of top-down or bottom-up estimates. Source data are provided as a Source Data file.

Substantial emission increases were also found over the western Siberian lowlands in 2020 and the Amazon Basin and the Orinoco floodplain in 2021, but reduced emissions were seen in the other years over these regions (i.e., the western Siberian lowlands in 2021; the Amazon Basin and the Orinoco floodplain in 2020), aligning with corresponding variations in temperature and/or LWE (Figs. 23e, f, and S5; Table S4). For the Amazon Basin in particular, the enhanced emissions in 2021 were associated with increased floodplain inundation and the record high water levels25, linked to intense rainfall over the northern Amazonia driven by the La Niña conditions, with strengthened Walker circulation and deep convection26,27. Conversely, the Pantanal and the Paraná floodplains in central and southeastern South America exhibited emission reduction in both 2020 and 2021, due to continued drier conditions and lower water tables (Figs. 23e, f and S5; Table S4). These two regions experienced severe and prolonged drought events since 2019, leading to record low water levels over the past decades28,29,30 and therefore persistent reduction in CH4 emissions.

The overall strong emission increases over tropical inundated areas were coincident with the occurrence of a prolonged La Niña event since 2020. La Niña is the cold phase of the El Niño—Southern Oscillation (ENSO) cycles—the dominant mode of interannual climate variability impacting weather patterns and ecosystem dynamics across the globe. La Niña is generally associated with a wetter climate and above-normal precipitation over tropical land areas31, which correlates with higher wetland methane emissions through expanded flooded areas. Previous studies suggested a broad anti-correlation between wetland CH4 emissions and ENSO at the global and pan-tropical scales, with less emissions during strong El Niño events and more emissions during strong La Niña events (Fig. S9; also see refs. 32,33,34,35). At the regional scales, however, the sign, magnitude and phasing of wetland responses to ENSO could vary32,34,35 depending on the spatiotemporal variability in climate anomalies during specific ENSO episodes36, which complicates the relationship between ENSO, wetland CH4 emissions, and CH4 growth rates37. Our results covering the most recent La Niña episode confirm the overall enhanced CH4 emissions over inundated areas during La Niña conditions, and further demonstrate the diverse response patterns across regions and between years that could be linked to the complex evolution of ENSO-related climate anomalies.

Bottom-up estimates of methane emission changes

Bottom-up inventories and process-based wetland emission models independently corroborate increased emissions over the northern tropics during 2020–2021, albeit with much smaller magnitudes than inversions (Figs. 2, 3 and S10; Table S3). Globally, the net emission changes were −0.7 ± 4.0 Tg CH4 yr−1 and 6.9 ± 4.1 Tg CH4 yr−1 in 2020 and 2021 relative to 2019 based on bottom-up methodologies, with emission increases of 3.2 ± 1.3 Tg CH4 yr−1 and 4.1 ± 1.9 Tg CH4 yr−1 estimated for the northern tropics (Figs. 3a, b and S10; Table S3). A breakdown into different processes showed that the anthropogenic CH4 emissions in 2020 and 2021 were higher than the 2019 level by 0.5 Tg CH4 yr−1 and 4.7 Tg CH4 yr−1, respectively (Fig. S10; Table S3). While fossil fuel CH4 emissions were reduced during COVID-19 lockdowns and rebounded partially afterwards (by −2.7 Tg CH4 yr−1 and −0.5 Tg CH4 yr−1 for 2020 and 2021 compared to 2019), the CH4 biogenic emissions from sectors of agriculture and waste treatment continued to increase (by 3.1 Tg CH4 yr−1 and 5.2 Tg CH4 yr−1 for 2020 and 2021 compared to 2019), leading to a net increase of anthropogenic emissions. The global fire CH4 emissions declined in 2020 by 6.5 Tg CH4 yr−1 relative to 2019, mainly contributed by reduced fire emissions of 5.1 Tg CH4 yr−1 in the southern tropics (30°S–0°)21. A similar reduction of fire emissions occurred again in the southern tropics in 2021, but the extreme fires in boreal North America and eastern Siberia during the hotter and drier summertime38 led to an increase in emissions of 4.7 Tg CH4 yr−1 in the northern extra-tropics and therefore a net global emission reduction of only 1.8 Tg CH4 yr-1 relative to 2019 (Figs. 3a, b and S10; Table S3). For wetlands, the results of process-based wetland model simulations from ORCHIDEE-MICT and LPJ-wsl driven by four different climate forcings (see Methods) reported an increase in global wetland emissions by 5.3 ± 4.0 Tg CH4 yr−1 in 2020 compared to 2019, which was dominated by enhanced wetland emissions in the northern extra-tropics and tropics as a result of warmer and wetter climate8. The updated simulations showed slightly smaller wetland emission increases by 4.0 ± 4.1 Tg CH4 yr−1 in 2021 compared to 2019, with similar latitudinal patterns (Figs. 3a, b, S10 and S11; Table S3).

Given differences of simulated wetland emissions by process-based biogeochemical models (Fig. S11), the global total emission increases from the bottom-up approach were much lower than those from inversions. When the maximum wetland emission increases are considered (i.e., 13.4 Tg CH4 yr−1 and 11.5 Tg CH4 yr−1 in 2020 and 2021 by ORCHIDEE-RFW with MERRA2; Fig. S11), the total emission increases amounted to 7.4 Tg CH4 yr−1 and 14.4 Tg CH4 yr−1 in 2020 and 2021, respectively, still 12.9 Tg CH4 yr−1 and 10.4 Tg CH4 yr−1 below the mean of the inversion ensemble estimates. While there was a general agreement in the broad spatial patterns of emission changes between the bottom-up and top-down approaches, large underestimation by the bottom-up approach was found over the tropics, especially in tropical Africa and Asia (Figs. 2, 3 and S10). Specifically, for the Niger River basin, the Congo basin, the Sudd swamp, the Ganges floodplains and the Southeast Asian deltas, the total emission increases from the bottom-up approach only accounted for 11% and 5% of the estimates from top-down inversions for 2020 and 2021, respectively (Fig. 3e, f; Table S5).

The lower bottom-up estimates of emission increases over tropical Africa and Asia inundated areas suggest that wetland biogeochemical models may misrepresent key processes in wetland dynamics and associated changes in CH4 emissions. A close look into the Sudd swamp, the wetland emission hotspot in the eastern Africa (3°N–17°N, 25°E–40°E; see also refs. 39,40,41) showed that the annual variations in top-down CH4 emissions and GRACE-FO LWE anomalies were highly correlated (Fig. S12), implying a strong impact of water-table depth on seasonal emissions over tropical inundated areas (Fig. S13; see also refs. 22,39,42). The annual peak of LWE occurs during September–October, about one month before the annual peak of emissions from top-down inversions. The seasonal peak of emissions during September–November in Sudd (about 0.3 Tg CH4 month−1 higher relative to April–June) is missed by most wetland model simulations (giving only a non-significant rise of 0.02 Tg CH4 month−1 during September–November relative to April–June; Fig. S12). This could be partly explained by the simulation of too small wetland extents of the Sudd swamp and weak intra- and inter-annual variations, as shown by the comparison between the simulated flooded areas and high-resolution observations from the CYGNSS satellites43,44. The too small seasonal changes in wetland extent and associated CH4 emissions were also identified previously in other process-based or data-driven wetland models over tropical regions33,41,45,46, resulting in smaller estimates of year-to-year emission anomalies. Besides, large uncertainties remain in the model representation of CH4 production, oxidation and vegetation-mediated transport processes for tropical wetlands47,48, where direct flux measurements are sparse49,50. For example, field studies have demonstrated the importance of vegetation-mediated CH4 fluxes in tropical wetlands even when the water table is below the soil surface51. These processes have been considered in recent development of process-based wetland models but not calibrated for tropical wetlands due to lack of field measurements52. More observations of tropical wetland hydrology and associated CH4 emissions are thus needed to calibrate and constrain model parameterizations, thereby improving models’ predictive capacity. It should be noted that other emissions from anthropogenic sources and biomass burning may co-exist in the inundated areas (Figs. S10 and S12). Estimation of these sources could also be uncertain and further contribute to the discrepancies between bottom-up and top-down estimates of emission anomalies.

In summary, the record high CH4 growth rates in 2020 and 2021 revealed heightened emissions from inundated areas in the tropical and boreal regions. Strong and persistent emission surges were found simultaneously over the Niger River basin, the Congo basin, the Sudd swamp, the Ganges floodplains, the Southeast Asian deltas, and Hudson Bay lowlands, coincident with elevated groundwater and warming in the north, and potentially linked to La Niña conditions prevailing since 2020. Our main findings on heightened emissions from both tropical and boreal wetlands, along with evidence from other bottom-up or top-down studies10,39,53,54,55,56, suggest recent intensification of wetland methane emissions and probable strong positive wetland climate feedback57. In a future with warming climate (unavoidable in the Arctic) and possibly increasing occurrences of extreme or prolonged La Niña events58,59, the synchronous emission rise across boreal and tropical wetlands, as seen in 2020 and 2021, may occur more frequently and has the potential to accelerate atmospheric CH4 growth. This would challenge the commitment of the Paris Agreement to limit global warming, and stresses the urgency of greater reduction in anthropogenic emissions to achieve climate mitigation goals57,60. Our study also adds key ‘natural experiments’ from La Niña conditions in 2020 and 2021 to better understand the complex spatial-temporal variations in wetland emissions linked to ENSO. A systematic evaluation of wetland responses and their asymmetry between El Niño and La Niña conditions will help identify fingerprint of ENSO-induced climate variability in wetland dynamics, and may eventually contribute to an improved characterization of future variations in wetland CH4 emissions and atmospheric CH4 growth rates. Current process-based wetland models fail to capture the strong emission surges over the tropics revealed by top-down inversions. Improved wetland extent mapping and modeling, informed by sustained and enhanced observations on wetland hydrology and biogeochemistry, will help constrain tropical wetland emissions under changing climate and reconcile CH4 budget estimation between bottom-up and top-down approaches.

Methods

We followed the methodologies described in ref. 8 to combine both top-down and bottom-up approaches for a synthesis study of recent methane growth during 2020–2021. The analyses of methane emission changes were extended to 2021 on the basis of ref. 8 using similar data sources and modeling tools. We further included satellite-based CH4 observations and flux inversions in addition to the analyses derived from surface CH4 networks, which allows for intercomparison of the emission change patterns that are informed by different datasets of atmospheric CH4 observations.

Atmospheric observations

For surface CH4 observations, in-situ continuous and flask-air CH4 measurements from a total of 121 stations for the inversion analyses were included (Fig. S2), most of which are operated and maintained by the NOAA61 and Integrated Carbon Observation System62. Observations from other networks, including those from Environment and Climate Change Canada (ECCC), Advanced Global Atmospheric Gases Experiment (AGAGE), Japan Meteorological Agency (JMA), etc., were obtained from the World Data Centre for Greenhouse Gases (https://gaw.kishou.go.jp) and the Global Environmental Database (https://db.cger.nies.go.jp). All observations are reported on or linked to the WMOX2004 calibration scale.

For satellite CH4 observations, we used two retrievals of GOSAT XCH4 provided by National Institute for Environmental Studies in Japan and the University of Leicester in the UK (denoted as “GSNIES” and “GSUoL” respectively). Launched by the Japan Aerospace Exploration Agency (JAXA) in early 2009, GOSAT achieves a global coverage every 3 days with a swath of 750 km and a ground pixel with a diameter of approximately 10.5 km at nadir. The Thermal And Near-infrared Sensor for Carbon Observation – Fourier Transform Spectrometer (TANSO-FTS) onboard enables the measurements of column-averaged dry-air CO2 and CH4 mole fractions by solar backscatter in the shortwave infrared (SWIR) with near-unit sensitivity across the air column down to the surface63,64. The two GOSAT XCH4 products used here differ in their algorithms to treat the scattering-induced issues in the retrieval of total column concentrations from spectral data. The GSNIES XCH4 retrieval was produced using a full-physics algorithm to infer CH4 column together with physical scattering properties of the atmosphere65,66. Alternatively, the GSUoL XCH4 retrieval employed a proxy algorithm that simultaneously retrieves CH4 and CO2 columns using the absorption features around the wavelength of 1.6 μm to minimize the scattering effect on the retrieval67,68. While the two conceptually different approaches have their own advantages and disadvantages69, the proxy retrieval is less sensitive to aerosol distribution and instrumental issues than the full-physics retrieval, therefore has much higher data density over geographic regions with substantial aerosol loading, such as in the tropics (Fig. S3). In this study, we used the GSNIES XCH4 retrieval version 2.95/2.9612 and the GSUoL XCH4 retrieval version 9.013, which are bias-corrected and in good agreement with ground-based XCH4 measurements from the Total Column Carbon Observing Network (TCCON) and aircraft-based CH4 profile measurements. These two retrievals have been widely used in global or regional methane inverse modeling to study recent trends and interannual variabilities7,9,54,55,70. Note that only retrievals over land (except those over Greenland; Fig. S3) were assimilated in our inversions in order to avoid potential retrieval biases between nadir and glint viewing modes.

To calculate atmospheric CH4 growth rates in Figs. 1 and S1, for surface observations, we used zonally-averaged marine boundary layer (MBL) references for CH4 constructed by NOAA’s Global Monitoring Laboratory (NOAA/GML) using measurements of weekly air samples from a subset of sites in the NOAA Cooperative Global Air Sampling Network11. Only sites that measure background atmospheric compositions are considered, typically at remote marine sea level locations with prevailing onshore winds (see the site map at https://gml.noaa.gov/ccgg/mbl/map.php?param=CH4). For GOSAT XCH4 observations, we used daily means of all valid land retrievals per 10° latitude band between 50°N and 50°S for subsequent growth rate calculation. Note that GOSAT XCH4 retrievals north of 50°N or south of 50°S were discarded for the calculation of growth rates due to data gaps over these regions during winter, but these data were included in atmospheric inversions (Fig. S3). The smoothed CH4 growth rates for each latitude band shown in Fig. 1 and S1 were extracted from time series of the zonally-averaged MBL CH4 references or GOSAT XCH4 observations over the period 2010–2021 following the curve fitting procedures of ref. 71.

Atmospheric 3D inversions

We used a variational Bayesian inversion system PYVAR-LMDZ-SACS (PYthon VARiational—Laboratoire de Météorologie Dynamique model with Zooming capability—Simplified Atmospheric Chemistry System)14,15 to optimize weekly CH4 surface fluxes at a spatial resolution of 1.9° in latitude by 3.75° in longitude over the period 2019–2021. An ensemble of six inversions was performed (Table S1), using three different observation datasets described above as constraints, combined with two physical parameterizations for the transport model Laboratoire de Météorologie Dynamique with zooming capability (LMDZ), the atmospheric component of the coupled IPSL climate model participating in IPCC Assessment Reports (AR). These setups allowed us to explore consistency of the emission change patterns informed by different observation datasets while accounting for some of the uncertainties in atmospheric transport. The two physical parameterizations, denoted here as the “classic” and “advanced” versions, represent two development stages of LMDZ for IPCC AR3 and AR672,73,74. The “classic” AR3 version72 uses the vertical diffusion scheme of ref. 75 to represent the turbulent transport in the boundary layer and the scheme of ref. 76 to parameterize deep convection. The “advanced” AR6 version74 combines the vertical diffusion scheme of ref. 77 and the thermal plume model by ref. 78 to simulate the atmospheric mixing in the boundary layer, and the deep convection is represented using the scheme of ref. 79 coupled with the parameterization of cold pools developed in ref. 80 While the “advanced” version showed overall improved representation of boundary layer mixing and large-scale atmospheric transport81,82,83, which would benefit trace gas transport simulations and inversions despite its comparatively larger computational costs, the “classic” version or its physical parameterization schemes are still widely used in the scientific community for methane studies (see Table S8 and S10 in ref. 84). A previous study85 showed that changing physical parameterizations would have small impact on the inverted methane emissions at the global scale (around 1%), but could lead to significant differences in the north–south gradient of emissions and the emission partitioning between regions (also see Table S6).

Depending on the observations assimilated and the physical parameterizations used in the inversion system, discrepancies in the derived emission changes do exist among inversions at global or regional scales. The emission growth inferred from surface observations was much lower than those from GOSAT-based inversions for 2020 (Fig. S5a), possibly because surface networks have limited spatial coverage over certain key source regions (e.g., the tropics), thus being blind to the methane growth there (Fig. S2). Among the four GOSAT-based inversions, the ones constrained by GSUoL retrievals always gave 15–30% higher global net emission increases than those constrained by GSNIES retrievals (Fig. S5a, b), consistent with the steeper rise of the methane growth rate seen from GSUoL (Figs. 1 and S1). For several important emitting regions such as eastern China, northern India and southern Africa, the directions of emission changes disagreed among inversions assimilating different observations (Figs. S6 and S7), reflecting uncertainties in flux solution related to sparse data density or GOSAT XCH4 retrieval algorithms. Note that we did not adjust GOSAT XCH4 retrievals in our inversions. We acknowledge the inconsistency between surface-based and GOSAT-based inversions regarding the inferred global total emissions and their latitudinal gradients (Table S6; also see refs. 55,85,86,87), likely due to systematic biases in GOSAT retrievals and/or model representation of the stratospheric methane. Several previous studies applied an empirical bias correction on GOSAT retrievals, making the optimized atmospheric CH4 more consistent with surface observations85,87,88,89. As our study focuses on the year-to-year changes of the optimized methane fluxes, which would be less sensitive to aforementioned systematic biases, we did not apply such an adjustment to GOSAT retrievals in our inversions (also see ref. 55).

Other configurations of the inversions followed the descriptions in ref. 8. The prior CH4 fluxes were built on bottom-up inventories or process-based land surface models for different categories (Table S7), consistent with the priors used for top-down inversions contributing to the new global methane budget (GMB) assessment84, which allows for comparisons with earlier studies and other modeling work within the GMB framework. The OH and O(1D) fields were prescribed from the simulation of a chemistry-climate model LMDZ-INCA (Laboratoire de Météorologie Dynamique model with Zooming capability – INteraction with Chemistry and Aerosols) with a full tropospheric photochemistry scheme16,17. The model was run at the resolution of 1.27° in latitude by 2.5° in longitude, driven by interannually varying horizontal winds from ECMWF ERA5 reanalysis18 and with natural and anthropogenic emissions of NOx, CO and hydrocarbons updated to 2021 (refs. 19,20,21). As countries relaxed and lifted COVID-19 restrictions, global NOx emissions started to rebound in spring 2021, but annual emissions were still lower by 3–4% compared to 2019 (Fig. S4a–c; Table S2). Variations in CO emissions are dominated by biomass burning. While global CO emissions declined by about 12% in 2020 compared to 2019 because of less fire emissions in the Southern Hemisphere, emissions in 2021 were quite similar to the 2019 level. The reduced CO emissions in the Southern Hemisphere were sustained in 2021 but offset by abnormally high boreal fire emissions during summertime (Fig. S4d–f; Table S2)38. Our chemistry transport model LMDZ-INCA simulated a reduction of global tropospheric OH by 3% in 2021 relative to 2019, more than the 1.6–1.8% decrease in 2020 reported by ref. 8. Note that the resulting oxidant fields were not adjusted in the inversions in order to keep the simulated OH changes from LMDZ-INCA. We acknowledge that the reduced atmospheric sink could play an important role in the record high CH4 growth rates during 2020–2021, although the attribution between changes in sources and sinks is still under debate4,7,8,9,90. While the changes in global tropospheric OH were broadly consistent with independent estimates and other studies for 20207,8,9,91, the simulated OH change for 2021 could be more uncertain due to e.g. inaccurate emission forcings92, thereby complicating attribution of the anomalous CH4 growth rates between sources and sinks. However, the dominant contribution from tropical regions to recent CH4 emission surges was confirmed by other top-down studies using different inversion systems and OH configurations7,9.

The preprocessing of surface CH4 observations and the assignment of observation errors were based on the protocol described in ref. 8. For continuous CH4 measurements in particular, we calculated daily afternoon averages (12:00–16:00 local sidereal time, LST) for stations situated below 1000 meters above sea level (a.s.l.) and morning averages (0:00–4:00 LST) for stations located above 1000 meters a.s.l., which were used in subsequent inversion analyses. The observation errors are comprised of measurement errors, representation errors, and model errors. The latter two were approximated using the synoptic variability of CH4 observations at each station, calculated as the residual standard deviation (RSD) of the de-trended and de-seasonalized observations85,93. For GOSAT XCH4 observations, the valid data were averaged into model grids for each time step (30 mins) to create “super-observations”, with the observation errors defined as the retrieval errors reported by the data product plus model errors whose standard deviations were empirically set as 1% (refs. 55,93).

Bottom-up estimates of methane emissions

Anthropogenic methane emissions were compiled from a combination of existing inventories. For the 42 Annex-I countries that report their national greenhouse gas inventories (NGHGIs) to UNFCCC each year, we used the reported anthropogenic methane emissions updated to 2021 from coal mining, oil and gas production, agriculture and waste sectors, respectively (https://unfccc.int/ghg-inventories-annex-i-parties/2023#fn2). For China, anthropogenic methane emissions were computed and updated to 2021 based on the activity data collected from national statistic books94 and specific emission factors at provincial levels95,96. For other countries, emissions from coal mining, oil and gas production, agriculture and waste were obtained from the Emissions Database for Global Atmospheric Research version 7.0 (EDGAR v7.0)97,98, with coal production and livestock data corrected by the activity data from International Energy Agency (IEA)99 and Food and Agriculture Organization of the United Nations (FAO)100. The national total emissions from the 42 Annex-I countries, China and other countries were distributed on 0.1° × 0.1° grid cells based on the spatial patterns of EDGAR v7.0. Note that the change in global anthropogenic methane emissions between 2020 and 2019 slightly differs from that reported by ref. 8, as different versions of EDGAR and IEA data were used in this study to estimate emissions.

The global fire methane emissions were obtained from the Global Fire Emissions Database version 4.1 including small fire burned area21,101. The data set produces monthly gridded burned area and fire emissions at a spatial resolution of 0.25° × 0.25°, based on satellite information on fire activity and vegetation productivity21.

For wetland emissions, we used two process-based wetland emission models (WEMs), ORCHIDEE-MICT102 and LPJ-wsl103, to simulate the global wetland CH4 emissions. Both ORCHIDEE-MICT and LPJ-wsl participated in the Wetland and Wetland CH4 Inter-comparison of Models Project104,105, and have contributed to wetland biogeochemical model simulations of the Global Methane Budget84,106,107. These two models differ in model structures and parameterizations, and represent part of the wide uncertainty range in characterizing CH4 emission processes104,105. An inter-comparison among GMB wetland model ensemble shows that ORCHIDEE-MICT and LPJ-wsl cover a considerable uncertainty range across the GMB wetland model ensemble in simulating global annual wetland CH4 emissions and emission changes from 2019 to 2020, and broadly agree with the interannual variability of model ensemble means (Fig. S15). Based on the simulation protocol in ref. 8, wetland methane emissions were updated to 2021 using these two WEMs with four climate datasets (i.e., CRU TSv4.06108, ERA518, MERRA2109, and MSWEPv2.8110,111) to account for uncertainties in simulating wetland CH4 emissions due to meteorological forcings34. The spatiotemporal dynamics of wetland areas were simulated by a TOPMODEL-based diagnostic model and applied to ORCHIDEE-MICT112,113 and LPJ-wsl114, respectively. For ORCHIDEE-MICT in particular, we utilized two wetland maps to calibrate the parameters in simulating the wetland area dynamics113. The static map of Regularly Flooded Wetlands (RFW)23 was applied for the grid-based calibration of the long-term maximum wetland extent, whereas the Global Inundation Estimate from Multiple Satellites version 2 (GIEMS-2)115 was applied to calibrate the yearly maximum wetland extent for each grid. Combining two wetland maps for ORCHIDEE-MICT and one for LPJ-wsl with four climate datasets, we included a total of twelve simulations of wetland CH4 in our analyses (Fig. S11).