Estimating methane emissions in the Arctic nations using surface observations from 2008 to 2019

The Arctic is a critical region in terms of global warming. Environmental changes are already progressing steadily in high northern latitudes, whereby, among other effects, a high potential for enhanced methane (CH₄) emissions is induced. With CH₄ being a potent greenhouse gas, additional emissions from Arctic regions may intensify global warming in the future through positive feedback. Various natural and anthropogenic sources are currently contributing to the Arctic's CH₄ budget; however, the quantification of those emissions remains challenging. Assessing the amount of CH₄ emissions in the Arctic and their contribution to the global budget still remains challenging. On the one hand, this is due to the difficulties in carrying out accurate measurements in such remote areas. Besides, large variations in the spatial distribution of methane sources and a poor understanding of the effects of ongoing changes in carbon decomposition, vegetation and hydrology also complicate the assessment. Therefore, the aim of this work is to reduce uncertainties in current bottom-up estimates of CH₄ emissions as well as soil oxidation by implementing an inverse modelling approach in order to better quantify CH₄ sources and sinks for the most recent years (2008 to 2019). More precisely, the objective is to detect occurring trends in the CH₄ emissions and potential changes in seasonal emission patterns. The implementation of the inversion included footprint simulations obtained with the atmospheric transport model FLEXPART (FLEXible PARTicle dispersion model), various emission estimates from inventories and land surface models, and data on atmospheric CH₄ concentrations from 41 surface observation sites in the Arctic nations. The results of the inversion showed that the majority of the CH₄ sources currently present in high northern latitudes are poorly constrained by the existing observation network. Therefore, conclusions on trends and changes in the seasonal cycle could not be obtained for the corresponding CH₄ sectors. Only CH₄ fluxes from wetlands are adequately constrained, predominantly in North America. Within the period under study, wetland emissions show a slight negative trend in North America and a slight positive trend in East Eurasia. Overall, the estimated CH₄ emissions are lower compared to the bottom-up estimates but higher than similar results from global inversions.

1 Introduction

The Arctic is an especially critical area in terms of global warming. As the near-surface air temperature has increased by approximately 3.1 ^∘C since the 1970s, 3 to 4 times as much as the global average (AMAP, 2021; Rantanen et al., 2022), environmental changes in that region are rapidly progressing (Serreze et al., 2009; Cohen et al., 2014; Jansen et al., 2020). Exceptional events like melting glaciers, a reduction in sea ice, thawing permafrost, an increasing occurrence of wildfires during summer and a shortening of the snow season have already been observed increasingly frequently during the most recent years (Hassol, 2004; Stroeve et al., 2007; Walker et al., 2019). Predictions assume that if the Arctic warming continues rising at this rate, by 2100 the temperature will have increased by 3.3 to 10.0 ^∘C (AMAP, 2021).

Short-lived climate forcers such as methane (CH₄) play a significant role in this framework (AMAP, 2015). Methane is globally the second most abundant anthropogenic greenhouse trace gas with a radiative forcing of about 0.56 W m⁻² (IPCC, 2022). The rising temperatures, at the global scale and particularly in the Arctic, influence the natural CH₄ sources in the Arctic, which may possibly intensify local emissions in the near future (IPCC, 2022). A positive feedback of the global – and regional – warming may therefore ensue.

Various CH₄ sources, both natural and anthropogenic, contribute to the Arctic methane budget. Today, the natural Arctic methane emissions are dominated by high-latitude wetlands, the extent of which is still highly uncertain however. Estimations on high-latitude wetland emissions show large discrepancies. Ito (2019) concluded from a process-based modelling study that the pan-Arctic (above 60^∘ N) wetland emissions in the 2000s were between 10.9 and 11.4 Tg CH₄ yr⁻¹. Estimates by Petrescu et al. (2010) of northern wetland emissions (defined as wetlands in regions with a yearly average temperature lower than 5 ^∘C) varied by a factor of 4 (between 38 and 157 Tg yr⁻¹), with the corresponding regions varying by a factor of 2 (2.2×10⁶–4.4×10⁶ km²). Uncertainties in the extent of high-latitude wetland areas are, among other factors, a reason for the large variations. Other natural CH₄ sources occurring in this area are freshwater emissions, e.g. from thermokarst lakes, as well as emissions from the Arctic Ocean and biomass burning due to wildfire events in the summer months (AMAP, 2015). As mentioned before, natural methane emissions are anticipated to increase with rising temperatures and overall changing conditions: in the Arctic, methane net emissions could possibly be twice as high by the end of this century (Schuur et al., 2015), in part related to the high sensitivity of CH₄ emissions to the state of the permafrost (Masyagina and Menyailo, 2020) and general atmospheric conditions (Chen et al., 2015). Indeed, the thawing and destabilization of permafrost lead to the exposure of large carbon pools that have so far been shielded by ice and frozen soil. Permafrost thaw is expected to influence at least four ways of carbon mobilization: (i) the deliberation of CH₄ reservoirs in the upper permafrost layers, (ii) retained activity from viable methanogens as well as (iii) the consumption of labile organic matter by these micro-organisms, and finally (iv) an increased production of CH₄ in the active zone (Rivkina and Kraev, 2008). Additionally, anthropogenic activities in high northern latitudes contribute to the global methane budget, with an estimated amount ranging between 2 and 18 Tg CH₄ yr⁻¹ (Saunois et al., 2020a). These emissions are mainly caused by the exploitation and distribution of fossil fuels and are especially predominant during the winter months (Thonat et al., 2017). Currently, five Arctic nations, Russia, Canada, Norway, Greenland and the United States of America, perform drilling activities in their territories and exclusive economic zones in neighbouring oceans. Decreasing the emissions from anthropogenic sources is an effective way to limit the overall methane emissions in the Arctic region. However, with an estimated 13 % of undiscovered mineral oil and 30 % of undiscovered gas resources north of the Arctic Circle (Gautier et al., 2009), the Arctic is of significant interest for the petroleum industry regarding future drilling campaigns.

Even though the CH₄ observation networks in northern high latitudes have been expanded since the early 2000s, the current stationary networks remain restricted, leaving vast areas uncovered due to the difficulties in carrying out measurements in such remote areas (Pallandt et al., 2022). Thus, obtaining accurate assessments of methane emissions in northern high latitudes remains challenging since their spatial distribution at the local scale is highly variable. Current estimations are primarily based on bottom-up studies which rely on upscaling local flux measurements or on process-based surface models and on emission inventories which combine emission factors with socio-economic activity data. These approaches are however subject to high uncertainties at the regional scale since they imply statistical approximations as well as simplifications of chemical, biological and physical processes (e.g. Saunois et al., 2020a).

Another approach is provided by top-down studies, in support of bottom-up products. Top-down studies optimally combine observations, provided by either ground-based or satellite measurements of atmospheric CH₄ mixing ratios; numerical transport modelling; and bottom-up emission data sets as prior emission estimates into the mathematical framework of data assimilation to retrieve emission fluxes and their uncertainties. The so-called atmospheric inversion method is therefore useful to reduce uncertainties in bottom-up estimates (used as priors) and thus gain a better understanding of the region's methane budget. Such studies have already been implemented for high-latitude regions at various scales and with regard to different sources. Inverse modelling approaches for methane emissions have for instance been carried out for the Canadian Arctic by Miller et al. (2016) (for the years 2005–2006), Ishizawa et al. (2019) (for the years 2012 to 2015), Chan et al. (2020) (for the years 2010 to 2017) and Baray et al. (2021) (for the years 2010 to 2015); for Scandinavia (Tsuruta et al., 2019, for the years 2004 to 2014); for high-latitude Eurasian regions (Berchet et al., 2015b, for the year 2010); for the Siberian Lowlands (Winderlich, 2012); and also for the whole region above 45^∘ N (Tsuruta et al., 2023, for the year 2018), above 50^∘ N (Thompson et al., 2017, for the years 2005 to 2013) and above 60^∘ N (Tan et al., 2016, for the year 2005).

In this study, we estimate methane emissions during the most recent years (2008 to 2019) through atmospheric inversion based on available in situ measurement data from observation sites located in the Arctic and sub-Arctic. In order to obtain a reliable assessment, we compute a large ensemble of possible posterior emission scenarios using different error estimations that are evaluated concerning their plausibility. The CH₄ emissions are subsequently analysed with particular regards to three different questions: (i) Is the available observation network sufficient to constrain all occurring CH₄ sources and sinks adequately? (ii) Do the different CH₄ sources and sinks show any significant trends between the years 2008 and 2019? (iii) Do the different CH₄ sources and sinks in the posterior state show any shifts in the seasonal cycle in comparison to the prior bottom-up estimates?

2 Methodology

To estimate the CH₄ fluxes in the Arctic region, a Bayesian inversion framework (Sect. 2.1) based on backward simulations of the Lagrangian particle dispersion model (LPDM) FLEXPART (FLEXible PARTicle dispersion model) is used (see details in Sect. 3.3). The inversion is based on all available observation sites in the Arctic and sub-Arctic region (see details in Sect. 2.3). Extensive sensitivity tests are carried out to evaluate the reliability of CH₄ estimates (see details in Sect. 2.2)

2.1 Inversion framework

We apply an analytical inversion which aims at explicitly and algebraically finding the optimal posterior state of a system x^a and the corresponding uncertainties P^a, which are given by

$$\begin{array}{}\text{(1)}& \left\{\begin{array}{rcl}{\mathit{x}}^{\mathrm{a}}& =& {\mathit{x}}^{\mathrm{b}}+\mathbf{K}({\mathit{y}}^{\mathrm{o}}-\mathbf{H}{\mathit{x}}^{\mathrm{b}})\\ {\mathbf{P}}^{\mathrm{a}}& =& \mathbf{B}-\mathbf{KHB}\end{array}\right.,\end{array}$$

with K the Kalman gain matrix given by

$$\begin{array}{}\text{(2)}& \mathbf{K}={\mathbf{BH}}^T(\mathbf{R}+{\mathbf{HBH}}^T{)}^{-\mathrm{1}}.\end{array}$$

In the present work, we apply the formula on a year-by-year basis. The control vector x^b refers to the prior knowledge on the system, in our case CH₄ surface fluxes from different sources (Sect. 3.2) and background mixing ratios (Sect. 3.3.3). The observation vector y⁰ contains the available observations of atmospheric CH₄ mixing ratios (detailed in Sect. 3.1.2). The observation operator includes not only the transport of the emitted methane (Sect. 3.3.2) in the domain and the import from outside the domain (Sect. 3.3.3) but also the filtering and other operations required to extract the simulated equivalents of the measurements (Sect. 3.3). We neglect chemical oxidation of CH₄ emitted in our regional domain by OH as further explained in Sect. 3.3, although oxidation by OH is still accounted for in the global model used for the background (Sect. 3.3.3).

Thus, all operations in the observation operator are linear, and we represent it by its Jacobian matrix H. The linear assumption is required to write Eq. (1) and solve the Bayesian system analytically.

The error covariance matrices in the observation and control spaces, R and B, define the weight of the mismatch between the modelled and the measured concentrations. R contains various types of errors: the error estimates of the differences between the observations and their simulated equivalents include uncertainties not only in the measurements but also in the transport in the model and in the discrete representation of the continuous world by a numerical model. The dimensions of R are equivalent to the number of elements in the observation vector per year; it varies between 217 and 384 as observations are aggregated by station and month (see Sect. 3.1.2). The covariance matrix B is composed of two parts: B^S which accounts for the uncertainties in the prior methane fluxes and B^B for the uncertainties in the background mixing ratios. B^S has a constant size of 10 164 × 10 164, following the number of emission regions, emission sectors and emission periods optimized in our system (see Sect. 3.2); the dimensions of B^B are, again, equivalent to the number of observations per year.

Defining the error covariance matrices can be challenging since only the measurement uncertainties can be determined with certainty, using rigorous calibration procedures (e.g Sasakawa et al., 2010). On the other hand, unrealistic error estimations can drastically distort the results of the posterior state (Berchet et al., 2013). Therefore, in this study an ensemble of $({\mathit{x}}_i^{\mathrm{a}}{)}_{i=\mathrm{1},\mathrm{500}}$ and $({\mathbf{P}}_i^{\mathrm{a}}{)}_{i=\mathrm{1},\mathrm{500}}$ using 500 realistic set-ups of the error matrices (R,B) is computed. The ensemble of $(\mathbf{R},\mathbf{B}{)}_{i=\mathrm{1},\mathrm{500}}$ pairs of matrices is described in Sect. 3.1.2 and in Sect. 3.2.2, respectively. To account for the uncertainties in the posterior state, from each vector ${\mathit{x}}_i^{\mathrm{a}}$, 10 random variations are generated with the corresponding covariance matrix ${\mathbf{P}}_i^{\mathrm{a}}$ following a multivariate normal distribution. Thus, we obtain a total of 5000 posterior states to assess the posterior uncertainties in the inversion.

For computational reasons, the 12-year period has been split into 12 independent 1-year inversion windows computed separately. The ensemble of 500 pairs of matrices $(\mathbf{R},\mathbf{B}{)}_{i=\mathrm{1},\mathrm{500}}$ is generated based on a limited number of parameters independent from the year j (see Sect. 3.1.2 and 3.2.2). Therefore, for a given member i of the ensemble, the yearly matrices {(R^j,B^j)_i} are built on the same set of underlying parameters. We then compute, for each year $j\in [\mathrm{2008},\mathrm{2019}]$, 500 independent inversions.

2.2 Framework evaluation

2.2.1 Log likelihood of samples

Though realistically chosen (see Sect. 3.2.2 and 3.1.2), the members of the Monte Carlo ensemble of (R,B) pairs are not equally plausible. To further compare and aggregate statistics on the ensemble, we weight each member i $\in [\mathrm{1},\mathrm{500}]$ for each year j $\in [\mathrm{2008},\mathrm{2019}]$ by its likelihood (see e.g. Michalak et al., 2005). It is defined by

$$\begin{array}{}\text{(3)}& \begin{array}{rl}\ln p_i^j& \left({\mathbf{R}}_i^j,{\mathbf{B}}_i^j|{\mathit{y}}_j^{\mathrm{o}},{\mathit{x}}_j^{\mathrm{b}},{\mathbf{H}}_j\right)=-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\mathrm{tr}\left({{\mathbf{S}}_{{\mathbf{R}}_i^j,{\mathbf{B}}_i^j}^j}^{-\mathrm{1}}{\mathbf{S}}^j\right)\\ & -{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\ln \left|{\mathbf{S}}_{{\mathbf{R}}_i^j,{\mathbf{B}}_i^j}^j\right|,\end{array}\end{array}$$

with ${\mathbf{S}}^j={\mathit{y}}_j^{\mathrm{o}}-{\mathbf{H}}_j{\mathit{x}}_j^{\mathrm{b}}$ and ${\mathbf{S}}_{{\mathbf{R}}_i^j,{\mathbf{B}}_i^j}^j={\mathbf{R}}_i^j+{\mathbf{H}}_j{\mathbf{B}}_i^j{\mathbf{H}}_j^T$, where $|\cdot |$ is the determinant operator and tr(⋅) is the trace function.

The estimation of the log likelihood provides a robust method to select the most reliable set-ups, with regards to the information provided by the observations and ideal statistics (e.g. Winiarek et al., 2011; Berchet et al., 2013, 2015a). For a given set-up, the higher the log likelihood, the more plausible the pair of covariance matrices. The log-likelihood estimator in a high-dimension problem like ours is extremely sensitive to any change in configuration.

The range of the log likelihood varies between the different years, due to variations in the number of available sites and measurements, as well as atmospheric conditions. Then, for each member of the Monte Carlo ensemble, we define the cumulative log likelihood as

$$\begin{array}{}\text{(4)}& \ln p_i=\sum _{j=\mathrm{2008}}^{\mathrm{2019}}\ln p_i^j.\end{array}$$

We use the cumulative log likelihood to define the most plausible posterior vector over the full period of interest from 2008 to 2019 as ${\mathit{x}}_{\max }^{\mathrm{a}}$ corresponding to the member i_max maximizing the cumulative log likelihood.

We also use the log likelihood to discard the less realistic members of the Monte Carlo ensemble. To do so, the most reliable pair $i_{\max }^j$ of error matrices $({\mathbf{R}}_{\max }^j,{\mathbf{B}}_{\max }^j)$ is determined for each year j separately. Then, each optimal member $i_{\max }^j$ for year j is used on all the years of interest $j^{\prime }\in [\mathrm{2008},\mathrm{2019}]$, so as to obtain the corresponding cumulative log likelihood $\ln p_{i_{\max }^j}$.

Since each cumulative log likelihood $\ln p_{i_{\max }^j}$ includes the most reliable configuration for year j, the lower threshold for the log likelihood ln p_min is defined as the minimum of the 12 thus computed cumulative log-likelihood values as $\underset{j\in [\mathrm{2008},\mathrm{2019}]}{min}\ln p_{i_{\max }^j}$. We define a sub-ensemble $\mathit{\left\{}{\mathit{x}}_{\max }^{\mathrm{a}}\mathit{\right\}}$ whose elements have a cumulative log likelihood greater or equal to this threshold: $\mathit{\left\{}{\mathit{x}}_i^{\mathrm{a}}\right|{\sum }_{j=\mathrm{2008}}^{\mathrm{2019}}\ln p_i^j>\ln p_{\min }\mathit{\}}$. This sub-ensemble contains 274 configurations, which correspond to 2740 posterior states, and is used in the following for a representative analysis of the posterior state.

2.2.2 Sensitivity and influence matrices

We use two other metrics to evaluate our system and the different set-ups: the influence and the sensitivity matrices. Both are calculated using the corresponding Kalman gain matrix K_max of the previously determined ${\mathit{x}}_{\max }^{\mathrm{a}}$. The influence matrix K_maxH (defined by Cardinali et al., 2006), also called the averaging kernel (Rodgers, 2000), contains diagonal terms between 0 and 1, which represent the sensitivity of each component of x to the inversion. The smaller the term K_maxH_r,r for emissions in region r is, the less constrained region r is by the inversion. The sensitivity matrix HK_max (Cardinali et al., 2006) gives the sensitivity of the inversion to a change in one component of the observation vector. An observation with a high sensitivity brings strong constraints on the inversion. The weight of each station in the inversion can be computed by summing up the corresponding diagonal elements of HK_max. The trace of these two matrices also gives the “degrees of freedom for signal” (Wahba et al., 1994; Cardinali et al., 2006), while the number of observations minus this number gives the “degree of freedom for noise”. This extra criterion informs how many observations are used to constrain fluxes (and background mixing ratios).

2.3 Area and period of interest

The area of interest, shown in Fig. 1a, for this study regarding the quantification of the methane fluxes includes the Arctic and sub-Arctic, with the southern boundary being roughly the southernmost border of the taiga. For the implementation of the inversion, only observation sites within the area of interest have been included in this study. To represent concentrations at these sites as properly as possible, we simulate the influence of fluxes from the area of interest and also from a buffer region from above 30^∘ N (see Sect. 3.3.1). Even though Arctic fluxes may influence observation sites in the buffer region, we do not include them in this study due to the increased computational costs this would induce; future work may inquire into the impact of using as many stations as possible.

Figure 1 Area of interest (a) and RECCAP regions above 30^∘ N (b). Measurement sites, listed in Table S1 (Supplement), are indicated with white stars.

The region above 30^∘ N is subsequently divided into sub-regions in order to better detect local differences. However, the sub-regions should not be too small and numerous, due to the limitation of available observations for constraining those areas. A more detailed description of the selected observation sites (indicated with white stars in Fig. 1b) can be found in Sect. 3.1.1. The sub-regions of this study are therefore selected following the proposition of the Regional Carbon Cycle Assessment and Processes (RECCAP; Ciais et al., 2022), which results in 121 regions within the area of interest (Fig. 1b).

The time period of interest is from 2008 to 2019. For the following years, no measurements were available for the majority of the measurement sites by the time this study was implemented.

The atmospheric sites in the area and time period of interest and the available observations are described in Sect. 3.1. CH₄ emissions in this area are described in Sect. 3.2.

3 Material

3.1 Atmospheric observations

3.1.1 Site description

For this study, both quasi-continuous measurements (35 observation sites providing hourly measurements) and discrete measurements (6 observation sites providing task samples two to four times per month) are used. The stations are exclusively located in seven Arctic nations (Canada, Russia, Finland, Norway, Iceland, Greenland and the USA), except for one site in Ireland that is used to constrain air masses from the Atlantic Ocean. The operators of these stations are Environment and Climate Change Canada, the Japan–Russia Siberian Tall Tower Inland Observation Network (JR-STATION) of the National Institute for Environmental Studies (NIES) (Sasakawa et al., 2010), the US National Oceanic and Atmospheric Administration Global Monitoring Laboratory (NOAA GML;  Dlugokencky et al., 2020), and the Finnish Meteorological Institute (FMI;  Hatakka et al., 2003; Aalto et al., 2007). The stations with their trigram identifications are shown in Fig. 2.

All the measurement sites are subsequently described briefly, sorted by their network operators. A summary of each station's characteristics is furthermore provided in the Supplement (Table S1).

ECCC

Environment and Climate Change Canada (ECCC) established its first 2 CH₄ measurement stations (ALT and FRD) at the end of the 1980s and has expanded its network to 22 sites to this date, with 12 of them being located in the Arctic or sub-Arctic. Alert (ALT) is often referred to as an Arctic background site since it is located remotely from any major methane emission sources on the northeastern tip of Ellesmere Island in Nunavut, where the land is covered with snow for approximately 10 months a year. Two additional sites are installed in Nunavut at slightly more southern latitudes: Cambridge Bay (CBY) and Baker Lake (BLK). The latter is located in the Arctic tundra around 320 km from Hudson Bay, surrounded by small lakes, whereas CBY lies on the southeastern coast of Victoria Island close to the largest port of the Northwest Passage of the Arctic Ocean.

Figure 2 Map of available observation sites in the area of interest. Crosses indicate quasi-continuous measurements; diamonds indicate discrete measurements. All sites were used for the study except CDL, where no measurements were available during the study period. Different network operators are marked with different colours. NOAA ESRL: NOAA Earth System Research Laboratories.

The measurement site Inuvik (INK) was established in the Arctic tundra of the Northwest Territories in the east channel of the Mackenzie delta. Further inland in the same Canadian province lies the station BCK, 10 km from the town of Behchoko and surrounded by mixed forests, lakes and ponds.

Three of the ECCC sites are located in British Columbia. FNE, which is located near the small town of Ford Nelson in the taiga, lies at the southern fringe of the Canadian permafrost region. Estevan Point (ESP) is located on the coast of the Pacific Ocean and surrounded by woodlands. The measurement station Abbotsford (ABT) lies close to the US border, 80 km from Vancouver, the largest city and main economic area in British Columbia.

The two sites in the province of Alberta are LLB at Lac La Biche, in a region of peatlands and forests, and Esther (EST), which lies in the open prairie with plenty of cattle ranches close by.

Two measurement stations have been established in Saskatchewan. East Trout Lake (ETL) in the centre of the province lies at the southern edge of a boreal forest region, and Bratt's Lake (BRA) is in the Canadian prairie.

Churchill (CHU) is located in Manitoba, north of the largest continuous boreal wetland region in North America on the western coast of Hudson Bay.

Four of the sites in the province of Ontario (EGB, DWN, HNP, TKP) are located relatively close to each other in the Mixedwood Plains Ecozone. Downsview (DWN) and Hanlan's Point (HNP) are urban stations in the north of Toronto and on the Toronto Islands in Lake Ontario, respectively. Egbert (EGB) lies around 80 km from Toronto, close to a rural village. The southernmost site, Turkey Point (TKP), is located on Lake Erie in a woodland area. Further north in Ontario lies the station Fraserdale (FRD) in the boreal forest, with extensive wetland coverage in the surrounding area.

The two sites located in Quebec, Chapais (CPS) and Chibougamau (CHB), are likewise established close to each other in an area dominated by boreal forest with many lakes.

Finally, the observation site Sable Island (WSA) is on a remote island in the North Atlantic Ocean, 175 km from the mainland. The island is uninhabited by people and covered with grass and low-growing vegetation.

NOAA GML

The two continuous measurement stations operated by NOAA GML are Barrow (BRW) and CARVE (Carbon in Arctic Reservoirs Vulnerability Experiment, CRV) in the USA (Dlugokencky et al., 2020). Methane measurements in BRW started in the late 1980s. The site is located in northern Alaska on the junction of the Chukchi Sea and Beaufort Sea, and the surrounding landscape is characterized by thermokarst lakes. The CRV tower is located in boreal Alaska with a surrounding landscape defined by evergreen forest, shrubland and some areas of woody wetlands (Karion et al., 2016).

The six discrete measurement sites operated by NOAA GML are ZEP, SUM, ICE, MHD, CBA and SHM. The Zeppelin Observatory (ZEP) is located near the village of Ny-Ålesund, which is surrounded by mountains and glaciers, on the island of Spitsbergen. From 2017, ZEP observations are available as continuous data via the Integrated Carbon Observation System (ICOS) Carbon Portal (Lund Myhre et al., 2022), but we did not include them as such to avoid perturbing the interpretation of the results for the last years. The sampling site Summit (SUM) was established on the Greenland ice sheet and is the highest measurement site in the Arctic Circle. Stórhöfði (ICE) lies in the south of Iceland at the top of a small cape with grassy slopes and cliffs to the sea close by. The sample site Mace Head (MHD) is located on the western coast of Ireland in a wet and boggy area. The surrounding landscape is characterized by small hills covered with grasses and sedges, with many exposed rocks. At the southern tip of the Alaska Peninsula, near the coast, lies the measurement site Cold Bay (CBA) within a wet tundra ecosystem consisting of a variety of sedges and grasses. Finally, the station SHM is located on the island of Shemya, which belongs to a cluster of small islands southwest of Alaska.

JR-STATION

The four JR-STATION locations were installed by NIES in 2004. Three are located in the Russian taiga forest surrounded by wetlands: Demyanskoe (DEM), Karasevoe (KRS) and Noyabrsk (NOY). Additionally, one station (Igrim, IGR) was installed in a small town close to the Ob River with around 10 000 inhabitants, likewise surrounded by wetlands.

The network has been extended by five stations in the upcoming years, incorporating different biomes. Three towers have been placed in steppe regions. Azovo (AZV) and Vaganovo (VGN) are located in the immediate vicinity of highly populated cities, whereas the SVV tower (Savvushka) is installed near a small village. Additionally, one tower is located in the middle of the taiga surrounded by boreal forest (Berezorechka, BRZ), and lastly, the YAK tower was placed close to Yakutsk in the East Siberian taiga (Sasakawa et al., 2010; Belikov et al., 2019). However, not all of the JR-STATION locations are currently still in operation: the dates of the beginning and end of operation are indicated in Table S1. Since the towers are provided with two to four different sampling heights up to 85 m a.g.l. (above ground level), only the measurements from the highest inlet are used in this study. The CH₄ measurements are reported on the NIES 94 scale and have been converted to the NOAA 2004 scale following Zhou et al. (2009).

FMI–NOAA

The Finnish station Pallas (PAL) is located close to the northern edge of the Scandinavian boreal zone, with a surrounding terrain of wetlands, lakes and patches of forest (Hatakka et al., 2003; Aalto et al., 2007). PAL data are available as FMI Global Atmosphere Watch (GAW) CH₄ data from 2004 onwards at the World Data Centre for Greenhouse Gases (WDCGG). PAL data from 2017 are also available from the ICOS Carbon Portal (Hatakka and Ri, 2022). Like PAL, the site Tiksi (TIK) is operated by the FMI in cooperation with NOAA GML and is installed on the shore of the Laptev Sea on the Lena River delta (Uttal et al., 2013, 2016).

3.1.2 Data selection and observation uncertainties

In regional inversions, concentration peaks carry a large part of the information content on local to regional fluxes. However, transport can be erroneous, and simulated peaks can be shifted in time compared to observed ones, although the magnitude can be well represented. Such errors heavily penalize Bayesian inversions (e.g. Vanderbecken et al., 2023), so we decided to aggregate observations at the monthly scale although it can have an impact on the number of usable pieces of information in the system. This focuses the inversion on emission trends and seasonal cycles.

In the observation vector y⁰ (Sect. 2.1), we use the monthly averages of the available CH₄ atmospheric measurements at each site. When hourly quasi-continuous data were available, only measurements between 12:00 and 16:00 local time were selected, assuming a well-mixed boundary layer, which is better simulated by the model (Sect. 3.3). The discrete observations are not filtered by the time of day the measurement was taken. However, the data sets contain several measurement outliers, mostly strong concentration peaks related to local emissions, which are difficult to simulate with our transport model. We excluded such peaks from the observations used for the inversion if they differed by more than 5 % (or 100 ppb) from the monthly average. Depending on the measurement site, between 8 % and 20 % of the observations are discarded this way.

Due to the discontinuity of measurement availability, the size of y⁰ for 1 year varies between 217 (2008) and 384 (2018). The number of observations per year used for the inversion (and thus the size of y⁰) can be found in Table S2. All the selected observations with the corresponding daily CH₄ concentrations are shown in Fig. 3.

Figure 3 Average daily methane concentration at each station. The observation sites are sorted by latitude.

The corresponding uncertainties in the observations are specified in the diagonal error covariance matrix R, of which an ensemble of 500 set-ups is generated (Sect. 2.1).

To generate a large number of different error set-ups, the first step consists in obtaining an estimate of the uncertainty for each station $s\in [\mathrm{1},\mathrm{41}]$ and each year $j\in [\mathrm{2008},\mathrm{2019}]$ which serves as a reference point. This is done by computing the differences between the monthly mean of the measured and corresponding modelled mixing ratios (see Sect. 3.3) in absolute values:

$$\begin{array}{}\text{(5)}& {\mathbf{\Delta }}_{m,j}^s=|{\mathit{y}}_{m,j}^{\mathrm{model}}-{\mathit{y}}_{m,j}^{\mathrm{obs}}|,\end{array}$$

with $\mathit{\left\{}m\right|\mathrm{0}\le m\le \mathrm{12}\mathit{\left\}}\right]$ as 1 month of a given year j. Then, the standard deviation of the ensemble of 12 monthly differences is computed for each year:

$$\begin{array}{}\text{(6)}& {\mathit{\sigma }}_{s,j}^R=\sqrt{{\displaystyle \frac{\mathrm{1}}{\mathrm{12}}}\sum _{m=\mathrm{1}}^{\mathrm{12}}{\left({\mathrm{\Delta }}_{m,j}^s-\stackrel{\mathrm{‾}}{{\mathrm{\Delta }}_{m,j}^s}\right)}^{\mathrm{2}}},\end{array}$$

with $\stackrel{\mathrm{‾}}{{\mathrm{\Delta }}_{m,j}^s}=\frac{\mathrm{1}}{\mathrm{12}}{\sum }_{m=\mathrm{1}}^{\mathrm{12}}{\mathrm{\Delta }}_{m,j}^s$. In the few cases when only one observation is available for a given station and a given year, no standard deviation can be computed so that the single difference between the modelled mixing ratio and measurement is used directly.

The obtained errors per station and year are subsequently varied following a log-normal distribution with ${\mathit{\sigma }}_{s,j}^R$ as its mode. This error distribution is chosen to include only a few very high outliers in the ensemble. To implement a log-normal distribution, a standard deviation ${\mathit{\sigma }}_{i,\mathrm{random}}^R$ must be provided, which is constant for each element i of the ensemble $i\in [\mathrm{1},\mathrm{500}]$. Thus, the random observation error for each station s is equal for all months within 1 year; however it varies between the different years of one element i of the ensemble. To ensure that the values of the observation errors do not vary to an unrealistic extent, a minimum of 0.5 ppb and a maximum of 150 ppb are set.

Finally, the elements of the diagonal of one error covariance matrix ${\mathbf{R}}_i^j(k,k)$ for $k\in \mathit{\{}s\in [\mathrm{1},\mathrm{41}\left]\mathit{\right\}}\times \mathit{\left\{}m\right|\mathrm{0}\le m\le \mathrm{12}\mathit{\}}$ and $i\in [\mathrm{1},\mathrm{500}]$ are defined as the variances ${\left({\mathit{\sigma }}_{s,j}^{R,i}\right)}^{\mathrm{2}}$, and the non-diagonal elements are zero.

Figure 4 shows an example of the frequency distribution of the observation errors at one of the selected sites (s= INK) for the year j=2012. The mode, and therefore the reference point, of the observation error for this year and station is around 8 ppb. To give an idea about the general magnitude of the computed uncertainties, the average of ${\mathit{\sigma }}_{s,j}^R$ over all the stations s including all years j in the period of interest is around 18 ppb.

3.2 Prior emissions

3.2.1 Emission scenarios

The emissions used as prior information are based on a set of various inventories and models. The different methane sources and sinks are described in Table 1 with their respective temporal resolution in the prior information. The natural CH₄ sources include emissions from wetlands, the Arctic Ocean and geological sources. Terrestrial freshwater systems other than wetlands are hereby not taken into account as a separate source since ponds in permafrost peatlands and thermokarst lakes are usually shallow (less than 1 m depth), which falls under the usual definition for wetlands with standing water up to a depth of 2 to 2.5 m (e.g. Tiner et al., 2015).

Natural methane emissions caused by biomass burning due to wildfire events are combined with anthropogenic biofuel activities for simplification. Since emissions caused by termites are negligible in the Arctic, they are not taken into account in this study. For the CH₄ sink, soil oxidation is included as negative emissions. To reduce the number of sectors to optimize, emissions related to the exploitation and distribution of mineral oil and gas have been combined into a single data set. The same applies to the emissions from agricultural activities and waste management.

Figure 4 Frequency distribution of the 500 random observation errors $\left\{{\mathit{\sigma }}_{s,j}^{R,i}\right\}$ for $i\in [\mathrm{1},\mathrm{500}]$ at s= INK for the year j=2012. The blue line marks the mode ${\mathit{\sigma }}_{\mathrm{INK},\mathrm{2012}}^R$.

Poulter et al. (2017)Weber et al. (2019)Etiope et al. (2019)

Table 1 Methane sources and sink taken into account in the prior emissions. The share of the global emissions of each source is based on the average fluxes between 2008 and 2019. For data sets with inter-annual differences, the range between the lowest and highest emissions is given. EDGARv6 is described in Crippa et al. (2021), and GFED4.1 is described in Randerson et al. (2017).

For the natural sources as well as the soil sink, monthly climatological data sets are used for the whole period so that their total fluxes do not differ between the years 2008 and 2019. The emissions from anthropogenic sources vary between the different years covered in this study, following the EDGARv6 emission trends (Crippa et al., 2021). Emissions caused by fossil fuel activities generally increase between 2008 (15.96 Tg yr⁻¹) and 2019 (17.31 Tg yr⁻¹), though the highest annual emissions occur in the years 2014 and 2015. Methane emissions from agricultural activities and waste management also increase slightly throughout the period of interest, however just by less than 0.18 Tg yr⁻¹. The combined biomass burning scenario also shows some inter-annual variability, though without any apparent tendencies. The lowest annual emissions occur in 2009 (1.87 Tg CH₄), and the highest are in 2012 (3.99 Tg CH₄).

At the intra-annual scale, in contrast to the other natural CH₄ sources, the wetland scenario has a clear seasonality in the Arctic, with higher emissions during the summer months. According to the data set used for this study, the highest wetland emissions occur in August (10.72 Tg CH₄ per month), and the lowest are in January (0.04 Tg CH₄ per month). The soil methane oxidation has a seasonal pattern symmetric to the wetland emissions, with the maximum uptake taking place in August (−1.02 Tg CH₄ per month) and a minimum in January (−0.01 Tg CH₄ per month). The combined biomass burning scenario shows a small seasonal variability with predominantly higher emissions during the summer. Between 2010 and 2016, the highest monthly CH₄ emissions occur in July and from 2017 to 2019 the peak emissions take place in August. Hereby, the maximum of the methane emissions ranges between 0.49 Tg CH₄ per month (2009) and 1.91 Tg CH₄ per month (2017). The first 2 years within the period of interest do not fall into this seasonal pattern with increased CH₄ fluxes during the summer months. Regarding the anthropogenic methane emissions, the agricultural and waste management fluxes also show a seasonal pattern with increased emissions during the summer. According to the inventory, the emissions are highest in June (around 0.80 Tg CH₄ per month) and lowest in January and December (around 0.67 Tg CH₄ per month). The methane emissions from oil and gas exploitation and distribution are nearly constant over the course of each year, with a maximum variation of 0.1 Tg CH₄ per month.

3.2.2 Prior uncertainties

As for the observation error, the elements of the prior error matrix B are obtained from a random sampling. The covariance matrix thereby contains both the uncertainties in the prior fluxes B^S and the uncertainties in the background mixing ratios B^B. In the following, only the methodology of the random sampling of the prior errors is explained; the details on B^B are described in Sect. 3.3.3.

For each CH₄ source or sink S, the mode σ^S is set following Baray et al. (2021):

• 50 % for S of anthropogenic emissions

• 60 % for S of wetland emissions

• 100% for S of other natural sources and soil oxidation.

A random sampling following a log-normal distribution with σ^S as its mode results in an ensemble of 500 prior errors per source or sink ${\left\{{\mathit{\sigma }}_i^{B,S}\right\}}_{i\in [\mathrm{1},\mathrm{500}]}$. We expect little sensitivity of our results to these prescribed reference values as we weight every sample according to its log likelihood (see Sect. 2.2.1).

These random errors remain identical for each region r and month of year m per element i of the ensemble. Exemplarily, Fig. 5 shows the frequency distribution of the random prior errors for S of wetlands emissions of all the set-ups.

Finally, the elements of the diagonal of one error covariance matrix ${\mathbf{B}}_i^{\mathrm{S}}(k,k)$ for $k\in \left\{S\in [\mathrm{1},\mathrm{7}]\right\}\times \left\{r\in [\mathrm{1},\mathrm{121}]\right\}\times \left\{m\in [\mathrm{1},\mathrm{12}]\right\}$ are defined as the variances ${\left({\mathit{\sigma }}_i^{B,S}\right)}^{\mathrm{2}}$. Hereby, ${\mathbf{B}}_i^{\mathrm{S}}$ is identical for each year.

Figure 5 Frequency distribution of the 500 random prior errors ${\left\{{\mathit{\sigma }}_i^{B,S}\right\}}_{i\in [\mathrm{1},\mathrm{500}]}$ for S wetland emissions. The blue line marks the mode σ^wetlands.

The off-diagonal elements ${\left({\mathit{\sigma }}_{m,n}\right)}^{\mathrm{2}}$ with m the row and n the column of the corresponding matrix ${\mathbf{B}}_i^{\mathrm{S}}$ are determined by applying spatial and temporal correlations. ${\mathbf{B}}_i^{\mathrm{S}}$ is hereby a symmetrical matrix so that ${\left({\mathit{\sigma }}_{m,n}\right)}^{\mathrm{2}}$ is identical to ${\left({\mathit{\sigma }}_{n,m}\right)}^{\mathrm{2}}$.

The off-diagonal errors are computed as follows:

$$\begin{array}{}\text{(7)}& \begin{array}{rl}{\left({\mathit{\sigma }}_{m,n}\right)}^{\mathrm{2}}& ={\left({\mathit{\sigma }}_{n,m}\right)}^{\mathrm{2}}=\left({\mathit{\sigma }}_{m,m}\times {\mathit{\sigma }}_{n,n}\right)\times \exp \left(-{\displaystyle \frac{\mathrm{\Delta }t}{t_{\mathrm{corr}}}}\right)\\ & \times \exp \left(-{\displaystyle \frac{\mathrm{\Delta }d}{d_{\mathrm{corr}}}}\right),\end{array}\end{array}$$

with Δt as the temporal difference between the rows/columns m and n and Δd as the spatial difference referring to the centres of the corresponding regions. For the spatial correlation d_corr a distance of 500 km is used, and the temporal correlation t_corr has a fixed value of 1 week.

3.3 Modelled CH₄ mixing ratios

As mentioned in Sect. 2.1, the simulated equivalents to the observations are included in the observation operator H. In this case, the elements of H consist of the monthly CH₄ mixing ratios sectioned into sub-regions and sectors as well as the monthly averages of the background mixing ratios by station. H is hereby linear since only emissions and transport of CH₄ are taken into account. The oxidation of methane by hydroxyl radicals (OH) is neglected inside the domain of simulation since the lifetime of CH₄ is ≈9 years (Prather et al., 2012) and the air masses remain in the domain up to 2 months (Berchet et al., 2020). The methane sink from the hydroxyl radical is however accounted for in the global simulations used to compute the background mixing ratios (Sect. 3.3.3).

3.3.1 Transport model set-up

The modelled CH₄ mixing ratios were obtained by using the Lagrangian atmospheric transport model FLEXPART (FLEXible PARTicle) version 10.3 (Stohl et al., 2005; Pisso et al., 2019). This model simulates numerous trajectories of infinitesimally small air parcels, called particles, and can be used either forward or backward in time. FLEXPART is an offline model that is driven by meteorological data from the European Centre for Medium-range Weather Forecast (ECMWF) ERA5 (Hersbach et al., 2020) with 3-hourly intervals and 60 vertical layers. ECMWF data are retrieved and formatted using the flex_extract toolbox (Hittmeir et al., 2018; Tipka et al., 2020). In this study, 2000 particles are released at each observation site and timestamp (receptor) and followed 10 d backwards in time. The horizontal resolution is $\mathrm{1}{}^{\circ }\times \mathrm{1}{}^{\circ }$, which is quite commonly used for inverse modelling set-ups using Lagrangian particle dispersion models in high northern latitudes (e.g. Thompson et al., 2017; Ishizawa et al., 2019).

3.3.2 Source contribution

By sampling the near-surface residence time of the various backward trajectories of the particles, the source-receptor sensitivity matrices, also called footprints, of each observation site can subsequently be determined. The near-surface residence time hereby corresponds to particles below 500 m instead of, for instance, the planetary boundary height (PBL) since simulated PBL in the Arctic can be unrealistically small, especially during the winter months. The thus obtained footprints define the connection between the fluxes discretized in space and time and the change in concentrations at the receptor (Seibert and Frank, 2004). To finally obtain a time series of modelled CH₄ mixing ratios, a time series of footprints is integrated with discretized methane emission estimates. Here, monthly averages of the footprints of each receptor are used to determine the mixing ratios for each sector (see Table 1 in Sect. 3.2.1) and sub-region (see Fig. 1b in Sect. 2.3).

The magnitude of the thus obtained total CH₄ mixing ratios, including all methane sources and the soil sink, ranges roughly between 3 and 90 ppb depending on the month of the year and location of the observation site, and the average standard deviation is around 14 ppb.

3.3.3 Background mixing ratios and uncertainties

Since CH₄ has a much longer lifetime than the released virtual particles, the previously obtained concentrations only display short-term fluctuations at the receptors. Therefore, in order to obtain a direct comparison to the measurements, the background mixing ratio needs to be taken into account.

The background mixing ratios are calculated by combining a CH₄ concentration field as the initial condition with the FLEXPART backward simulations nudged to the observations of the corresponding site (e.g. Thompson and Stohl, 2014; Pisso et al., 2019). The background thus obtained represents the average of the mixing ratios in the grid cells where each particle trajectory terminated 10 d before the observation. The initial concentration field is provided by the Copernicus Atmospheric Monitoring Service (CAMS): a CH₄ mixing-ratio field from CAMS global reanalysis EAC4 (ECMWF Atmospheric Composition Reanalysis 4) with 60 vertical layers and a 3 h temporal and a $\mathrm{0.75}{}^{\circ }\times \mathrm{0.75}{}^{\circ }$ spatial resolution has been used (Inness et al., 2019). The practical computation of background time series from background fields and backward trajectories is carried out using the Community Inversion Framework (CIF;  Berchet et al., 2021b).

The thus computed background mixing ratios show a gradual increase over the period of interest, with mean annual concentrations over all sites ranging between 1842 ppb (2008) and 1974 ppb (2019). At the intra-annual scale, the monthly background mixing ratios vary from the corresponding annual average by around 8 %. Figure S5 shows the average background mixing ratios at each station as well as their average standard deviation.

As stated previously, the background is the major share of the total modelled mixing ratios and, in this study, makes up approximately 97.6 % at continental observation sites and 99.5 % at stations located remotely. A summary of the proportion of source contribution and background mixing ratios for each station can be found in Table S4.

As mentioned before (e.g. Sect. 2.1), the uncertainties in the background mixing ratios B^B are included in the error covariance matrix B. In contrast to the uncertainties in the prior emissions B^S, which are given by region, month and CH₄ source/sink, the uncertainties in the background mixing ratios are given by observation site and month. Therefore, the size of B^B is equivalent to the number of available observations per year.

The elements of B^B are composed in a manner similar to the elements of R (Sect. 3.1.2), by first computing a reference error for each station and year and varying these values randomly to obtain and ensemble of 500 set-ups.

In this case, the standard deviations of the monthly background mixing ratios ${\mathit{y}}_{s,m,j}^{\mathrm{back}}$ per station $s\in [\mathrm{1},\mathrm{41}]$ and year $j\in [\mathrm{2008},\mathrm{2019}]$ serve as reference errors:

$$\begin{array}{}\text{(8)}& {\mathit{\sigma }}_{s,j}^B=\sqrt{{\displaystyle \frac{\mathrm{1}}{\mathrm{12}}}\sum _{m=\mathrm{1}}^{\mathrm{12}}{\left(y_{s,m,j}^{\mathrm{back}}-\stackrel{\mathrm{‾}}{y_{s,m,j}^{\mathrm{back}}}\right)}^{\mathrm{2}}},\end{array}$$

with $\stackrel{\mathrm{‾}}{y_{s,m,j}^{\mathrm{back}}}=\frac{\mathrm{1}}{\mathrm{12}}{\sum }_{m=\mathrm{1}}^{\mathrm{12}}{y}_{s,m,j}^{\mathrm{back}}$ and $m\in [\mathrm{1},\mathrm{12}]$.

Subsequently, the computed errors per station are varied following a log-normal distribution with a mode of ${\mathit{\sigma }}_{s,j}^B$. Again, in order to achieve a log-normal distribution, a random standard deviation ${\mathit{\sigma }}_{i,\mathrm{random}}^B$ must be set which is consistent per element $i\in [\mathrm{1},\mathrm{500}]$ of the ensemble. Similar to the observation errors, this means that each observation site s has identical values of background errors for every month m within 1 year, but each station may have unequal errors for the different years j of one element i of the ensemble. The lower and upper limits of the background mixing-ratio uncertainties are, hereby, 0.5 and 150 ppb.

The diagonal elements of one error covariance matrix ${\mathbf{B}}_i^{B,j}(k,k)$ for $k\in \mathit{\{}s\in [\mathrm{1},\mathrm{41}\left]\mathit{\right\}}\times \mathit{\left\{}m\right|\mathrm{0}\le m\le \mathrm{12}\mathit{\}}$ and $i\in [\mathrm{1},\mathrm{500}]$ are finally defined as the variances ${\left({\mathit{\sigma }}_{s,j}^{B,i}\right)}^{\mathrm{2}}$.

Other than the observation error covariance matrix R, B^B is not a diagonal matrix, and the non-diagonal elements are defined by applying correlations in space and time. The computation of the non-diagonal errors ${\left({\mathit{\sigma }}_{m,n}\right)}^{\mathrm{2}}$ with m as the corresponding row and n as the corresponding column of the symmetrical matrix ${\mathbf{B}}_i^{\mathrm{B}}$ is similar to the implementation of correlations for the prior error covariance matrices ${\mathbf{B}}_i^{\mathrm{S}}$ (Sect. 3.2.2):

$$\begin{array}{}\text{(9)}& \begin{array}{rl}{\left({\mathit{\sigma }}_{m,n}\right)}^{\mathrm{2}}& ={\left({\mathit{\sigma }}_{n,m}\right)}^{\mathrm{2}}=\left({\mathit{\sigma }}_{m,m}\times {\mathit{\sigma }}_{n,n}\right)\times \exp \left(-{\displaystyle \frac{\mathrm{\Delta }t}{t_{\mathrm{corr}}}}\right)\\ & \times \exp \left(-{\displaystyle \frac{\mathrm{\Delta }d}{d_{\mathrm{corr}}}}\right),\end{array}\end{array}$$

with Δt as the temporal difference between the rows/columns m and n and Δd as the spatial difference referring distance between the two corresponding measurement sites. The correlation lengths are d_corr=500 km for spatial correlations and t_corr= 1 week for temporal correlations.

Drawing Monte Carlo members by varying d_corr and t_corr is perfectly possible and would impact the number of degrees of freedom relative to the background (the higher d_corr and t_corr, the lower the number of degrees of freedom). In the present work, we limit our analysis to a single set of d_corr and t_corr. The chosen configuration is a compromise to account for the consistent influence of the background between neighbouring stations and successive time steps while avoiding forcing unrealistic isotropic correlations when close sites are influenced by different background air masses. A future work will include the four-dimensional background fields in the control vector as the background.

4 Results

4.1 Performance of the inversions in the observation space

To evaluate the performance of the inversion, the prior and posterior CH₄ mixing ratios are compared to the observations. Figure 6 shows the Taylor diagrams indicating the Pearson correlation coefficient to determine similarities between the observations and simulations as well as the normalized standard deviation (SD) displaying how well the variability in the modelled mixing ratios is captured. Thus, a shorter distance to the reference point indicates a closer fit to the measured mixing ratios. In Fig. 6, we split the results for the full data set and de-trended data. The performance of the simulations for the full data set is mostly driven by the long-term trend. The de-trended data exhibit the performance in terms of seasonal cycle.

In general, and as expected, the posterior results show better agreement with the observations compared to the prior mixing ratios of the corresponding observation site. This is more distinctive for the trended (Fig. 6a–c) than for the de-trended time series (Fig. 6d–f), although in both cases the majority of the posterior mixing ratios are closer to the measurements than the prior ones. This confirms that the climatological priors are not good enough as year-to-year changes were present and the inversion can realistically improve the flux trends. Both the normalized standard deviation and the correlation coefficient should ideally be close to 1. The prior trended SD range between 0.19 and 1.62, and the correlation coefficients are between 0.20 and 1.0. For the posterior results the values lie between 0.19 and 1.00 (standard deviation) and 0.29 and 1.0 (correlation coefficient). Regarding the de-trended time series, the normalized SD lies between 0.19 and 2.61 (prior) and 0.02 and 0.99 (posterior), and the correlation coefficient ranges between 0.20 and 1.41 (prior) and 0.10 and 1.00 (posterior).

Figure 6 Examples of Taylor diagrams for various site categories. Raw (a, b, c) and de-trended (d, e, f) mixing ratios over the whole period of interest. Prior simulated mixing ratios are indicated with circles; posterior ones are indicated with diamonds.

The improvement in the posterior results is quite evident for observation sites which are remote from methane emission sources, such as ALT or ZEP (Fig. 6a), where the posterior results are nearly equal to the observations. Here, the standard deviations have a maximum deviation of 0.10 from the observation point, whereas the difference between the correlations is ≤0.02. It is however noteworthy that the prior CH₄ concentrations show already a good agreement with the observations at those remote stations, which are often referred to as background observation sites. This good fit can be explained by the fact that background mixing ratios are computed using global mixing-ratio fields generated by systems optimized using these same remote sites.

A much larger improvement can be observed at sites close to the North American coast, such as INK, BLK and CHU (Fig. 6b). In general, the majority (8 out of 11) of the measurements from the coastal stations have a lower standard deviation than their modelled equivalents, which implies that the variability in the modelled mixing ratios is overestimated. The magnitude of the prior modelled CH₄ mixing ratios is overall higher than the measurements at North American observation sites in high northern latitudes (up to approximately 80 ppb and on average 50 ppb).

The simulations of continental observation sites further south in North America (e.g. BRA) as well as continental sites in Russia such as NOY and IGR (Fig. 6c) show, in general (8 out of 12 sites), a normalized standard deviation which is lower than the observations. This indicates an underestimation of the variability in the simulated CH₄ mixing ratios, both in the prior as well as in the posterior results. However, the correlation with the observations could still be improved by the inversion.

The only observation site where the posterior results show less agreement with the corresponding measurements is ABT (Fig. 6b), where, in contrast to most other sites in North America, the observations are significantly higher than the simulated CH₄ mixing ratios by up to approximately 100 ppb. Local fluxes (mostly from urban environments) and complex topography (mountain range surrounding the flat area around Vancouver, Canada) are likely to influence the observations at this site and are poorly represented by the model at a coarser resolution. A higher transport resolution and finer-scale inversion regions could solve this issue in a future study; stations too close to urban centres and in too complex topographical configurations could also be discarded altogether to pan-Arctic studies focusing on large-scale patterns.

4.2 Distribution of information in the inversion system

4.2.1 Impact of observations on the inversion system

To further analyse the network efficiency, the sensitivity matrices HK (see Sect. 2.2.2) are calculated for each year and averaged over the whole period of interest (Fig. 7). The percentages indicate how much of the theoretically available observations at each site are actually used by the inversion. The observation sites, which are located remotely from any other station, mostly along the Arctic, Atlantic and Pacific shores, show values of almost 100 %, which means that the information provided by the measurements is almost entirely used, mostly to constrain background concentrations. This is confirmed by the amplitude of the background at these sites, as shown in Fig. S5, where the ratio between the standard deviation of the simulated signal from the background and from emissions has been computed and shows similar patterns of the sensitivity to observations. In areas where the observation network is much denser (e.g. in the southeast of Canada and, to a lesser extent, in the Siberian lowlands), most observations contribute to less than 50 % to the inversion. Lower constraints in dense continental areas are caused by redundant constraints from neighbouring sites in the same emission areas and/or higher noise due to transport errors from nearby emissions. The latter has the largest impact if the site is located close to CH₄ emission sources.

Figure 7 Sensitivity of the inversion to observation sites as derived from by the sensitivity matrix HK. Larger and darker circles indicate a higher usage of available observations [%] by the inversion. The percentage thereby shows the share of used and theoretically available measurement data.

4.2.2 Noise and information content in the inversions

The trace of the influence matrix tr(KH) (equal to the trace of the sensitivity matrix tr(HK)) indicates how much noise is contained in the provided observations and how the information content is used by the inversion (Sect. 2.2.2). The closer the value of tr(KH) is to the number of available observations, the more useful each given observation is for the inversion. Furthermore, the ratio between the number of observations used to constrain the emissions and those used to constrain the background mixing ratios can be determined by separately calculating tr(KH_emis) and tr(KH_back), using only the corresponding elements of KH. The obtained traces for each year are given in Table S2, and Fig. 8 shows the ratio between tr(KH_emis) and tr(KH_back).

In total, tr(KH) ranges between approximately 60 % and 75 % of the number of available observations, with the majority constraining the background mixing ratios. Only around 10 % of the available observations are used for constraining the emissions, whereby the share remains relatively constant through the years. Moreover, it is noticeable that the trace of KH is closer to the number of observations during the years in which the smallest numbers of measurements are provided (e.g. 2008 and 2019).

Figure 8 Traces of influence matrices divided by the number of available measurements of the corresponding year. See Sect. 2.1 for details on the computation of the influence matrix. The closer tr(KH) is to 1, the more observations are used in the inversion.

With this limited availability of data, a higher percentage of the observations is used as information for the inversion. By contrast, in years during which more observations are available (e.g. 2015), a higher share is identified as noise and hence redundant information, similarly to spatial redundancy in regions where the observation network is denser.

Figure 9 Seasonal variation in tr(KH) averaged over the period of interest (2008–2019). The monthly traces are divided by the number of available observations for the corresponding month.

The fraction of useful information in the available observations follows a seasonal variability as shown in Fig. 9 (see Fig. S3 for seasonal variations in the individual years).

The constraints on the emissions during winter are relatively small not only since the CH₄ emissions are comparatively smaller than during the summer months but also because meteorological conditions (in particular a stratified cold boundary layer) make the comparison of observations with simulations more challenging. During the summer months, a higher fraction of observations (up to 20 %) is used to constrain emissions. In general, the total trace tr(KH) is higher during the summer months, which means that fewer of the observations are identified as noise. However, additional constraints on the emissions during summer do not ensure constant constraints on the background. Instead, a share of the constraints on the background mixing ratio is transferred to constrain the emissions during the summer months.

By construction tr(KH_back) is proportional to B_back and H_back. H_back cannot be reduced due to the physics of the atmospheric transport (see Table S4). One way to reduce the share of the information constraining the background in the inversion set-up would be to decrease the uncertainties in the background mixing ratios in B_back. This relies, however, also on the performances of simulations of global CH₄ concentration fields. Even though in recent years those applications have already improved, they still do not provide a sufficient level of precision that would allow for reducing the uncertainties for the implementation of our inversion set-up (Inness et al., 2019).

Moreover, the limited transport backwards in time in FLEXPART (10 d in our case) is much smaller than the average residence time of air masses in the Arctic (typically a few weeks; see e.g. Berchet et al., 2020). Hence, part of the influence of Arctic fluxes on observations is diluted in the background in our system. One way of mitigating this issue would be to dramatically increase the backward transport time of virtual particles up to a few weeks; but to limit numerical artefacts, multi-week backward simulations need a very large number of particles to be accurate, at the expense of much higher computational costs. Another way of solving the issue would be to fully couple FLEXPART within a global circulation model, thus accounting for the influence of fluxes on observations indefinitely backwards in time; this is what is done in e.g. Maksyutov et al. (2021) or could be done in the Community Inversion Framework with one of the available global models (LMDZ, Laboratoire de Météorologie Dynamique Zoom model, or TM5, Transport Model; Berchet et al., 2021b).

Another option to increase the ratio of information used to constrain emissions instead of the background would be to use higher temporal resolution for the observations. A compromise would be to use daily afternoon averages instead of monthly averages. However, as illustrated in Belikov et al. (2019), Arctic observations, even at the daily scale, can see strong daily peaks due to emissions in combination with meteorological conditions. Such peaks may be poorly reproduced by the model and could be shifted in time and magnitude, posing a double-penalty effect to the inversion (Vanderbecken et al., 2023). This can have a critical impact on the inversion conclusion, especially with diagonal observation covariance matrices R, as done in the present work, consistent with the general practice of the inversion community.

4.2.3 Spatial distribution of constraints on regions and sources

The influence matrix KH defines how well each emission sector is constrained by the inversion in each sub-region. The majority of the CH₄ sources are quite poorly constrained in the sub-regions defined in Sect. 2.3, with the elements of the influence matrix being less than 10 %. In comparison to that, the wetland emissions are relatively well constrained, as shown in Fig. 10. Hereby, the figure on the left shows the average constraints over all years, while the middle and right figures show 2 exemplary years (2011 and 2014) to highlight inter-annual differences. The remaining years are shown in Fig. S1 in the Supplement.

Figure 10 Regional constraints on wetland emissions as derived from the influence matrix KH. Darker areas thereby indicate higher constraints. The percentages of the areas refer to the corresponding summed elements of KH. The observation sites are marked as grey circles.

The average values of the annual influence matrices (Fig. 10a) indicate that the current observation network is able to constrain wetland emissions well for most North American sub-regions. In Eurasia, on the other hand, most areas are unseen by the inversion and the well-constrained areas are predominantly limited to certain parts of Siberia (e.g. the West Siberian Plain). This is partly due to the distribution of the observation network (the denser the network, the better the constraints) and to the heterogeneity of data collection within the period of interest (some years have many more available observations than others, especially towards the end of the period). As shown in Fig. 10b and c, the extent of the constraints strongly varies between the different years due to the availability of observations in Eurasia. Those variations can also be noticed in North America; however, the well-constrained areas remain relatively identical over the whole period.

Another cause of the limited constraints on the emissions is that the available observations in Russia are rather used to constrain the background mixing ratios (see Sect. 4.2.2). In North America, where a larger number of observation sites are established and more evenly distributed over the area, the observations of certain stations are used to provide the information on the background.

Installing additional observation sites in high northern latitudes in Eurasia would therefore be useful to better constrain local emissions in the future. However, measurement stations in lower latitudes at the sub-Arctic boundary would also be necessary to better constrain transport from CH₄ hotspots such as China, India and the Middle East.

4.3 Analysis of posterior fluxes

4.3.1 Total methane fluxes

In order to compare the prior to the posterior fluxes, the area of interest is divided into four different supra-regions, North America, East Eurasia, West Eurasia and the Arctic (including the High and Low Arctic), as shown in Fig. 11.

Figure 11 Supra-regions for analysis of posterior CH₄ fluxes.

Since most emission sources do not show large differences between the prior and the posterior state and are also poorly constrained by the inversion (Sect. 4.2.3), the sectors described in Sect. 3.2.1 are combined to wetlands and other natural (including the CH₄ sink from soil oxidation) and anthropogenic emissions. In particular, geological fluxes from the ocean do not deviate significantly from the prior and are not further commented on here. Thereby, the combined natural and anthropogenic fluxes from biomass burning are included in the natural emission sources for simplification since the natural emissions exceed the anthropogenic ones.

The mean annual prior and posterior CH₄ emissions in each region are shown in Fig. 12 and, in more detail, in Table S3 for the set-up ${\mathit{x}}_{\max }^{\mathrm{a}}$ with the highest log likelihood (Sect. 2.2.1) together with the corresponding uncertainties obtained from the P^a matrix. As expected, the poorly constrained anthropogenic and other natural emissions do not show significant changes between the prior and posterior fluxes for either of the regions, neither in their magnitude nor in their uncertainties. The wetland emissions are decreased in the posterior state, except in West Eurasia. The largest decrease is found in North America, which is also the region best constrained by the inversion. Here, the prior wetland emissions have a magnitude of around 30±26 Tg CH₄ yr⁻¹, whereas the posterior emissions amount to 19±13 Tg CH₄ yr⁻¹. Even though the uncertainties in the posterior wetland fluxes are still high, at around 69 %, they are reduced by around 17 % in comparison to the prior uncertainties. In East Eurasia, the wetland emissions are decreased from approximately 14±12 to 12±10 Tg CH₄ yr⁻¹ and in the Arctic from 13±11 to 10±8 Tg CH₄ yr⁻¹ with an uncertainty reduction of 8 % and 6 %, respectively.

Figure 12 Mean annual CH₄ emissions by sector for ${\mathit{x}}_{\max }^{\mathrm{a}}$ with corresponding uncertainties.

Comparison to global inversion set-ups

In order to compare this study to other inversion set-ups, the prior and posterior emissions are set against five different posterior states obtained with variational inversion frameworks used for the Global Carbon Project (GCP). The comparative CH₄ fluxes are hereby an updated version of the results from Saunois et al. (2017, 2020a). The variational inversions are performed globally with two different inversion systems, CIF–LMDZ using surface observations (Thanwerdas et al., 2021) and PYVAR–LMDZ (PYthon VARiational) using satellite observations from GOSAT (Greenhouse Gases Observing Satellite) (Zheng et al., 2018). Inversion set-ups 1 and 2 use the prior fluxes distributed for the Global Methane Budget and Atmospheric Tracer Transport Model Intercomparison Project (TransCom) chemical fields, with the latter including OH inter-annual variability from Patra et al. (2021). The third set-up is a sensitivity test where freshwater fluxes are added in the prior state. The mean annual total CH₄ emissions in the different regions are shown in Fig. 13. Since the GOSAT observations are not available for the years 2008 and 2009, the PYVAR–LMDZ posterior results are averaged over the remaining period of interest.

Figure 13 Total mean annual CH₄ emissions in comparison to different inversion set-ups from the GCP.

In general, the total fluxes of the variational inversion set-ups are all lower than the posterior results of this work. The largest discrepancies are found in the Arctic, where the total posterior fluxes are up to 59 % higher than the results from the GCP and only inversion set-up 3 lies within the posterior uncertainty range of our inversion set-up. The CH₄ fluxes of the variational inversion set-ups are lower by between 14 % and 44 % in North America, 38 % and 51 % in East Eurasia, and 18 % and 38 % in West Eurasia in comparison to our posterior emissions. In all of the regions, the results from the inversions using satellite data (PYVAR–LMDZ) are the least consistent with the posterior CH₄ emissions obtained in this work. The smallest difference to our results is given by the inversion set-up in which the freshwater emissions are added in the prior state (set-up 3).

As our system explicitly provides posterior uncertainties, contrary to many other inversion systems, it is possible to assess the consistency of our results with other inversions. The discrepancies between the posterior methane emissions from our study and the global variational inversions could be due to the fact that global inverse systems do not perform as well in high latitudes. This has already been identified in Saunois et al. (2017) and can be tracked back to the following.

• i.

  Global inversions use fewer observation sites in the Arctic.

• ii.

  Global inversions constrained by satellite measurements (GOSAT; Infrared Atmospheric Sounding Interferometer, IASI) provide fewer data points above 30^∘ N compared to regions in lower latitudes because of factors such as the solar zenith angle, the surface albedo and the limited light during the polar night. More recent inversions have started using TROPOspheric Monitoring Instrument (TROPOMI) data (Tsuruta et al., 2023), with possibly a higher density of data points at high latitudes, but these data were not publicly available at the time of the present work.

• iii.

  Global models with very low resolution cannot reproduce the Arctic atmosphere properly.

However, the discrepancies should be further investigated.

Comparison to previous Arctic studies

In comparison to previous studies using inverse modelling to assess methane emissions in high-northern-latitude regions, our results lie roughly in the same magnitude. Thompson et al. (2017) determined the total CH₄ emissions between 2005 and 2013 to lie between 16.6 and 17.1 Tg yr⁻¹ in North America (above 50^∘ N), and Baray et al. (2021) estimated the combined natural and anthropogenic emissions in Canada at 16.6 and 18.2 Tg yr⁻¹ (between 2010 and 2015). Both values are within the lower limit of the uncertainty range of our ensemble of posterior states in North America (31±15 Tg yr⁻¹). Berchet et al. (2015b) estimated the methane fluxes in the Siberian lowlands to be between 5 and 28 Tg yr⁻¹ in the year 2010 (comparable to region East Eurasia in this study at 34±18 Tg yr⁻¹). In Eurasia, the total CH₄ emissions obtained by Thompson et al. (2017) are between 55.2 and 59.5 Tg yr⁻¹, which is at the upper limit of the uncertainty range of the results from our study for the combined areas of East Eurasia and West Eurasia (43±23 Tg yr⁻¹).

Due to the differences in the spatial extent of the regions covered in those studies, it is, however, difficult to obtain reliable comparisons of the estimated methane emissions.

4.3.2 Trends of emission sources

In a changing climate, detecting changes in trends of regional emissions in high northern latitudes is critical. Therefore, the trends of all 5000 possible posterior fluxes from the ensemble (see Sect. 2.1) have been calculated by sector and region. The results for wetland emissions, which is the only source well constrained by the inversion, are shown in Fig. 14 for North America and East Eurasia.

• The mean annual CH₄ emissions are displayed on the horizontal axis, and the corresponding trend of the annual wetland fluxes is on the vertical axis.

• The associated probability density functions (PDFs) are shown next to the corresponding axes.

• The darker-shaded segments show the range of the ensemble $\mathit{\left\{}{\mathit{x}}_{\max }^{\mathrm{a}}\mathit{\right\}}$ with the most plausible error configurations (Sect. 2.2.1), which make out 55 % of the total ensemble.

• The posterior result with the maximum log likelihood ${\mathit{x}}_{\max }^{\mathrm{a}}$ is highlighted as well as the trend and the mean annual emissions of the prior flux estimates.

Since the data set of wetland emissions is equal for each year within the period of interest, there is no trend in the prior state. The trend of the posterior wetland emissions in North America (Fig. 14a), including all possible uncertainty configurations, ranges approximately between −7.3 and 12.2 % yr⁻¹ with corresponding mean annual emission between around 15 and 30 Tg CH₄ yr⁻¹. The trends of the corresponding ensemble of $\mathit{\left\{}{\mathit{x}}_{\max }^{\mathrm{a}}\mathit{\right\}}$ range between −1.4 and 1.2 % yr⁻¹, with 65 % of the 2740 posterior results showing a negative trend. The most plausible of all set-ups x_max, according to the log likelihood, also has a decreasing trend of −1.4 % yr⁻¹. Thus, according to our system, although small (less than 20 % per decade), there is a plausible negative (although uncertain) trend on wetland emissions in North America between the years 2008 and 2019.

The trend of the posterior results of the wetland emissions in East Eurasia shown in Fig. 14b ranges between −7.5 and 11.7 % yr⁻¹, and the mean annual amount of CH₄ emissions ranges between 10 and 15 Tg CH₄ yr⁻¹. Here, the elements of $\mathit{\left\{}{\mathit{x}}_{\max }^{\mathrm{a}}\mathit{\right\}}$ do not include any negative trends, with values ranging between 0 and 2.1 % yr⁻¹, and ${\mathit{x}}_{\max }^{\mathrm{a}}$ shows increasing trend of 0.8 % yr⁻¹. The results point to a very small but statistically significant positive trend in East Eurasia. A positive growth rate in CH₄ mixing ratios between 2009 and 2019 in West Siberia was detected by Someya et al. (2020) and attributed to increased wetland emissions in this area, which is compatible with our conclusion.

Figure 14 Trend and mean annual fluxes of wetland emissions for the ensemble of posterior results with corresponding density distributions. Brighter colours of scatters indicate a higher density.

Figure 15 Trend and mean annual fluxes of biomass burning emissions for the ensemble of posterior results with corresponding density distributions.

To give an example of a CH₄ emission source with a trend in the prior state, Fig. 15 shows the emissions from biomass burning of the two previously discussed regions. Since the uncertainties in the emissions from biomass burning have been chosen to be higher in comparison to the wetland emissions (Sect. 3.2), the posterior results contain several negative CH₄ fluxes, which are not included in the figures. In both regions, the prior state shows an increasing trend of 30 % yr⁻¹ in North America and 9 % yr⁻¹ in East Eurasia, with corresponding mean annual emissions of 1.15 and 2.01 Tg CH₄ yr⁻¹.

In North America (Fig. 15a), the posterior trends show a large variability, with values of the total ensemble lying between −11 and 92 % yr⁻¹. The trend of the most plausible configuration ${\mathit{x}}_{\max }^{\mathrm{a}}$ is around 21 % yr⁻¹, which still indicates a positive trend, though one that is lower than the prior estimate. The magnitude of the CH₄ emission from biomass burning in North America is predominantly decreased in the posterior state, with ${\mathit{x}}_{\max }^{\mathrm{a}}$ being 0.7 Tg yr⁻¹ lower than the prior mean annual emissions. In East Eurasia (Fig. 15b), the mean annual emissions of the total ensemble lie between 0.35 and 3.47 Tg CH₄ yr⁻¹, and the trends show, similarly to North America, high variability in the posterior state (between −8 and 89 % yr⁻¹). The ensemble $\mathit{\left\{}{\mathit{x}}_{\max }^{\mathrm{a}}\mathit{\right\}}$ on the other hand only shows minor deviations, and both the trend and the fluxes are close to the prior state.

It has to be mentioned that, even though the fluxes from biomass burning are partially well constrained in some regions and years, the emissions are poorly constrained throughout the whole period of interest, and as described in Sect. 4.2.3, the results are highly uncertain.

For most of the methane emission sources and sinks, the results of the posterior ensemble do not show major deviations from the prior state at all, independent of the prior and observation error.

In addition to that, the obtained trends of the posterior fluxes could as well be influenced by the varying data availability, and the results are therefore still highly uncertain.

Figure 16 Normalized total CH₄ fluxes of the prior (dash-dotted lines) and posterior (continuous lines) state per month in North America. The coloured lines display the different years; the black lines show the average over all years.

Figure 17 Normalized wetland CH₄ fluxes of the posterior state per month in North America. The dash-dotted black line shows the average prior; the continuous line shows the average posterior state.

4.3.3 Seasonal variability

Subsequently, the seasonal cycles of the prior and posterior CH₄ fluxes are examined. Figure 16 shows the seasonal cycle of the prior and posterior fluxes for the total CH₄ emissions in North America. Hereby, the displayed posterior fluxes are the median values of the ensemble $\mathit{\left\{}{\mathit{x}}_{\max }^{\mathrm{a}}\mathit{\right\}}$. To achieve a better comparison between the different years, the monthly values are divided by the maximum methane fluxes of the prior state.

Figure 18 Total seasonal prior (purple) and posterior (blue) CH₄ emissions between 2008 and 2019.

Over the course of the period of interest, the prior emissions show greater consistency in the annual seasonal cycles. The peak of the emissions is predominantly in August, though sporadically already in June. The sectors that contribute to the seasonal differences in the prior state are emissions from wetlands and biomass burning as well as the soil oxidation (Sect. 3.2.1). Since the data set of the wetland emissions and the soil oxidation are consistent for each year, the differences in the seasonality of the prior fluxes is entirely driven by the CH₄ fluxes from biomass burning. In comparison to that, the average of the posterior state still reaches the maximum emissions in August; however the peak is less pronounced, and the emissions decrease more gradually during the autumn months. The annual seasonal cycles of the posterior fluxes are more divergent from each other. The majority of the years still show the highest methane emissions in August, although some years (e.g. 2012 and 2015) have a local minimum during that month. Unlike the prior state, the differences in the seasonal fluxes of the different years are not exclusively influenced by emissions from biomass burning. As shown in Fig. 17, most of the changes in the seasonal cycle of the total CH₄ emissions arise from adjustments of the monthly wetland fluxes since the local high and low points are predominantly during the same time of the year.

The other three regions assessed in this study show similar minor changes in the seasonal cycles, which are predominantly influenced by wetland emissions and do not show any clear seasonal pattern between the different years.

Figure 19 Seasonal CH₄ emissions from biomass burning between 2008 and 2019.

4.3.4 Inter-annual variability

Finally, we inquire into the inter-annual differences in the methane emissions. Figure 18 shows a time series of the total CH₄ fluxes in the prior and posterior state. The posterior emissions are again obtained by calculating the median of the most plausible results $\mathit{\left\{}{\mathit{x}}_{\max }^{\mathrm{a}}\mathit{\right\}}$, with the corresponding minimum and maximum values of this ensemble marking the uncertainty range. Therefore, in the following section, the given quantification of the monthly CH₄ emissions refers to the median values.

In West Eurasia (Fig. 18c) there is effectively no inter-annual variability, in neither the prior nor the posterior state, with peak values of 1.2 Tg CH₄ per month each year; a small deviation from the prior state comes from very small constraints on that region. Similarly, the emissions in the Arctic (Fig. 18d) show very few deviations in the prior state, with the maximum fluxes only deviating by 0.1 Tg. The posterior fluxes, on the other hand, gradually decrease from a peak value of 4.0 Tg CH₄ per month in 2008 to 3.1 Tg CH₄ per month in 2015, which marks the year with the largest difference from the prior state. In the following years, the emissions slightly increase again to a maximum value of 3.3 Tg per month and stay constant until the end of the period of interest. In East Eurasia, shown in Fig. 18b, the inter-annual variations are fairly low both in the prior and posterior state. According to the prior estimates, the highest methane emissions occur in 2012 with a peak value of 5.6 Tg, CH₄ per month. Also, the fluxes are slightly lower at the beginning of the period of interest (around 4.9 Tg) than during the last 4 years (around 5.4 Tg). The posterior emissions mostly follow that pattern with a maximum reduction of 0.5 Tg in 2012. The most prominent variabilities are, unsurprisingly, found in North America (Fig. 18a), which is best constrained by the inversion. As mentioned before, the inter-annual differences in the prior state are caused by variations in the emissions from biomass burning. The largest prior total methane emissions occur in 2017 with a maximum of 9.5 Tg CH₄ per month; the year with the lowest peak emissions (7.9 Tg CH₄ per month) is 2008. Like the emissions in East Eurasia, the peak values of the prior state are generally higher during the last years of the period of interest, whereas the years with the lowest CH₄ emissions can be found at the beginning. The posterior fluxes in North America vary between peak emissions of 5.0 Tg (2015) and 6.6 Tg (2018) without showing a clear pattern. In fact, the year 2008, which shows the lowest emissions in the prior state, has the second highest methane fluxes of all years, whereas the emissions in the year 2019 are almost as low as in 2015.

An explanation for the large discrepancies in the inter-annual variabilities between the prior and the posterior state in North America is the reduction in fluxes from biomass burning in the posterior state which are shown in Fig. 19. The prior estimates show a large variability with exceptionally high emissions in the years 2013 to 2015 and, most evidently, in the year 2017 with up to 1.4 Tg CH₄ per month. These increased emissions during certain years do not agree with the observations though, which is why the peak emissions in the posterior state are up to 50 % lower.

Another factor that has to be taken into account when analysing the posterior inter-annual variabilities is the inconsistent availability of observations throughout the different years. For years during which only limited information from the measurements is available, the posterior state is likely to be closer to the prior.

5 Conclusions

We designed an inversion system to constrain CH₄ surface fluxes in high northern latitudes based on surface observations in Arctic regions and an atmospheric transport model. Extensive sensitivity tests were carried out to comprehensively assess the methane emissions and uptake, respectively, from different CH₄ sources and sinks during the years 2008 to 2019. We aimed to reduce uncertainties in current bottom-up estimates and thereby gain a more accurate understanding of the extent, seasonality and inter-annual trends of methane emissions in Arctic and sub-Arctic regions. In order to achieve that, we computed a total of 5000 posterior states of posterior methane fluxes with varying uncertainties in the observations, background mixing ratios and prior flux estimates and evaluated their plausibility to get a reliable assessment of the methane emissions in high northern latitudes.

The atmospheric observations used for this study included both quasi-continuous and discrete measurements from 41 observation sites in different Arctic nations. We found that this observation network is not sufficient to satisfactorily constrain most CH₄ sources and sinks in high northern latitudes over the whole period of interest. Only wetland emissions are adequately constrained in North America and sporadically in Russia, although with inter-annual variabilities. It is therefore not possible to reduce uncertainties in most CH₄ emission sources and sinks occurring in high-northern-latitude regions to a substantial extent. Besides, a considerable share of the observations is used by the inversion to constrain the background mixing ratios. This share could be reduced by improved initial CH₄ mixing-ratio fields, which would allow for lower uncertainties in the background. Moreover, additional stations at the sub-Arctic boundary would be necessary to better constrain transport from CH₄ hotspots such as China, India and the Middle East. Whether the additional stations within the buffer zone, which in this case included latitudes from 30^∘ N, would improve the constraints on high-northern-latitude regions should be investigated further. Increased observation resolution instead of the monthly means used in the present work should be explored as well to potentially improve the constraints on fluxes, although at higher temporal resolution, missed emission peaks in simulations compared to observations could lead to depreciated emission estimates.

The obtained posterior CH₄ fluxes were, in comparison, predominantly lower than the prior estimates, though still higher than comparable posterior results from variational inversion set-ups. In North America, the average total methane fluxes were reduced by around 11 Tg yr⁻¹ with a corresponding uncertainty reduction of 26 %. In East Eurasia and the Arctic, fluxes were reduced by 2 and 3 Tg yr⁻¹, respectively, with uncertainty reductions of 13 % (East Eurasia) and 12 % (Arctic).

Significant changes in the seasonal cycles of the methane emissions could not be observed in either of the regions studied. Minor shifts in the seasonal cycles in certain years were exclusively influenced by CH₄ emissions from high-latitude wetlands.

Inter-annual differences were most significant in North America, where the largest discrepancies between the prior and posterior state could be observed since the region is well constrained. Whereas the highest peak emissions in the prior state took place at the end of the period of interest, raised methane emissions in the posterior state were predominantly observed at the start. The CH₄ peak emissions were hereby also up to 4.3 Tg per month lower. Those differences can be explained by a reduction in emissions from biomass burning, which appear to have been overestimated for certain years (e.g. 2014 and 2017) in the prior estimate.

The wetland emissions in North America showed a small decreasing trend between 2008 and 2019, whereas the CH₄ emissions from wetlands in East Eurasia were slightly increasing within the period of interest. Since most regions in the Arctic and sub-Arctic were poorly constrained by the inversion, most methane emission sources as well as soil oxidation did not show any significant trends in the period under study.

To get a conclusive understanding about the magnitude, long-term trend and seasonal variability in methane emissions in the entire Arctic region, it would be beneficial to expand the observation network, especially in Eurasia, to better constrain the area for future works. Complementary approaches bringing fixed and mobile platforms (ships, aircraft, trains, etc.) together should also be explored to refine our understanding of the regional Arctic budget (e.g. Pankratova et al., 2022; Berchet et al., 2020; Pisso et al., 2016; Thornton et al., 2020). New satellite platforms may also in the future expand our coverage of Arctic methane emissions, even though technical difficulties (albedo, clouds, etc.) hamper our capability of using high-latitude satellite retrievals.

This is particularly important since ongoing environmental changes due to rising temperatures in high northern latitudes are affecting natural sources and sinks of CH₄, further complicating the estimation and prediction of Arctic methane emissions, their contribution to the global budget, and the resulting potential climate feedback.