High-resolution inversion of methane emissions over Europe using the Community Inversion Framework and FLEXPART

Constraining methane (CH₄) emissions at high spatial and temporal resolution is critical for accurate European greenhouse gas budgets and mitigation policy. We use the Community Inversion Framework to estimate monthly CH₄ fluxes across Europe (2017–2022) at 0.2°×0.2°, coupling FLEXPART and assimilating observations from 46 in situ stations, including Integrated Carbon Observation System (ICOS) and non-ICOS sites. Prior emissions combine anthropogenic inventories, biomass burning estimates, wetland models, and climatological natural sources. The inversion substantially improves agreement with atmospheric observations (r²=0.87, RMSE=24.4 ppb, $\text{mean bias}=-\mathrm{2.1}$ ppb), performing best at northern European stations. For the European Union countries, together with the UK, Norway, and Switzerland, we estimate total methane emissions of 23.3±2.3 Tg CH₄ yr⁻¹, which is 6.6 % higher than the prior. For these countries, we found an average anthropogenic emissions of 17.6 Tg CH₄ yr⁻¹, with a decreasing trend of 0.3 Tg yr⁻¹. This estimate exceeds the prior by 11 %, EDGARv8 by 4 %, and UNFCCC NGHGI (2023) by 3 %, while remaining consistent with recent studies. Country-level differences are notable, with higher emissions estimated for the Netherlands and Germany, and lower for Romania and Italy. Sectoral changes mainly reflect agricultural increases, alongside reductions in northern wetlands and southern geological sources. Sensitivity tests highlight the influence of spatial spread of emissions distribution assumptions and the importance of dense observational networks for refining regional CH₄ budgets.

1 Introduction

Methane (CH₄) mole fractions in Earth's atmosphere have nearly tripled since 1750, significantly contributing to a 0.6 °C rise in global temperatures since the pre-industrial era (Chen et al., 2022; Shen et al., 2023; Saunois et al., 2020, 2025). The World Meteorological Organization (WMO) has reported record increases in global methane levels from 2020 to 2021, with further rises observed in 2022 (WMO, 2023). These increases pose substantial challenges, such as accelerating climate change, disrupting ecosystems, and complicating efforts to meet international climate goals. In response to these challenges and the urgent need to mitigate climate change, countries have pledged to reduce methane emissions by 30 % from 2020 levels by 2030, aiming for a 0.2 °C reduction in global temperatures by 2050 (Cael and Goodwin, 2023). Similarly, the European Commission aims to achieve a climate-neutral Europe by 2050, as outlined in the European Climate Law (https://eur-lex.europa.eu/legal-content/EN/TXT/?uri=celex%3A32021R1119, last access: 30 October 2025), which mandates net-zero greenhouse gas (GHG) emissions (Rayner and Jordan, 2016). To effectively meet this goal, strategies must address reductions in both methane and carbon dioxide emissions, as both are crucial for mitigating global temperature rise. Methane is a significantly more potent greenhouse gas than carbon dioxide, with a global warming potential 28 times greater over a 100-year period and 84 times greater over 20 years (IPCC, 2023b; Myhre et al., 2013; Saunois et al., 2020). As a result, reducing methane emissions offers an especially effective short-term strategy for mitigating overall greenhouse gas emissions (Dlugokencky et al., 2011; Kikstra et al., 2022). Achieving this, however, requires accurate estimates of methane emissions. Despite substantial research, considerable uncertainty remains in identifying the geographic and temporal sources of these emissions.

Current national methane estimates reported to the United Nations Framework Convention on Climate Change (UNFCCC) predominantly rely on bottom-up methodologies, which apply emission factors to activity data, often supplemented by facility-specific information. However, these inventories are affected by significant uncertainties, often varying by a factor of two or more (Saunois et al., 2020; Solazzo et al., 2021), primarily due to the considerable variability in emission intensity across sources such as landfills, gas production facilities, and distribution networks (Leip et al., 2018). This variability cannot be fully captured by the use of generic emission factors, leading to substantial uncertainty in the resulting estimates. The reliance on uncertain and sparse input data, combined with poorly characterized emission factors, further undermines the accuracy of these estimates, especially when lacking comprehensive characterization. Moreover, since NGHGIs by design include only anthropogenic emissions, they do not account for natural methane sources and sinks. Therefore, complementary top-down approaches are essential to provide a more complete understanding of the total methane budget and to better constrain both anthropogenic and natural components.

To address these challenges, top-down approaches have increasingly been employed to generate independent, optimized emissions estimates. These methods use inverse techniques that assimilate observational data from in situ and/or satellite observations. Such estimates help to refine and constrain the data from bottom-up inventories. In recent years, the growing availability of greenhouse gas measurements and advancements in regional monitoring networks, particularly in Europe and North America, have significantly bolstered the effectiveness of top-down approaches (Bergamaschi et al., 2015). Numerous atmospheric inverse modelling studies have demonstrated the effectiveness of these approaches in quantifying methane emissions across regional to national scales. Such analyses have provided essential constraints on source magnitudes and spatial distributions, thereby informing climate mitigation strategies and policy development (e.g. Bocquet and Sakov, 2013; Saunois et al., 2020, 2025; Qu et al., 2021; Bergamaschi et al., 2022; Chen et al., 2022; Ernst et al., 2024; Petrescu et al., 2021, 2023; Steiner et al., 2024; Ioannidis et al., 2026). Complementing these scientific efforts, international initiatives such as the Global Carbon Project (GCP) (Friedlingstein et al., 2022), VERIFY (Petrescu et al., 2021), and EYE-CLIMA (https://eyeclima.eu, last access: 30 October 2025) have played a central role in developing coordinated frameworks for independent, observation-based assessments of greenhouse gas emissions and sinks. These programs integrate atmospheric observations, bottom-up inventories, and Earth system modelling to enhance transparency, support policy evaluation, and improve the robustness of emission estimates at regional to global scales. VERIFY, in particular, aimed to develop scientifically robust tools for verifying national emission inventories by integrating atmospheric observations, satellite data, emission inventories, and ecosystem models. Building on these efforts, the EYE-CLIMA project has taken a step further by focusing explicitly on reducing uncertainties in inversion-based estimates of methane and other greenhouse gases. Similarly, GCP provides comprehensive global assessments of carbon sources and sinks, supporting transparency and consistency in reporting and informing international climate agreements. These coordinated efforts have laid the groundwork for standardized, policy-relevant verification systems that bridge the gap between scientific research and national reporting. These studies highlight the essential role of top-down approaches in enhancing the reliability of national methane inventories through independent verification. For example, Chen et al. (2022) identified a 21 % upward correction needed for the Chinese methane inventory reported to the UNFCCC in 2019, demonstrating how top-down methods can uncover significant underestimations in bottom-up data. Similarly, Bergamaschi et al. (2022) reported elevated methane emissions for Germany, France, and the BENELUX countries in 2018 compared to those reported to the UNFCCC, further showcasing the ability of top-down approaches to reveal higher emissions than those recorded by bottom-up inventories. In contrast, Bergamaschi et al. (2022) also showed a close alignment between top-down estimates and both anthropogenic and natural bottom-up inventories for the UK and Ireland, illustrating how these methods can also validate bottom-up data. Although atmospheric inverse modeling is widely recognized as a valuable tool for verifying bottom-up estimates (IPCC, 2006), its incorporation into national reports faces several challenges. These include the limited availability of high-quality atmospheric measurements and uncertainties in transport models. Additionally, most past top-down studies have focused either on total emissions or on specific regions, often lacking the spatial resolution and source attribution needed to provide robust estimates at the national or sectoral level. As a result, while progress has been made in constraining total methane fluxes, there remains a clear research gap in applying inverse modeling to optimize emissions from individual sectors such as agriculture, waste, fossil fuels, wetlands, and geological sources. This gap limits our ability to disentangle source-specific contributions, thereby constraining the design of targeted and effective mitigation strategies. Furthermore, the complexity of atmospheric processes and the resource-intensive nature of inversion systems present additional barriers to integrating top-down estimates into national inventory frameworks.

In this study, we address these gaps by presenting high-resolution atmospheric inversion estimates for CH₄ emissions over Europe, covering the domain between 12° W, 37° E and 35, 73° N. Using the Community Inversion Framework (CIF; Berchet et al., 2021a), we apply a 4-dimensional variational optimization approach (4D-Var) driven by footprint estimates from FLEXPART v10.4 (Stohl et al., 1995; Pisso et al., 2019). The inversion assimilates data from 46 Integrated Carbon Observation System (ICOS) and non-ICOS in situ CH₄ observation sites across Europe (ICOS RI et al., 2023), providing monthly CH₄ emission estimates for the years 2017–2022 at a resolution of 0.2°×0.2°. In addition to optimizing total methane emissions, our framework performs sector-specific optimization by explicitly partitioning the control vector into major source categories, including the energy sector, agriculture, waste management, wetlands, and geological sources. While sectors are optimized separately, atmospheric CH₄ observations alone cannot fully distinguish emission sources without isotopic or co-emitted tracer information Mikaloff Fletcher et al., 2004; Turner et al., 2019. Nonetheless, the inversion yields a first-order, observation-constrained estimate of sectoral emissions that remains useful for interpreting regional variability and supporting national reporting. Due to computational expenses of the calculations, sensitivity and uncertainty analyses were conducted using a representative sample month, serving also as a verification step to identify discrepancies in the optimized fluxes arising from inversion setup choices. By integrating both total and sector-resolved inversion estimates, our study provides a more detailed spatial and source-specific characterization of methane emissions over Europe. These advances are crucial for informing climate policy, guiding sector-targeted mitigation strategies, and deepening scientific understanding of methane dynamics. The findings have the potential to substantially improve national reporting accuracy and support global efforts to reduce methane emissions in line with climate goals.

2 Datasets and Methodology

2.1 CIF-Flexpart Inversion Framework

We perform total and sector-specific methane inversions using the CIF coupled with a Lagrangian transport model FLEXPART and a 4D-Var optimization scheme, hereafter referred to as CIF-FLEXPART. CIF provides a unified platform for atmospheric inversions, supporting multiple transport models and enabling consistent assessments of greenhouse gas fluxes and their uncertainties. The posterior estimate is obtained by 4D-Var approach that seeks the optimal state vector x through minimizing the cost function:

$$\begin{array}{ll}{\displaystyle }J\left(\mathit{x}\right)=& {\displaystyle }{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}(\mathit{x}-{\mathit{x}}_b{)}^{\mathrm{T}}{\mathbf{B}}^{-\mathrm{1}}(\mathit{x}-{\mathit{x}}_b)\\ \text{(1)}& {\displaystyle }& {\displaystyle }+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}({\mathit{y}}_{\mathrm{o}}-\mathbf{H}(\mathit{x}\left){)}^{\mathrm{T}}{\mathbf{R}}^{-\mathrm{1}}\right({\mathit{y}}_{\mathrm{o}}-\mathbf{H}\left(\mathit{x}\right)),\end{array}$$

where x is the control vector of surface fluxes, x_b the prior estimate with covariance B, y_o the observed CH₄ mole fractions, and R the observation error covariance. Efficient minimization of the cost function necessitates the computation of its gradient with respect to the control vector. This gradient is evaluated using the adjoint operator H^⊤ of FLEXPART:

$$\begin{array}{}\text{(2)}& {\displaystyle }{\displaystyle }\mathrm{\nabla }J\left(\mathit{x}\right)={\mathbf{B}}^{-\mathrm{1}}(\mathit{x}-{\mathit{x}}_b)-{\mathbf{H}}^{\top }{\mathbf{R}}^{-\mathrm{1}}({\mathit{y}}_{\mathrm{o}}-\mathbf{H}(\mathit{x}\left)\right),\end{array}$$

which propagates mismatches between observed and modeled mole fractions back into the flux space. The minimization is performed iteratively using the M1QN3 quasi-Newton algorithm (Gilbert and Lemaréchal, 2009), a limited-memory variant of the variable metric method, which updates the state vector with each iteration. The iterations continue until the norm of the gradient falls below a preset convergence threshold, which in this study is set to 0.01 % of its initial value or a maximum of 30 simulations.

Unlike analytical Bayesian solutions, which provide a closed-form expression for the posterior error covariance, variational methods, such as 4D-Var, primarily yield the maximum a posteriori (MAP) flux estimate in high-dimensional settings. The analytical solution

$$\begin{array}{ll}{\displaystyle }{\mathit{x}}_a& {\displaystyle }={\mathit{x}}_b+\mathbf{K}({\mathit{y}}_{\mathrm{o}}-\mathbf{H}({\mathit{x}}_b\left)\right),\\ {\displaystyle }\mathbf{K}& {\displaystyle }={\mathbf{BH}}^{\top }({\mathbf{HBH}}^{\top }+\mathbf{R}{)}^{-\mathrm{1}},\\ \text{(3)}& {\displaystyle }\mathbf{A}& {\displaystyle }=({\mathbf{B}}^{-\mathrm{1}}+{\mathbf{H}}^{\top }{\mathbf{R}}^{-\mathrm{1}}\mathbf{H}{)}^{-\mathrm{1}},\end{array}$$

represents the exact Bayesian solution, where x_a is the posterior (optimized) flux state, K is the Kalman gain matrix, and A is the posterior error covariance matrix. In high-resolution applications, iterative numerical methods such as M1QN3 provide a computationally tractable approximation of the MAP estimate, but do not yield the full posterior covariance A in closed form. In this study, we adopt a variational 4D-Var framework implemented in CIF-FLEXPART to compute the maximum a posteriori (MAP) flux estimate. To quantify posterior uncertainties, we employ a Monte Carlo ensemble approach in which both prior fluxes and observations are perturbed, and the posterior variance is derived from the ensemble spread. Specifically, posterior uncertainties are estimated by approximating the posterior error covariance (P_a,ii) as the sample covariance of the optimized flux ensemble (Bocquet and Sakov, 2013):

$$\begin{array}{}\text{(4)}& {\displaystyle }{\displaystyle }{P}_{a,ii}={\displaystyle \frac{\mathrm{1}}{N-\mathrm{1}}}\sum _{k=\mathrm{1}}^N{\left(x_{a,i}^{\left(k\right)}-{\overline{x}}_{a,i}\right)}^{\mathrm{2}},\end{array}$$

where $x_{a,i}^{\left(k\right)}$ denotes the optimized state for ensemble member k, and ${\overline{x}}_{a,i}$ is the ensemble mean. The posterior uncertainty is then calculated as ${\mathit{\sigma }}_{\text{post}}=\sqrt{P_{a,ii}}$. Thus, although CIF-FLEXPART does not directly provide the full posterior covariance matrix, posterior uncertainties can be effectively approximated using ensemble-based methods.

The prior covariance matrix B is modeled as:

$$\begin{array}{}\text{(5)}& {\displaystyle }{\displaystyle }\mathbf{B}={\mathbf{D}}^{\mathrm{T}}\cdot \mathbf{C}\cdot \mathbf{D}\end{array}$$

where D is a diagonal matrix representing relative flux uncertainties, set to 50 % of the prior fluxes for all optimized categories. The state vector is defined in physical flux units, with each element corresponding directly to the flux of a given category or grid cell. Accordingly, the variances in D are expressed relative to the prior flux magnitudes, representing proportional uncertainties in the physical flux values rather than dimensionless scaling factors. The matrix C contains the correlation structure, with off-diagonal elements modeled using a Gaussian function that decays exponentially with spatial distance r as $\exp (-r^{\mathrm{2}}/l^{\mathrm{2}})$ (Gaspari and Cohn, 1999; Peters et al., 2005). The spatial correlation length is set to l_land=200 km over land and l_ocean=500 km over oceans, while the temporal correlation length is prescribed as 90 d.

Observational errors are represented by R, which accounts for instrument precision, model representativeness, and transport model errors. A minimum uncertainty of 5 ppb is imposed on all station observations to account for ∼3 ppb of transport model error. To avoid biasing surface flux estimates, we jointly optimize the total CH₄ fluxes (f), which include contributions from both the nested domain and the global background fluxes, together with the initial atmospheric mole fractions (c₀), representing the model's initial atmospheric state. Thus, the background contribution to the observations is implicitly represented through the combination of optimized initial mole fractions and flux contributions from outside the domain. Initial condition uncertainties of 0.5 % (∼10 ppb) were applied, consistent with previous studies that used values in the range of 0.05 %–1 % (Thompson and Stohl, 2014; Szénási et al., 2021; Steiner et al., 2024; Ioannidis et al., 2026). We also tested a lower value of 0.05 % in sensitivity experiments.

In our setup, we optimize the initial mole fractions together with the fluxes, and the extended control vector is:

$$\begin{array}{}\text{(6)}& {\displaystyle }{\displaystyle }\mathit{x}=\left[\begin{array}{c}\mathit{f}\\ {\mathit{c}}_{\mathrm{0}}\end{array}\right],\end{array}$$

where f represents the total CH₄ fluxes from both the nested domain and the global background to be optimized, and c₀ represents the initial mole fraction. Optimizing c₀ ensures the transport model starts from a consistent atmospheric state, which is crucial for species like CH₄ where early-time conditions affect the entire simulation. The correlation structure of c₀ is defined consistently with the fluxes, using a Gaussian function that decays with spatial distance r as $\exp (-r^{\mathrm{2}}/l^{\mathrm{2}})$ (Gaspari and Cohn, 1999; Peters et al., 2005). The spatial correlation length is set to l_land=200 km over land and l_ocean=500 km over oceans, while the temporal correlation length is prescribed as 30 d. This formulation allows CIF to adjust both the emissions and the initial atmospheric state to best match the observed mole fractions, while respecting the prior uncertainties and correlation structures.

Sector-specific inversion, conducted separately from the total-flux inversion, all emission sectors are explicitly included in the state vector together with the initial mole fractions, and optimized simultaneously using the same concentration observations. In this case, the control vector in Eq. (6) becomes:

$$\begin{array}{}\text{(7)}& {\displaystyle }{\displaystyle }\mathit{x}=\left[\begin{array}{c}{\mathit{f}}_{\text{Wetlands}}\\ {\mathit{f}}_{\text{Agriculture}}\\ {\mathit{f}}_{\text{Energy}}\\ {\mathit{f}}_{\text{Waste}}\\ {\mathit{f}}_{\text{Industrial}}\\ {\mathit{f}}_{\text{Fire}}\\ {\mathit{f}}_{\text{Geological}}\\ {\mathit{f}}_{\text{Termites}}\\ {\mathit{f}}_{\text{Ocean}}\\ {\mathit{c}}_{\mathrm{0}}\end{array}\right]\end{array}$$

Here, each sub-vector f_s represents emissions from sector s. Unlike the total-flux inversion, where only the aggregated flux is optimized, the sector-specific inversion optimizes each sector independently using the same concentration observations. The prior covariance B encodes both uncertainties and cross-sector correlations, thereby guiding how flux adjustments are redistributed among sectors. By preserving the correlation structure in B, the inversion allows sector-specific adjustments. Although this approach does not replace source-specific tracers (e.g. isotopes or co-emitted species), it provides a first-order quantification of sectoral contributions to methane variability. The prior covariance B for the sector-specific inversion is constructed with diagonal elements set to 50 % of the sectoral prior flux. A corresponding correlation matrix C, consistent with the total flux inversion, is also assumed. The matrix C is defined using isotropic Gaussian spatial correlations, $\exp (-r^{\mathrm{2}}/l^{\mathrm{2}})$, with l_land=200 km and l_ocean=500 km, together with a temporal correlation length of 90 d. These assumptions are applied uniformly across all sectors, without sector-specific covariance tuning.

2.2 FLEXPART in the Inversion Framework

We employ the Lagrangian Particle Dispersion Model FLEXPART v10.4 (Pisso et al., 2019), a widely used open-source transport model for simulating the dispersion and turbulent mixing of atmospheric tracers (Stohl et al., 1995). Within our inversion framework, FLEXPART is used to quantify source–receptor relationships (SRRs) between surface methane fluxes and atmospheric observations.

Meteorological input fields were taken from the ERA5 reanalysis of the European Centre for Medium-Range Weather Forecasts (ECMWF). We use hourly data on 137 vertical levels, regridded to a horizontal resolution of 1°×1°. These data are pre-processed for FLEXPART using the FlexExtract toolbox (Tipka et al., 2020). FLEXPART is operated in backward mode: for each observation, 10 000 pseudo-particles are released at the receptor location and traced backward in time for 10 d. The resulting surface flux footprints are archived at an hourly temporal resolution with a spatial resolution of 0.2°×0.2° for the nested domain and 2°×2° for the global domain. The dominant atmospheric sink of CH₄ through oxidation by OH is explicitly represented in FLEXPART as a first-order chemical loss process along particle trajectories (Pisso et al., 2019). This loss is computed using temperature-dependent reaction rate coefficients and prescribed three-dimensional OH concentration fields from the GEOS-Chem model. The OH fields are provided as monthly mean distributions with global coverage and vertical structure, consistent with standard GEOS-Chem simulations of tropospheric oxidant chemistry. The particle residence time within each grid cell is proportional to the observation's sensitivity to fluxes in that grid cell, yielding the source–receptor sensitivity matrix

$$\begin{array}{}\text{(8)}& {\displaystyle }{\displaystyle }{\mathbf{S}}_{ij}={\displaystyle \frac{\partial y_i}{\partial {\mathit{q}}_j}},\end{array}$$

where S_ij represents the sensitivity of observation y_i to a surface flux q_j in grid cell j. These SRRs explicitly account for advection, turbulence, convection, and deposition processes (Stohl et al., 1995; Seibert and Frank, 2004; Pisso et al., 2019).

For long-lived greenhouse gases such as CH₄, it is essential to represent not only local and regional flux contributions but also the large-scale background concentration (Thompson and Stohl, 2014). CIF accounts for this by employing two types of sensitivities within the observation operator: (i) surface-flux sensitivities, separated into nested-domain and outside-nested contributions, and (ii) three-dimensional (3D) concentration sensitivities that are multiplied by external 3D mole-fraction fields to represent the large-scale background. Sensitivities are computed by FLEXPART and include methane loss due to oxidation by OH. The modeled methane mixing ratio at receptor i (${\mathit{y}}_i^m$) is then given by:

$$\begin{array}{}\text{(9)}& {\displaystyle }{\displaystyle }{\mathit{y}}_i^m=\sum _j{\mathbf{S}}_{ij}^{\text{nest}}{\mathit{f}}_j^{\text{nest}}+\sum _j{\mathbf{S}}_{ij}^{\text{out}}{\mathit{f}}_j^{\text{out}}+\sum _k{\mathbf{S}}_{ik}^{\mathrm{3}\mathrm{D}}{\mathit{c}}_k,\end{array}$$

where ${\mathbf{S}}_{ij}^{\text{nest}}$ and ${\mathbf{S}}_{ij}^{\text{out}}$ denote the FLEXPART sensitivities of observation i to surface fluxes inside and outside the nested domain, respectively, and ${\mathbf{S}}_{ik}^{\mathrm{3}\mathrm{D}}$ denotes the sensitivity to the CH₄ mole fraction in 3D grid cell k. For the sector-specific inversion flux vectors f_j can be further decomposed into sectoral contributions f_js Eq. (7) (e.g. agriculture, wetlands, biomass burning), such that ${\mathit{y}}_i^m$ becomes:

$$\begin{array}{ll}{\displaystyle }{\mathit{y}}_i^m=& {\displaystyle }\sum _s\sum _j{\mathbf{S}}_{ij}^{\text{nest}}{\mathit{f}}_{js}^{\text{nest}}+\sum _s\sum _j{\mathbf{S}}_{ij}^{\text{out}}{\mathit{f}}_{js}^{\text{out}}\\ \text{(10)}& {\displaystyle }& {\displaystyle }+\sum _k{\mathbf{S}}_{ik}^{\mathrm{3}\mathrm{D}}{\mathit{c}}_k.\end{array}$$

Thus, each modeled observation combines (i) local and global flux enhancements, optionally partitioned by sector, and (ii) large-scale background contributions from external 3D mole-fraction fields.

2.3 In situ observations of CH₄ mixing ratio

We utilise surface methane concentration measurements primarily from the Integrated Carbon Observation System (ICOS), which provides a harmonized European network of atmospheric CH₄ observations widely used in inverse modelling. The dataset includes both ICOS observations, which follow standardized ICOS protocols for calibration, quality control, and data processing, and non-ICOS observations, which are collected outside the ICOS quality-controlled may follow different quality-assurance procedures. In the base inversion setup, we assimilate data from 44 in situ stations obtained from the obspack_CH₄_466_GVeu_v10.0_20240729 European CH₄ time series (ICOS RI et al., 2023). Additionally, we included data from two Finnish Meteorological Institute (FMI) stations, Kumpula and Sodankylä, which have been identified as reliable observation sites (Tsuruta et al., 2019; Tenkanen et al., 2025). To ensure robust observational constraints on CH₄ emissions, we applied the site selection criteria separately for each year, retaining only stations with at least 30 d of data coverage in each year. A preliminary screening was performed to exclude stations for which the model did not adequately reproduce observed variability or magnitude. This screening removed only one site, Ispra (IPR, Italy), which exhibited consistently low correlation and comparatively high bias relative to observations, consistent with previous studies (Steiner et al., 2024). In addition, stations TAC and ZSF were excluded because they fell within the same model grid boxes as WAO and HPB, respectively. Since the two stations in each grid box likely experience identical modelled conditions, only one was retained, with preference given to ICOS and low-altitude stations. Furthermore, several mountain stations in the Alps were removed due to the model's limited ability to represent complex topography at its resolution.

For the sensitivity analysis, an additional set of 10 stations was included, comprising: (i) the station with poor model–observation agreement (IPR), (ii) stations located in grid boxes shared with another site, and (iii) the previously omitted Alpine mountain stations. These additional stations were included in the sensitivity analysis to assess the impact of observation network density.

Figure 1a shows the locations of all observation sites and the nested modelling domain. The study domain spans 12° W, 37° E and 35, 73° N and is hereafter referred to as Domain/Europe. It is further subdivided into subregions: the 27 European Union member states (EU27) plus the UK, Norway, and Switzerland (EU27+3); and four subregional groupings (Northern, Central, Western, and Southern). A detailed list of countries within each subregion is provided in Appendix A Table A1. Figure 1b displays the daily-averaged concentration time series for each station, calculated over the assimilated hours only (see Appendix A Table A2 for a complete list of stations). For a detailed description of the measurement procedures and data processing, see Ramonet et al. (2020), Hazan et al. (2016).

Figure 1 Distribution of the CH₄ observational network used in the inversion. Black markers indicate low-altitude stations (<1000 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$), while blue markers represent high-altitude stations (>1000 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$). Stations included in the sensitivity analysis are highlighted with star symbols. The EU27+3 countries are grouped into four sub-regions, as indicated by the color-coded legend. (b) Overview of daily mean CH₄ mole fractions from each assimilated station, including the timing of ICOS labelling. Data points after the red marker are ICOS-labelled, indicating that measurements have been processed according to ICOS protocols from the indicated date onward. White gaps indicate periods with missing data. (c) Relative variability reduction (%) for each station, defined as $(\mathrm{1}-\frac{{\mathit{\sigma }}_{\text{window}}}{{\mathit{\sigma }}_{\text{full}}})\cdot \mathrm{100}$, where σ_full is the standard deviation over the full period (2017–2022), and σ_window is computed over 14:00–16:00 LT for low-altitude stations and 02:00–04:00 LT for high-altitude stations.

In cases where multiple intake heights were available, such as at the Cabauw station, where intakes were positioned at 27, 67, 127, and 207 m a.g.l., we chose to use data exclusively from the highest intake height. This approach was adopted to mitigate the difficulties that transport models encounter in accurately representing concentration gradients close to the ground, where such gradients tend to be steep and variable. By selecting the highest intake height, we ensured that the assimilated data more closely represent the characteristics of a well-mixed boundary layer, thereby improving the consistency of the modeled methane mole fractions (Vermeulen et al., 2011). Vertical gradients of CH₄ near the surface are notoriously difficult to simulate accurately due to complex local meteorology, surface interactions, and diurnal variability (Peltola et al., 2015).

To further reduce these challenges, we use recommended quality-controlled observations during periods when vertical gradients are expected to be weaker, thereby minimizing potential transport model errors. For low-altitude sites (≤1000 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$) over relatively flat terrain, we assimilate 3 h afternoon averages (14:00–16:00 LT), when the boundary layer is typically well mixed and vertical mixing is strong. For high-altitude or mountainous sites (>1000 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$), we instead use 3 h nighttime averages (02:00–04:00 LT), which are more representative of free-tropospheric conditions and less affected by daytime valley winds that are difficult to resolve at coarse model resolution. By focusing on these time windows, we aim to reduce biases associated with shallow nocturnal boundary layers and unresolved slope and valley circulations, thereby improving the reliability of the assimilation.

This strategy is consistent with common practices in atmospheric inverse modeling, where daytime observations are generally preferred for lowland sites to avoid stable nocturnal conditions, while nighttime observations are often used in mountainous regions to limit the influence of unresolved diurnal circulations (Steiner et al., 2024; Monteiro et al., 2024). To further characterize temporal variability, we computed the standard deviation of CH₄ mole fractions within these periods and quantified the relative reduction compared to the full daily variability (Fig. 1c). The mean standard deviation decreases from 39.4 ppb (range: 18.9–125.7 ppb) to 37.2 ppb (range: 18.8–103.8 ppb), corresponding to a mean reduction of 4.0 % (range: −4.2 % to 23.6 %). This indicates that the selected windows generally capture more stable atmospheric conditions, supporting their use in the inversion, despite some site-specific deviations.

2.4 Prior fluxes

High-resolution monthly data are employed for key contributors (anthropogenic, biospheric and fire) to capture fine-scale spatial variability and to identify localized hotspots. For sources with sparse data availability and limited temporal variability (ocean, geological and termites), coarser climatological estimates are used. However, all datasets are regridded to match the inversion spatial resolution (0.2°×0.2° for the nested European domain and 2°×2° for the global background). The inversion is performed at monthly temporal resolution, optimizing fluxes independently for each month. Within the nested European domain, agriculture is the dominant source, contributing approximately 11.7 Tg yr⁻¹ (estimates for 2020, 34 % of the annual regional total). Other major contributors include waste (9.5 Tg yr⁻¹, 28 %), energy (5.5 Tg yr⁻¹, 16 %), geological emissions (5.1 Tg yr⁻¹, 15 %) and wetlands (2.9 Tg yr⁻¹, 9 %), smaller contributions arise from ocean (0.7 Tg yr⁻¹, 2 %), termites (0.6 Tg yr⁻¹, 2 %), industrial processes (0.03 Tg yr⁻¹, <1 %) and fires (0.03 Tg yr⁻¹, <1 %), with all estimates in brackets corresponding to the year 2020. These individual source contributions are combined to derive the total prior CH₄ emissions used in the total-flux inversion, while the sector-specific inversion optimizes each sectoral flux component separately within the expanded control vector. A comprehensive overview of the prior data is provided in Table 1, with the corresponding map of mean fluxes for 2017–2022 shown in Appendix A Fig. A1.

Table 1 Overview of CH₄ emission sources, their native resolutions, and average total emissions for the global and nested domains (Tg yr⁻¹). The reported values correspond to 2020 prior emissions.

Monthly anthropogenic CH₄ emissions are derived from the Emission Database for greenhouse gas and Atmospheric Research (EDGARv8.0) (Crippa et al., 2020) and Greenhouse Gas and Air Pollution Interactions and Synergies (GAINS) (Höglund-Isaksson et al., 2020) inventories, both providing high-resolution (0.1°×0.1°) global estimates of greenhouse gas emissions across multiple sectors. For EU27+3, we use emission estimates from the GAINS inventory, which incorporates country-specific activity data and mitigation assumptions, providing regionally consistent estimates that reflect national circumstances more closely than the globally uniform EDGARv8.0 inventory. On average, GAINS estimates are 3 Tg yr⁻¹ (ranging from 1 to 4 Tg yr⁻¹) lower than EDGARv8.0, amounting to approximately 10 % of the total anthropogenic CH₄ emissions in the EU27+3 region during the study period. For categorising anthropogenic emission sources, we adopted the IPCC (2006) Common Reporting Format (CRF) to classify the emissions into four main sectors: Energy, Industrial Processes and Product Use, Agriculture, and Waste. This sectoral breakdown is consistent with UNFCCC NGHGI reporting guidelines and used by major inventories such as GAINS and EDGAR, allowing for direct comparison with national reports.

Biospheric emissions in our prior estimates are sourced from the JSBACH model, which simulates key ecosystem processes such as photosynthesis, respiration, carbon and water cycling, vegetation dynamics, and land ecosystem processes (Reick et al., 2021). Specifically, we use the version of JSBACH coupled with HIMMELI (Raivonen et al., 2017) (JSBACH-HIMMELI), which provides hourly methane emissions at high spatial resolution (0.1°×0.1°) over Europe. The modeled CH₄ fluxes include contributions from peatlands, water-saturated and inundated soils, and net mineral soils, thereby representing the dominant diffuse terrestrial and freshwater-influenced methane sources. As the JSBACH-HIMMELI product used in this study provides data only for the European domain, we complement it with global biospheric emissions from the LPX-Bern DYPTOP v1.4 dataset (Lienert and Joos, 2018) to represent the global background.

Monthly prior fire emissions are obtained from the Global Fire Emissions Database version 4 (GFED4s) described in van der Werf et al. (2017) at a spatial resolution of 0.25°×0.25°, excluding methane emissions from biomass and agricultural waste burning, which were already included in GAINS and EDGAR. Oceanic CH₄ emissions are represented by climatological estimates from Weber et al. (2019), while geological emissions are based on the dataset by Etiope et al. (2019), globally scaled to 23 Tg yr⁻¹, following the Intergovernmental Panel on Climate Change (IPCC) AR6 Working Group I report (IPCC, 2023a; Tsuruta et al., 2023). Termite-related emissions are included according to the estimates provided by Saunois et al. (2020). The total flux constructed using the GAINS inventory for the EU27+3, combined with the other data sources listed in Table 1, is hereafter referred to as the “GAINS-based” prior.

2.5 Sensitivity Tests

We conducted eight sensitivity experiments (S1–S8) for the total-flux inversion to evaluate how key elements of the inversion setup influence the posterior CH₄ flux estimates and the associated error reduction. The base inversion (S1) serves as a reference for comparison. In S2, the initial CH₄ mixing ratio fields from CAMSv22r2 (Bergamaschi et al., 2013) are replaced with those from CTE-CH₄ (Tenkanen et al., 2025), both of which are outputs of global atmospheric inversion systems based on the TM5 transport model and available at daily resolution. In this study, we use these datasets for the 2017–2022 period. These fields are used to compute sensitivities to the mixing ratio at the endpoints of FLEXPART back-trajectories. This substitution allows us to assess the sensitivity of the inversion to the choice of initial mixing ratio and associated large-scale background representation. Although both products rely on TM5-based transport (Krol et al., 2005), they differ in their inversion configurations, including data assimilation methods, flux resolution, prior constraints, and assimilated observations. Inversions S3 and S4 address uncertainty assumptions. In S3, the prior flux uncertainty was increased from 50 % to 100 %, allowing more flexibility for the inversion to adjust emissions. In S4, the uncertainty assigned to initial mixing ratios was reduced from 0.5 % to 0.05 %, thereby increasing the influence of observed enhancements on the posterior flux adjustments. In S5, the observational network was expanded from 46 to 56 sites, including a denser station distribution along the Swiss Alps (see Fig. 1), to examine the sensitivity to observational coverage. In S6, the prior emissions were replaced from GAINS-based estimates to EDGARv8, which are approximately 10 % higher, allowing assessment of how differences in emission inventories affect the posterior flux estimates. Finally, S7 and S8 address spatial error correlations in the prior. The horizontal correlation length was reduced from 200 to 50 km (S7) and 20 km (S8), enabling finer spatial adjustments in the inversion. All inversions were performed using a 4D-Var ensemble approach, in which 30 Monte Carlo realizations were generated by perturbing both prior fluxes and observational data. Posterior uncertainties were estimated from the ensemble by approximating the posterior error covariance discussed in Sect. 2.1, Eq. (4). The error reduction was computed as $\mathrm{1}-\frac{{\mathit{\sigma }}_{\text{post}}}{{\mathit{\sigma }}_{\text{prior}}}$, where σ_post and σ_prior denote the posterior and prior uncertainties, respectively. Due to computational constraints, the sensitivity experiments were carried out for July 2021 only. This summer-month analysis is intended as a diagnostic evaluation of inversion-framework sensitivity and does not capture the full seasonal variability. A summary of the inversion configurations is provided in Table 2. The base inversion (S1) configuration is applied consistently throughout the study period, and the results presented in the total-flux analysis section are based on this setup unless otherwise stated. The sectoral inversion shares the same general configuration as the total-flux inversion but employs an extended control vector that explicitly decomposes emissions into individual source sectors, enabling simultaneous optimization of sector-specific fluxes within CIF-FLEXPART.

Table 2 Inversion setups for the sensitivity analysis. Deviations from the base inversion (S1) are indicated with an asterisk (*).

3 Results and discussion

3.1 Comparison of Simulated and Measured Methane Mole Fractions

Standard metrics such as mean bias, root-mean-square error (RMSE), and correlation play a crucial role in evaluating the performance of atmospheric inversion models by enabling direct comparison between simulated and observed CH₄ mole fractions. In this section, these metrics are applied to assess the CIF-FLEXPART inversion simulated mole fractions against the assimilated in situ observations. Figure 2 displays time series of both prior and posterior mole fractions from the base inversion set-up, alongside observations from selected European in situ stations in 2021, which cover a range of geographical settings and emission regimes. The prior methane mixing ratios generally underestimate the observations, resulting in a mean negative bias at most sites, whereas the posterior estimates substantially reduce this bias.

Figure 2 Time series of observed and modeled CH₄ mole fractions at selected stations, averaged at the assimilated hour for 2021: BIS (Biscarrosse, France), CBW (Cabauw, Netherlands), CRP (Carnsore Point, Ireland), ERS (Ersa, coastal Mediterranean site), HPB (Hohenpeissenberg, mountain station in Germany), KIT (Karlsruhe, Germany), PAL (Pallas, Finland), and TRN (Trainou, France), where the numbers in parentheses indicate the station's latitude and longitude. Shown are observations (black), prior simulations (orange), posterior simulations (green), and background contributions (cyan for prior, purple for posterior). Text boxes in each panel report mean bias (MBias) in ppb, root mean square error (RMSE) in ppb, and correlation coefficient (r²) for the prior and posterior relative to the observations.

These improvements are clearly visible in the time series shown in Fig. 2, illustrating the discrepancies between observed and prior-simulated mole fractions and their correction in the posterior. For instance, at stations BIS, CRP, HPB, KIT, and TRN, the prior exhibits a pronounced negative bias (−11.1 to −22.4 ppb) and systematically underestimates sharp observational peaks when regional fluxes are elevated. The prior initial concentration contribution is overly smoothed and lacks the sharp variability of the observed and simulated mole fractions, though it still shows meaningful variations. In contrast, the posterior reduces these biases (mean bias −3.2–3.8 ppb) and improves the representation of both the initial condition and regional flux contributions, providing a more realistic depiction of the regional signal. This adjustment is further reflected in the corresponding flux contributions in Appendix A Fig. A2, where regional fluxes are increased in the posterior to better match the observations. At CBW, where local influences dominate, the prior fails to reproduce the full variability (see Fig. 2, CBW); the posterior reduces the overall bias, although some high-frequency discrepancies persist. The regional flux contributions show substantial adjustments in the posterior, consistent with the improved match to observations. At PAL, prior mole fractions generally match observations, but pronounced mismatches occur during the first months of 2021. These are largely corrected by adjustments to the initial concentration contribution, while the global flux contribution remains relatively more influential and unchanged in the posterior (see Appendix A Fig. A2). At ERS, unlike most other stations, the prior exhibits a mean positive bias of 8.5 ppb, which is reduced to 4.1 ppb in the posterior. The background contribution is adjusted upward, reinforcing the positive bias, while the regional flux contribution is lowered in the posterior. Together, these adjustments improve the overall agreement between simulated and observed mole fractions, illustrating the inversion's ability to correct both the initial condition and the regional flux components. Overall, these results demonstrate that the inversion effectively corrects biases in both the initial conditions and the regional fluxes, improving the representation of CH₄ dynamics across stations.

These station-level mole fractions are further illustrated in Fig. 3, which evaluates model–observation agreement across all stations and years (2017–2022). The coefficient of determination (R²) increased from 0.54 to 0.76, RMSE decreased from 33.4 to 24.4 ppb, and mean bias was reduced from −12.2 to −2.1 ppb. The posterior regression slope (0.75) is closer to unity than the prior (0.57), and Pearson correlation increased from 0.78 to 0.87 (p<0.001), reflecting enhanced responsiveness to observed changes. Modelled vs. observation performance metrics for all assimilated stations are summarized in Appendix A Table A3. Our posterior statistical results closely align with the findings of Steiner et al. (2024), Bergamaschi et al. (2022). By assimilating data from 28 stations, Steiner et al. (2024) reported a Pearson correlation of approximately 0.7 for 25 of the stations. Figure 3 also reveals that the prior simulation systematically underestimates high CH₄ mole fractions, with many data points falling below the regression line. This underestimation is particularly pronounced at stations such as LUT, LMU, STE, and CBW, all of which exhibit strong negative biases (below −35 ppb). Although the posterior simulation substantially reduces this underestimation, with biases improving to around −10 ppb, some discrepancies persist at the highest mole fractions, likely reflecting unresolved local emissions or sub-grid variability. These results highlight the importance of targeted improvements in emission inventories and model resolution for stations with persistent negative bias.

Figure 3 Modeled vs. observed CH₄ mole fractions for all assimilated in situ observations. Left panels: posterior vs. observations; right panels: prior vs. observations. Statistical metrics are shown in each panel. Black dashed line: one-to-one; red line: best-fit regression.

To further assess the level of agreement across subregions, we compare posterior and prior simulations for different categories. This includes geographical regions (Central, Northern, Western, Southern) consistent with Petrescu et al. (2021), as well as elevation (mountain vs. non-mountain) and station type (coastal vs. non-coastal). Figure 4 presents the distributions of residuals using violin plots, with the 25th–75th interquartile range indicated by a vertical line. Posterior residuals exhibit systematically narrower spreads and improved alignment with observations. Across regions, RMSE and interquartile range (IQ75) decrease by 26 %–37 %, with median residuals moving closer to zero, except in the Southern region. Mountain stations show a smaller spread reduction (IQ75: 21.5 vs. 27.6, RMSE reduction 26 %), reflecting higher variability likely due to complex transport, whereas non-mountain, coastal, and non-coastal stations show 28 %–32 % RMSE reductions. Overall, these analyses demonstrate that the inversion effectively reduces uncertainty and bias across all categories, enhancing model fidelity and observational agreement (see Table 3).

Figure 4 Residuals between observed (obs) and modelled (sim) CH₄ concentrations for assimilated in situ stations, grouped by subregional category (see Appendix A TableA2) The horizontal dashed line indicates zero residual, and the vertical bar shows the 25th–75th interquartile range. Negative values indicate that simulated concentrations are lower than observed.

Table 3 Prior and posterior performance metrics of CH₄ mole fraction simulations from the total flux inversion setup grouped by sub-regional category (see in Appendix A Table A2). Reported statistics include regional mean bias, the interquartile ranges: the 25th–75th (IQ75), and the 5th–95th range (IQ95), root mean square error (RMSE), relative RMSE reduction (ΔRMSE), number of observation stations, and the total number of data points (N-data) used in the evaluation.

3.2 Posterior Flux

In Fig. 5, we presented the six-year mean methane fluxes (2017–2022) over the study domain. Figure 5a–c display prior estimates, posterior estimates, and their differences, while Fig. 5d–e present monthly time series and inland average estimates for subregions. The spatial maps highlight regions of substantial flux adjustments, and the time series illustrate temporal variability and regional contrasts between prior and posterior estimates. The inversion reveals pronounced regional corrections, with posterior emissions being systematically higher over BENELUX, France, and Germany, pointing to underestimation in the prior, whereas reductions occur over the UK, Italy, and Romania, indicating prior overestimation. These patterns align well with earlier inversion studies (Bergamaschi et al., 2022; Steiner et al., 2024). In regions with sparse observational data, such as the Iberian Peninsula and eastern Europe, the prior fluxes change only slightly.

Figure 5 Top row: Mean methane (CH₄) fluxes over the EU27+3 domain for the study period (a) prior fluxes, (b) posterior fluxes after inversion, and (c) flux increments (posterior – prior). Middle rows: Monthly CH₄ flux time series averaged over the entire domain, within EU27+3, and outside EU27+3. Bottom row: Monthly CH₄ flux time series for selected EU27+3 subregions (Northern, Central, Southern, and Western), showing prior (dashed lines) and posterior (solid lines) estimates. Bars on the right summarize mean annual prior and posterior fluxes estimates. Shaded areas denote the corresponding prior and posterior uncertainties in the time series, while error bars on the bars represent the associated annual uncertainties.

Despite these pronounced regional adjustments, the mean emission exhibits only a modest net change at the domain scale, with average total methane emissions shifted 36.1±4.5 Tg yr⁻¹ in the prior to 36.7±4.0 Tg yr⁻¹ in the posterior, a rise of 1.7 % (see Fig. 5d), reflecting strong subregional contrasts. Notable posterior CH₄ flux increases are observed in Central Europe (+32 %) and Western Europe (+19 %), in contrast, Southern Europe shows a decrease of 17 % (see Fig. 5e). The monthly time series reveals a clear seasonality in northern regions, with a peak in summer months. The southern region exhibits a slight opposing seasonality, while other regions show no clear seasonal cycle. Additionally, it displays that posterior emissions from Central and Western Europe are consistently higher than the prior throughout the study period, while in the southern region, the posterior emissions are generally lower than the prior. At the domain and EU27+3 scales, however, there is no clear or persistent difference between posterior and prior emissions, and no significant trend is evident over the 2017–2022 period at either the domain or subregional scale. Over the EU27+3, the average posterior flux for 2017–2022 is 23.3±2.3 Tg yr⁻¹, representing a 6.6 % increase relative to the prior (21.7±2.7 Tg yr⁻¹). Our results falls within the 22–26 Tg yr⁻¹ range reported by Petrescu et al. (2021), Petrescu et al. (2023) for 2006–2017 based on three inversion setups. Although our estimate is at the lower end of this range, it provides a meaningful comparison given the difference in time periods, as EU27+UK emissions show no strong trend over the interval considered.

While spatial patterns in total methane flux adjustments indicate where the inversion has modified prior estimates, analysing these changes by sector reveals the dominant sources driving regional flux adjustments and highlights systematic biases in prior inventories. Figure 6 presents the spatial distribution of sector-specific methane flux increments (posterior – prior) over Europe for 2017–2022. Panels (a)–(d) correspond to anthropogenic sectors: agriculture (AGR), energy (ENG), waste (WST), and industry (IND); whereas panels (e)–(g) depict natural sources, including wetlands (WET), geological seepage (GEO), and biomass burning (FIR). Posterior agricultural emissions increase markedly across the BENELUX region, suggesting systematic underestimation in the prior inventory. Adjustments in other anthropogenic sectors are generally more spatially localized. Energy-related emissions decrease over the North Sea, but localized increases are observed in Belgium, Luxembourg, Ukraine, and western Russia. Waste emissions show localized hotspots in BENELUX, France, and Poland, while decreasing in Italy and the UK. Industrial emissions rise predominantly over BENELUX and Germany. CH₄ emissions from wetland are reduced, particularly in northern Europe and the UK, largely reflecting a reduction in summer fluxes. Reductions in the GEO sector over Italy and Romania are consistent with previously reported overestimations arising from the global scaling of geological emission factors. Steiner et al. (2024) likewise identified strong posterior reductions over Italy and Romania, which they attributed to inflated geological emissions in these regions. Although prior geological emissions were harmonized to 23 Tg yr⁻¹ at the global scale, the inversion results highlight the need for regionally differentiated emission estimates. Footprint sensitivities from FLEXPART (Appendix, Figs. A3–A4) reveal a heterogeneous observational constraint, with moderate sensitivities in Italy and Romania. Nevertheless, substantial flux corrections are inferred due to the combination of moderate footprint sensitivities and relatively large prior geological emissions, yielding sufficient emission sensitivity to drive posterior adjustments.

Figure 6 Spatial distribution of sector-specific methane flux increments (posterior – prior) over Europe averaged in the years (2017–2022). Panels show adjustments for agriculture (AGR), energy (ENG), waste (WST), industry (IND), wetlands (WET), geological sources (GEO) and fires (FIR). Agricultural emissions exhibit the largest posterior increase, particularly across central and western Europe, whereas wetlands show strong reductions in northern Europe. Geological emissions decrease mainly in southern Europe and Romania. Note that the colour scale is different.

Figure 7 summarizes the prior and posterior fluxes by European subregion. Agricultural emissions show the largest positive adjustments. In central EU27+3, emissions increase from a prior estimate of 2.4 to a posterior estimate of 3.2 Tg yr⁻¹ (+33 %), and in western EU27+3 they rise from 3.5 (prior) to 4.5 Tg yr⁻¹ (posterior), which is a 29 % rise. The strongest relative increase occurs in the BENELUX region, where mean agricultural emissions nearly double from 0.7 to 1.3 Tg yr⁻¹ (+86 %). At the domain scale, posterior agricultural CH₄ emissions amount to 13.6 Tg yr⁻¹, representing a 16 % increase compared with the prior. EU27+3 posterior emissions are estimated at 10.8 Tg yr⁻¹, corresponding to a 21 % increase relative to the GAINS prior (Fig. 7a). Agricultural emissions account for 46 % of the total EU27+3 CH₄ budget in our estimate, representing the largest sectoral share. This aligns with consolidated assessments and UNFCCC NGHGI reports, which attribute 52.4 % (±8.7 %) of EU27+UK CH₄ emissions in 2019 to agriculture (Petrescu et al., 2021, 2023). The dominant role of agriculture in total anthropogenic CH₄ emissions also holds globally (IPCC, 2019). Emissions from energy-related sectors show modest increases in central Europe, rising from 0.8 to 0.9 Tg yr⁻¹ (+12 %), and in western Europe, the prior estimate of 0.3 Tg yr⁻¹ remains effectively unchanged in our posterior at this precision. In contrast, emissions decrease over oceanic regions from 1.4 to 1.2 Tg yr⁻¹ (−14 %), largely reflecting reduced fossil fuel contributions (Fig. 7b). Waste-sector emissions show relatively minor adjustments at the regional and subregional scales, reflecting compensating regional contrasts. Posterior corrections are nevertheless evident in specific hotspot regions (Fig. 7c), with increases in BENELUX, France, and Poland, and decreases in the UK and Italy, highlighting substantial country-level heterogeneity. Wetland emissions in northern Europe were obtained to be 1.1 Tg yr⁻¹ with no discernible change at this precision. However, notable summer reductions occur in northern Europe, as well as in the UK, Ireland, Italy, and coastal areas, indicating minor overestimations in these regions. In contrast, increases are observed in Eastern Europe and BENELUX (see Fig. 7d). CH₄ emissions from geological seepage decline sharply in southern Europe, from 2.8 Tg yr⁻¹ in the prior estimate to 1.5 Tg yr⁻¹ in the posterior, representing a 46 % reduction, primarily driven by substantial prior overestimated emissions from Italy and Romania (Steiner et al., 2024) (see Fig. 7e). Details of regional and sectoral posterior fluxes estimates are provided in Appendix A Table A4. Total CH₄ emissions increase at the EU27+3 level, mainly in western and central Europe, while decreases occur in southern and oceanic regions (Fig. 7f). Overall, these posterior flux adjustments highlight that anthropogenic activity, especially emissions from agriculture and waste, drives most of the net increase in CH₄ emissions, while geological emissions decrease. These heterogeneous corrections underscore the importance of refining prior emission patterns at both sectoral and regional scales to develop robust methane budgets for Europe.

Figure 7 Mean estimates of prior and posterior CH₄ fluxes averaged over Europe and its subregions for the years 2017–2022. The panels display boxplots for total CH₄ emissions (Total), agriculture (AGR), energy (ENG), waste (WST), wetlands (WET), and geological sources (GEO). The Total emissions are derived from the total-flux inversion, whereas the sectoral emissions are obtained from the sector-specific inversion. Boxplots within the lightblue-shaded area of each subpanel (Ocean–Scandic) correspond to smaller regions and are plotted against the right-hand y axis. The boxplot ranges represent temporal variability over the study period.

3.3 Temporal Variation of CH₄ Flux

Figure 8 illustrates the seasonal patterns of CH₄ emissions from the Agriculture, Energy, Waste, and Wetlands sectors, spatially aggregated across the EU27+3 subregions. Among these, the Wetlands sector in Northern Europe exhibits the most pronounced seasonality, with emissions peaking during the summer months, particularly in July. This pattern is primarily driven by warmer temperatures and consistently sufficient soil moisture, which together enhance microbial activity and methane production in wetland ecosystems (Bechtold et al., 2025; Aalto et al., 2025). The accelerated decomposition of organic matter under these conditions further contributes to elevated emissions during this period (Voigt et al., 2023). The red hatched shading in the figure indicates a posterior reduction in summer wetland emissions, suggesting a downward correction of potentially higher prior estimates in summer. The Waste sector shows notable seasonal patterns, particularly in Central EU27+3 countries, where emissions increase significantly during summer. This is likely due to intensified microbial decomposition of organic waste in warmer temperatures, leading to greater methane generation in landfills and wastewater treatment facilities. In contrast, the Southern EU27+3 countries display a distinct summer minimum in waste emissions, highlighting regional differences that may be influenced by climate or waste management practices. Emissions from the Energy sector exhibit an opposite seasonal trend, with a clear peak during the winter and a minimum in summer. This pattern is largely attributed to increased heating demand in colder months, resulting in higher fossil fuel combustion and methane emissions. Additionally, winter energy use often involves greater reliance on natural gas, which can lead to methane leakage from pipelines and storage systems, further contributing to wintertime emission peaks. Recent studies support this seasonal behaviour. For example, Hu et al. (2025) found that oil and gas methane emissions can be substantially higher in winter, with estimates indicating increases of up to 40 % compared to summer. Similarly, Varon et al. (2025) reported pronounced seasonal variability in methane emissions from oil and gas production regions based on satellite inversions. These findings provide independent observational evidence supporting the wintertime enhancement of energy-sector methane emissions inferred in this study. These contrasting seasonal dynamics between sectors may counterbalance one another, contributing to the lack of a clear overall seasonal pattern in total CH₄ emissions at the regional scale. The Agriculture sector, by comparison, does not exhibit a distinct seasonal pattern across any of the regions. However, hatched blue shading in the figure indicates posterior increases in agricultural emissions over Central and Western EU27+3 countries, suggesting model adjustments based on observational constraints. These regional and sector-specific variations underscore the complexity of CH₄ emission dynamics and highlight the importance of disaggregated analyses to improve understanding and model representation of seasonal fluxes.

Figure 8 Seasonal patterns of CH₄ emission from major sectors, aggregated over the subregions for 2017–2022: (a) northern, (b) western, (c) central, and (d) southern Europe. The solid line represents the multi-year monthly mean, while the shaded area indicates the interannual variability (standard deviation). The hatched area shows the difference from the prior.

3.4 EU27+3 CH₄ Emission Estimates and Comparison with UNFCCC Reports

We calculated country-level anthropogenic CH₄ emission estimates for the EU27+3 and compared them with EDGARv8, GAINS, and UNFCCC NGHGI (2023) reports. Our posterior estimates averaged 17.6 Tg yr⁻¹ for 2017–2021 (range: 17.0–18.2), with a mean annual decline of 0.3 Tg yr⁻¹. These posterior estimates are 11 % higher than GAINS, 4 % higher than EDGARv8, and 3 % higher than UNFCCC, indicating good overall consistency at the EU27+3 scale. They are also in close agreement with previous inversion studies; for example, Steiner et al. (2024) reported 17.4 Tg yr⁻¹ for EU27+UK in 2018, compared to our estimate of 17.6 Tg yr⁻¹ for the same region and year. The top three emitting countries (France, Germany, and the UK) accounted for 39 % of EU27+3 emissions, similar to inventories but with differing country rankings. The UK stood second in the UNFCCC report. However, larger discrepancies appear at the national level (see Fig. 9 and Appendix A Table A5). Compared with UNFCCC, posterior emissions are higher for BENELUX (+54 %), Germany (+37 %), and France (+10 %), but lower for Romania (−25 %), Poland (−16 %), and Italy (−11 %). Notable differences also occur in the Nordics, with decreases in Norway (−39 %), Finland (−16 %), and Denmark (−13 %), but an increase in Sweden (+10 %). These national-scale differences are consistent with the findings of East et al. (2025), who compared inversion-derived methane emissions with UNFCCC inventories and similarly reported both positive and negative deviations at the country level across Europe. Their results underline the value of atmospheric inversion approaches for identifying potential biases and uncertainties in national methane emission reporting. Relative to EDGARv8, posterior estimates are higher for Portugal (+40 %), Italy (+18 %), and Germany (+10 %), but substantially lower in the UK (−18 %), Poland (−7 %), Romania (−18 %), and especially the Nordics. Adjustments relative to GAINS are upward for Germany (+29 %), France (+16 %), and BENELUX (+40 %), but downward for Switzerland (−13 %) and Italy (−7 %).

Figure 9 Annual mean anthropogenic CH₄ emissions (2017–2021) for selected EU27+3 countries, shown as bar plots for the posterior estimates, GAINS, UNFCCC, and EDGARv8. Subpanels display CH₄ emissions from major anthropogenic sectors: (a) total anthropogenic emissions excluding land-use and land-cover change, (b) agriculture, (c) waste, and (d) energy. For each country and sector, bars of the same color represent a given dataset and are stacked by year, from bottom (2017) to top (2021).

In addition to country-level discrepancies in total anthropogenic estimates, systematic patterns emerge across sectors, as displayed in Fig. 9. In agriculture, inventories generally underestimate emissions in large agricultural countries such as France, Germany, and the Netherlands, whereas our posterior consistently suggests higher values. This implies that bottom-up activity data and emission factors may not fully capture agricultural methane sources. In the waste sector, EDGARv8 tends to overestimate emissions in countries like Germany and France but underestimates them in others, including the UK, Italy, and Poland. Our posterior estimates often fall closer to GAINS and UNFCCC values. For the energy sector, inconsistencies are more pronounced. In the UK, EDGARv8 clearly overestimates emissions, whereas in Poland and Romania, both EDGARv8 and UNFCCC report higher values than our posterior. These differences are likely driven by variations in fossil fuel activity data and the emission factors employed by different reporting frameworks. Overall, the EU27+3 total appears relatively robust across inventories, yet sectoral and country-level comparisons reveal significant discrepancies. Our inversion estimates provide an essential independent constraint in this context, helping to reconcile inconsistencies in bottom-up reporting systems and strengthening confidence in national and sectoral methane emission estimates.

3.5 Sensitivity Experiments and Their Impact on Optimised Flux Estimates

To assess the sensitivity of the inversion framework, we conducted eight inversions (S1–S8) by varying key parameters, including background mole fractions, prior flux and background mole fraction error assumptions, observational coverage, choice of prior inventories, and correlation assumptions detailed in Sect. 2.5. Figure 10 presents the resulting posterior flux increments and error reductions for July 2021; this one-month analysis is intended to compare sensitivity across different setups rather than to provide a detailed uncertainty estimate. Inversion setups from S1 to S5 exhibit broadly consistent spatial patterns and magnitudes of posterior flux increments, though the degree of error reduction varies among cases. The base inversion setup (S1) achieves substantial reductions of 14 % for Europe and 21 % for EU27+3, with particularly strong constraints in Central Europe (38 %) and Western Europe (30 %) (see Table 4). S2 is nearly indistinguishable from S1 at both regional and subregional scales, showing only slightly weaker reductions in northern Europe, suggesting that posterior adjustments are relatively insensitive to the choice of background concentration. S3 and S4 yield somewhat weaker mean error reductions of about 16 % for EU27+3, indicating moderate sensitivity of posterior adjustments to prior uncertainty scaling and background error assumptions. Larger assumed background errors allow observations to exert a stronger influence, whereas smaller errors suppress corrections and increase reliance on the prior. Although these effects are modest in magnitude (−2 to +2 ppb), they highlight the critical role of background error settings in modulating inversion flexibility and posterior responsiveness. The influence of background error settings on concentration fields is illustrated in Fig. 11.

Table 4 Prior, Posterior, and Mean Error Reduction (%) for different regions and sensitivity inversion setups (S1–S8) as described in Sect. 2.5.

Figure 10 Flux increments (posterior – prior) and error reduction for July 2021 across eight sensitivity inversion setups (S1–S8). Each panel illustrates the spatial distribution of flux adjustments relative to the prior, highlighting the influence of varying prior emissions, background treatment, observational network density, and error correlation assumptions. The last subpanel displays the stations added for inversion setup S5. Mountain stations are shown as black dots, while low-altitude stations are shown as blue dots.

Figure 11 Impact of background error assumptions on optimized methane mole fractions and fluxes for 2021. Top row: (a) Prior mole fractions from CAMS and posterior increments with background error assumptions of (b) 0.5 % and (c) 0.05 %. Bottom row: (d) Time series of posterior mean increments under both settings and (e) relative annual flux difference between inversions with background errors of 0.5 % and 0.05 %. To quantify the impact of background-error assumptions, the relative difference was defined as $\mathrm{1}-\frac{F_{\mathrm{0.5}\mathit{\%}}}{F_{\mathrm{0.05}\mathit{\%}}}$, where F_0.5 % and F_0.05 % denote the flux computed using a background-error assumption of (0.5 %) and (0.05 %) respectively.

Network densification in S5 enhances error reduction across most regions, particularly in Western Europe, where reductions exceed 35 %. This improvement reflects the added influence of stations such as TAC in the UK and PDM and OHP in France, which help better constrain emissions over the UK and southern France. The additional mountain station in the Alps, however, does not show a clear enhancement in posterior uncertainty reduction, suggesting that its effectiveness may depend on representativeness or transport model resolution. Although the spatial pattern of posterior flux increments remains consistent with earlier inversion setups, the expanded network clearly strengthens the overall robustness of the inversion. These findings are consistent with previous studies demonstrating that increased observational density improves posterior reliability in well-sampled regions (Villani et al., 2010; Thompson and Stohl, 2014).

More pronounced differences emerge from the S6–S8 inversion setup. Using EDGARv8 as the prior (S6) increases the prior flux magnitude to 26.1 Tg yr⁻¹ for EU27+3, compared to 24.5 Tg yr⁻¹ in S1, reflecting EDGARv8's higher reported emissions. However, the posterior flux in S6 decreases substantially to 23.0 Tg yr⁻¹, closely aligning with the posterior estimate from the base inversion (S1). The 19 % mean error reduction in S6 is also comparable to S1, highlighting the sensitivity of inversion results to the choice of emission inventory and the convergence of posterior estimates despite differing priors. In contrast, reducing correlation lengths leads to diminished error reduction. At 50 km (S7), reductions are minimal (1 %–3 %), while at 20 km (S8) they are negligible (<0.2 %). This reflects the trade-off between allowing fine-scale flux variability and weakening the effective observational constraint per degree of freedom. Such behaviour is consistent with Thompson and Stohl (2014), who showed that extreme localization increases small-scale heterogeneity at the cost of inversion robustness.

Despite these sensitivities, the dominant spatial patterns of posterior corrections remain robust across sensitivity inversion setups: decreases over Italy, Romania, and the UK, and increases over BENELUX, Germany, and France. This consistency with previous European inversion studies (Saunois et al., 2020; Steiner et al., 2024) underscores the resilience of large-scale flux signals, even under varying methodological assumptions. Overall, S1 emerges as a balanced configuration, providing a reasonable compromise between observational constraint and inversion flexibility. Alternative setups highlight sensitivities to priors, background errors, and correlation lengths, but the preservation of large-scale spatial patterns across experiments supports the robustness of the main conclusions.

Figure 11 illustrates the background concentration changes and their sensitivity to background error assumptions. Under the larger background error scenario (0.5 % uncertainty), we found a systematic decrease in background CH₄ over the Northern and Mediterranean Seas, suggesting that the prior background fields were biased high in these regions (see Fig. 11b). The spatially averaged time series of the posterior increment exhibits a clear seasonal pattern, occurring from May to October (see Fig. 11d), reflecting the correction of this high bias through the inversion. Conversely, assuming a smaller error (0.05 %) resulted in much weaker corrections, both spatially (see Fig. 11c) and temporally (see Fig. 11d), as seen in the posterior mean increment time series. Under a low-error assumption, the inversion relies more strongly on the initial concentration estimates, limiting adjustments to the posterior background concentration field. Conversely, higher error assumptions allow the observations to exert greater influence, producing stronger corrections in regions where the prior background fields were biased.

Figure 11e shows the relative annual mean difference of the posterior fluxes under background error assumptions of 0.5 % and 0.05 %. The relative difference is defined as $\mathrm{1}-F_{\mathrm{0.5}\mathrm{\%}}/F_{\mathrm{0.05}\mathrm{\%}}$, where F_0.5 % and F_0.05 % denote the annual mean fluxes computed using background-error assumptions of 0.5 % and 0.05 %, respectively. Flux adjustments under these scenarios generally follow the spatial patterns of the mean map (see Fig. 5c), with additional corrections of approximately −20 % in regions showing decreases and +20 % in regions showing increases. An extreme reduction (−40 %) occurs over Italy, likely contributing to the observed posterior flux decrease. Overall, background-related adjustments remain modest (−2 to +2 ppb) and are secondary to the dominant sectoral flux corrections.

4 Conclusions

In this study, we have presented a high-resolution, top-down estimates of European CH₄ emissions for the period 2017–2022, offering robust new evidence on their magnitude, spatial structure, and temporal evolution across the continent. By combining the CIF-FLEXPART inversion framework with a 4D-Var data assimilation approach, we demonstrate the value of advanced atmospheric inverse modelling for resolving emission patterns in high-resolution inversion settings, particularly where large state vectors make conventional analytical inversion approaches computationally challenging. The use of FLEXPART-derived source–receptor sensitivities, together with the optimisation of a large set of control variables, enables the inversion of emissions at a monthly temporal resolution and a spatial resolution of 0.2°×0.2°. The assimilation of observations from a dense network of 46 in situ monitoring stations substantially strengthens the observational constraints over Europe, leading to improved confidence in the inferred emissions and their sectoral attribution. By separately optimising emissions from major source categories, including agriculture, energy-related activities, waste management, industrial processes, wetlands, geological sources, and biomass burning, this work provides actionable information that is directly relevant to national inventories and mitigation strategies.

The posterior concentration in our study shows a great improvement with atmospheric observations (r²=0.87, RMSE=24.4 ppb, $\text{mean bias}=-\mathrm{2.1}$ ppb) and captures most of the observed CH₄ variability. Our results reveal substantial regional-scale adjustments relative to prior inventories. Posterior emissions are consistently higher over the BENELUX region, France, and Germany, suggesting underestimation in current inventories, whereas posterior reductions over the UK, Italy, and Romania indicate possible prior overestimation. These regional patterns are broadly consistent with earlier inversion studies (Bergamaschi et al., 2022; Saunois et al., 2020; Petrescu et al., 2021; Steiner et al., 2024), which lends confidence to their robustness. At the aggregated EU27+3 scale, our total posterior methane emissions are corrected from 21.7±2.7 to 23.3±2.3 Tg yr⁻¹, representing a modest 6.6 % rise compared with the prior. Although this is a modest overall increase, there are strong regional and sectoral differences. Anthropogenic CH₄ emissions are averaged 17.6 Tg yr⁻¹ during 2017–2021 (range: 17.0–18.2 Tg yr⁻¹), corresponding to a mean annual decrease of 0.3 Tg yr⁻¹. These values are 11 % higher than GAINS, 4 % higher than EDGARv8, and 3 % higher than UNFCCC NGHGI (2023), indicating good overall consistency with bottom-up inventories and recent inversion-based estimates (e.g. Steiner et al., 2024). However, country-level comparisons reveal notable discrepancies: emissions are higher than reported for BENELUX (+54 %), Germany (+37 %), and France (+10 %), but lower for Romania (−25 %), Poland (−16 %), and Italy (−11 %). Notable decreases are also found for Norway (−39 %), Finland (−16 %), and Denmark (−13 %), while Sweden exhibits a modest increase (+10 %), pointing to areas where inventory improvements are needed. Furthermore, our sectoral analysis indicated that agriculture is the dominant source of European methane emissions and the primary driver of the posterior increase, accounting for nearly half of total EU27+3 emissions and corrected upward by 21 % relative to the prior estimate. Energy-sector emissions show a small increase over land regions, while decreasing in coastal and ocean areas, and waste emissions remain broadly stable, despite pronounced national variability. In contrast, emissions from wetlands and geological sources show decreases, particularly reductions in wetland emissions observed during the summer months in Northern Europe and decreases in geological emissions in Italy, Romania, and the UK. In addition, we conducted eight sensitivity experiments by varying key parameters in the inversion setup. The results indicated that posterior fluxes show limited sensitivity to the choice of initial mole fractions and prior fluxes, underscoring the dominant role of atmospheric observations in constraining emissions. In contrast, assumptions regarding error correlation lengths and observational density strongly influence uncertainty reduction and the magnitude of posterior fluxes, emphasising the importance of realistic error characterisation and the continued expansion of atmospheric measurement networks.

Our study demonstrates that dense atmospheric monitoring combined with high-resolution inversion modelling provides a robust top-down framework to evaluate national inventories, identify regional emission hotspots, and support verification of methane mitigation efforts in Europe. By improving both spatial resolution and sectoral attribution, these results advance understanding of European methane sources and their variability, with direct relevance for climate policy and mitigation verification. Nevertheless, some limitations remain. Observational coverage is sparse in parts of eastern and southeastern Europe, limiting the ability to resolve emissions at finer spatial scales in these regions. Uncertainties in the transport model and assumptions about background CH₄ levels may also affect the posterior estimates. In addition, this study didn't systematically assess the impacts of different resolutions, grid-based uncertainty, different natural priors and OH fields. Thereby, future work should address these limitations through dedicated sensitivity analyses and by integrating satellite observations and isotopic signatures (e.g. ¹³CH₄) to improve source attribution, expand spatial coverage, and further strengthen European methane emission estimates.

Appendix A: Additional Figures and Tables

Figure A1 Mean spatial distribution of prior methane (CH₄) fluxes over the inversion domain for 2017–2022. Prior emissions are shown by sector: (a) Agriculture (AGR), (b) Energy (ENG), (c) Industrial processes (IND), (d) Waste (WST), (e) Ocean (OCE), (f) Wetlands and soil sinks (WET), (g) Biomass burning (BBR), (h) Geological (GEO), and (i) Termites (TER). Anthropogenic prior fluxes were derived monthly from GAINS (for EU27+3) and EDGARv8 (for the rest of the world). The GAINS sectors contributing to the AGR, ENG, and WST categories include: Energy: A_PublicPower, B_Industry, D_Fugitives, F_RoadTransport, C_OtherStationaryComb; Agriculture: K_AgriLivestock, L_AgriOther; Waste: J_Waste. The EDGARv8 sectoral inputs correspond to IPCC categories: Energy: BUILDINGS, TRANSPORT, IND_COMBUSTION, POWER_INDUSTRY, FUEL_EXPLOITATION; Industrial processes: IND_PROCESSES; Agriculture; AGRICULTURE; Waste: WASTE.

Table A1 Subregional grouping of European countries used in this study. ISO codes follow ISO 3166-1 alpha-3.

Figure A2 Time series of CH₄ flux contributions from the nested and global domains at selected stations, averaged at the assimilated hour for 2021: BIS (Biscarrosse, France), CBW (Cabauw, Netherlands), CRP (Carnsore Point, Ireland), ERS (Ersa, coastal Mediterranean site), HPB (Hohenpeissenberg, mountain station in Germany), KIT (Karlsruhe, Germany), PAL (Pallas, Finland), and TRN (Trainou, France).

Figure A3 FLEXPART footprint sensitivity to the assimilated observations, shown as the mean sensitivity field for July 2021, illustrating the spatial heterogeneity of observational constraints across Europe. Regions with denser measurement coverage exhibit stronger sensitivities.

Figure A4 Emission sensitivity, defined as the product of the footprint sensitivity and the monthly mean prior CH₄ flux for each emission sector, indicating where prior emissions and observational sensitivity jointly contribute to constraining the inversion.

Table A2 List of CH₄ concentration observation sites assimilated in this study. The “Alt” represents the sum of the surface elevation and the intake height a.g.l.. The “Usage” column uses “M” to denote stations included in the base inversion and “S” to denote stations used only in the sensitivity test. An asterisk (*) in the three-letter station code indicates the station is classified as a mountain station, while a double asterisk (**) indicates the station is classified as a coastal station. Time series of CH₄ concentration are obtained from the ICOS European ObsPack compilation, except for SOD and KMP, which are provided by FMI.

All links were accessed on 20 June 2026.

Table A3 Prior and posterior model performance statistics at all measurement stations used in this study. The table reports the root-mean-square error (RMSE), mean bias (MBias), and Pearson correlation coefficient (r²) for both the prior and posterior simulations at the assimilated hour. RMSE and MBias are given in ppb.

Table A4 Six-year mean methane fluxes (2017–2022) by region and sector. Values are given in Tg yr⁻¹ with ranges in parentheses.

Table A5 Comparison of posterior anthropogenic CH₄ emissions (Tg yr⁻¹) for EU27+3 countries with GAINS, UNFCCC, and EDGARv8 inventories, including relative differences (%).