Flow-dependent observation errors for greenhouse gas inversions in an ensemble Kalman smoother

Atmospheric inverse modeling is the process of estimating emissions from atmospheric observations by minimizing a cost function, which includes a term describing the difference between simulated and observed concentrations. The minimization of this difference is typically limited by uncertainties in the atmospheric transport model rather than by uncertainties in the observations. In this study, we showcase how a temporally varying, flow-dependent atmospheric transport uncertainty can enhance the accuracy of emission estimation through idealized experiments using an ensemble Kalman smoother system. We use the estimation of European CH₄ emissions from the in situ measurement network as an example, but we also demonstrate the additional benefits for trace gases with more localized sources, such as SF₆. The uncertainty in flow-dependent transport is determined using meteorological ensemble simulations that are perturbed by physics and driven at the boundaries by an analysis ensemble from a global meteorology and a CH₄ simulation. The impact of direct representation of temporally varying transport uncertainties in atmospheric inversions is then investigated in an observation system simulation experiment framework in various setups and for different flux signals. We show that the uncertainty in the transport model varies significantly in space and time and that it is generally highest during nighttime. We apply inversions using only afternoon observations, as is common practice, but also explore the option of assimilating hourly data irrespective of the hour of day using a filter based on transport uncertainty and taking into account the temporal covariances. Our findings indicate that incorporating flow-dependent uncertainties in inversion techniques leads to more accurate estimates of GHG emissions. Differences between estimated and true emissions could be reduced more effectively by 9 % to 82 %, with generally larger improvements for the SF₆ inversion problem and for the more challenging setup with small flux signals.

1 Introduction

Atmospheric greenhouse gas (GHG) inversions use observed atmospheric GHG concentrations to estimate surface fluxes. This helps to verify emission reduction targets or the fulfillment of the Paris Agreement, monitor substances whose emissions are prohibited or regulated by the Montreal Protocol, and understand global and regional carbon cycles. The independence of such top-down inversions could also support the development of future emission inventories. However, GHG inversions are subject to various uncertainties, which are expressed by discrepancies among inversion results (e.g., Brunner et al., 2017; Monteil et al., 2020; Petrescu et al., 2023). Current inversion systems are mostly based on Bayes' theorem and solve the inversion problem by minimizing a cost function with two components: one component penalizes deviations from the a priori state, while the other penalizes differences between simulated and observed concentrations. The weighting of the two terms is determined by corresponding error covariance matrices, which define the magnitude and correlation structure of errors related to these two components. Errors in the second term, commonly denoted by R and referred to as the observation representation error or model–data mismatch error (MDM), include all processes that contribute to discrepancies between observations and model, such as aggregation and representation errors, measurement errors, and atmospheric transport errors (Kaminski et al., 2001; Engelen et al., 2002). A major source of uncertainty arises from the transport uncertainty in the atmospheric transport model (Munassar et al., 2023). Errors in the atmospheric transport lead to corresponding errors in modeled tracer concentrations, which may ultimately lead to erroneous flux estimates if not properly accounted for. Constructing the MDM is therefore a crucial yet challenging task, especially given the limited understanding of the errors associated with the representation of atmospheric transport in numerical models. Although these errors are recognized as significant, a gap still remains in adequately addressing them within inversion systems, mainly due to the fact that determining atmospheric transport uncertainties typically requires expensive meteorological ensemble simulations.

The importance of transport errors in inverse emission estimation was recognized already by Enting (1993), but systematic studies of their impacts and possible solutions were published only much later. Lin and Gerbig (2005) tested the effects of unaccounted wind errors, directly determined from radiosonde data, on inversion results and demonstrated that these uncertainties can result in significantly biased flux estimates. In a subsequent study, the authors focused on potential biases in inversion results associated with a misrepresentation of vertical mixing in the atmospheric boundary layer (Gerbig et al., 2008). Lauvaux et al. (2009) investigated transport errors in the context of inversions, determining covariances between two stations in southwest France from a meteorological ensemble simulation, and showed that the effective information content of the observations in an inversion is significantly reduced when considering these covariances. Several subsequent studies investigated transport errors and their characteristics for greenhouse gases, especially CO₂, without directly examining the impacts in an inversion framework. Using global simulations, Liu et al. (2011) demonstrated that the CO₂ transport uncertainty is highest in the tropics and in regions with the highest emissions from fossil fuels. Miller et al. (2015) examined the importance of transport errors compared to emission signals, analyzing monthly biases and associating them with meteorological conditions. Further studies explored the predictability of CO₂ (Polavarapu et al., 2016; Kim et al., 2021) or the significance of errors in transport and boundary conditions for estimating the terrestrial carbon sink in limited-area simulations (LAM; Feng et al., 2019). Two recent studies provided more detailed investigations of the transport error and its characteristics. Chen et al. (2019) focused on a sensitivity analysis of transport errors to errors in emission fluxes and initial and boundary data in LAM simulations, while McNorton et al. (2020) focused on the biogenic feedback to transport errors in global simulations. A comprehensive study on the impacts of transport uncertainties in inverse CO₂ estimation in an urban context was recently presented by Ghosh et al. (2021). They tested various transport uncertainties in synthetic CO₂ flux inversions with pseudo-observations in an urban area with a dense observation network. They compared conventional parameterizations with transport uncertainties sampled from an ensemble of simulations, both with and without spatiotemporal covariances. They demonstrated the importance of considering covariances in transport uncertainties when using data from a dense network. A method to partially account for the effects of different meteorological situations on the transport uncertainty was applied in Bergamaschi et al. (2022). In their European CH₄ inversions they parameterized MDM as a function of wind speed. Their approach assigns larger uncertainties to observations at low wind speeds, assuming that under such conditions, local emissions may have a greater influence on the observed concentrations but are not well-represented in the model.

This paper presents a comprehensive examination of flow-dependent transport uncertainties and proposes a simple method for including them in atmospheric inversions. The method is showcased for the estimation of European CH₄ inversions using the ensemble Kalman smoother system CTDAS-ICON-ART introduced by Steiner et al. (2024). Making use of meteorological ensemble simulations, we investigate the spatial, vertical and temporal characteristics of the transport error for CH₄ and investigate the impact of realistic transport uncertainties in inversion results in an observation system simulation experiment (OSSE). We use synthetically generated CH₄ observations (pseudo-observations) mimicking the observations from the current network of European in situ stations but also demonstrate the advantages for trace gases characterized by a more irregular emission distribution, such as SF₆. We evaluate the new approach in different setups and for various flux signals. Furthermore, we demonstrate how information on flow-dependent uncertainties may be used to assimilate hourly observations as opposed to afternoon observations only, as is frequently done. For this we use a filter based on transport uncertainty and take into account the temporal covariances. Finally, we also demonstrate the impacts of the new transport uncertainty on European CH₄ emission estimates using real observations. All inversions presented in this study are listed in Tables 1 and 2.

2 Model description and methodology

We conducted atmospheric transport simulations of CH₄ using the Icosahedral Nonhydrostatic (ICON) atmospheric modeling framework (Wan et al., 2013; Zängl et al., 2015; Pham et al., 2021) with the ART extension for passive and reactive tracers (Rieger et al., 2015; Weimer et al., 2017; Schröter et al., 2018). Our simulations were configured following the setup described in Steiner et al. (2024), employing an R3B6 limited-area grid (Δx≈26 km) covering Europe with 60 vertical levels. In the simulations, we use a time step of 120 s and nudge the simulation weakly towards the driving reanalysis data. In the inversion step, the 21 344 emission regions (grid cells) in the state vector are optimized using the fixed-lag ensemble Kalman smoother implemented in CTDAS with an assimilation window length of 9 d and a lag of 2 assimilation windows. The number of members in the emissions ensemble is 192. Unlike our previous study with a single forward simulation, we created a meteorological ensemble (see Sect. 2.1) with 10 ensemble members driven by perturbed meteorological boundary conditions and model physics, as well as by perturbed CH₄ boundary conditions from the same global ensemble simulation that provided the meteorological boundary conditions. Each ensemble member contains two CH₄ tracers: a background tracer (${\mathrm{CH}}_{\mathrm{4}}^{\text{bg}}$) representing the perturbed CH₄ boundary conditions and an emission tracer (${\mathrm{CH}}_{\mathrm{4}}^{\text{emis}}$) representing the additional CH₄ emitted within our European model domain. These ${\mathrm{CH}}_{\mathrm{4}}^{\text{emis}}$ tracers experienced the same (unperturbed) CH₄ emissions, such that the ensemble of 10 ${\mathrm{CH}}_{\mathrm{4}}^{\text{tot}}={\mathrm{CH}}_{\mathrm{4}}^{\text{bg}}+{\mathrm{CH}}_{\mathrm{4}}^{\text{emis}}$ tracers represents only the effect of transport uncertainty.

2.1 ICON-ART ensemble simulations

To generate a meteorological ensemble, we ran the model with 10 members, each driven by the output of an experimental Ensemble of Data Assimilation (EDA) simulation conducted at ECMWF (experiment ID hyfd), which included the global surface emissions inversion and transport of CO₂ and CH₄ (McNorton et al., 2020). The experiment was conducted for the month of July 2019 and was carried out using the Integrated Forecasting System (IFS) cycle 48r1 with 10 ensemble members at a 25 km resolution (Tco399), in which model physics (in the form of Stochastically Perturbed Parameterization Tendencies (SPPT; Leutbecher et al., 2017)), observations (observations of both weather and GHG) and sea surface temperature were perturbed. Additionally, the GHG emissions were also perturbed, with the perturbations sampled from a log-normal distribution with an a priori uncertainty of 100 %. The perturbations were not spatially correlated, but a correlation length of 100 km was assumed in the background covariance matrix for CH₄ emissions. The background error covariance matrix for GHG concentrations was static, that is, it was not updated during the experiment and was based on differences between forecasts with different lead times obtained from a climatological sample (see NMC method Parrish and Derber, 1992). The background error covariance matrix for the numerical weather prediction (NWP) fields was based on the archived IFS Ensemble of Data Assimilations (EDA) o-suite (experiment 0001) for the same period. Greenhouse Gases Observing Satellite (GOSAT), Infrared Atmospheric Sounding Interferometer (IASI) and TROPOspheric Monitoring Instrument (TROPOMI) retrievals of CH₄ were assimilated. The emissions were optimized independently in each 12 h window.

In the ICON ensemble simulations, in addition to the perturbed driving data (meteorological variables and CH₄ concentrations to drive the ${\mathrm{CH}}_{\mathrm{4}}^{\text{BG}}$ tracer), we also applied perturbations to model physics tuning parameters as implemented in ICON for the ensemble data assimilation scheme of the German weather service. Together, these perturbations are expected to represent the typical level of uncertainties present in state-of-the-art meteorological analysis products.

The perturbed CH₄ concentration fields of the driving data were used to initialize and drive the background tracers, ${\mathrm{CH}}_{\mathrm{4}}^{\text{bg}}$, in our simulations. In this study, we did not optimize the background concentrations, but the perturbations of the different background fields were part of the artificial transport error. In addition to the ${\mathrm{CH}}_{\mathrm{4}}^{\text{bg}}$ tracer, we also incorporated a ${\mathrm{CH}}_{\mathrm{4}}^{\text{emis}}$ tracer into our ensemble simulations using ICON-ART, transporting the emitted signal with (unperturbed) emissions introduced via the online emission module (OEM; Jähn et al., 2020; Steiner et al., 2024). The computational cost, measured in node-hours, for the 10-member ensemble simulation was about 1.7 times the cost of regular inversions. This means that the total computational cost (a priori meteorological ensemble and inversion) was 2.7 times the cost of a regular inversion, which is a considerable but not prohibitive increase.

2.2 CTDAS inversion setup

The setup of the CTDAS-ICON-ART inversion system aligns with Steiner et al. (2024). We optimize the emission scaling factors in each grid cell of the R3B6 grid while applying exponential decaying correlations with a length scale of 200 km in the a priori error covariance matrix. However, optimization is restricted to a single emission category, namely the total CH₄ emissions, which is the sum of anthropogenic and various natural emissions.

Since we compare different approaches in representing the MDM uncertainty, the question of how these uncertainties should be scaled to enable a fair comparison arises. Since the innovation chi-squared statistics are a common diagnostic to judge the validity of the uncertainty assumptions made in an inversion (Berchet et al., 2015; Michalak et al., 2017), we contend that comparability is best achieved when the uncertainties in each inversion are scaled such that the innovation chi-squared value is 1. The chi-squared metric delineates the ratio of a priori residuals between observed and simulated concentrations to the total a priori variance in the observation space, accounting for the projected a priori variance from the observation space alongside the model–data mismatch:

$$\begin{array}{}\text{(1)}& {\displaystyle }{\mathit{\chi }}_{\text{innov}}^{\mathrm{2}}={\displaystyle \frac{\mathrm{1}}{n}}({\mathit{y}}^{\mathit{o}}-\mathbf{H}\mathbf{(}{\mathit{x}}^{\mathit{b}}\mathbf{\left)}{)}^T\right({\mathbf{HP}}^{\mathbf{b}}\mathbf{H}+\mathbf{R}{)}^{-\mathrm{1}}({\mathit{y}}^{\mathit{o}}-\mathbf{H}\mathbf{(}{\mathit{x}}^{\mathit{b}}\mathbf{\left)}\right).\end{array}$$

Here, y^o denotes observed data, H the observation operator, x^b the a priori (background) state, P^b the a priori error covariance matrix, R the observation error covariance matrix (or MDM) and n the number of observations. The requirement of a chi-squared value close to 1 results in some differences in the magnitude of the MDM (indicated as α in Tables 1 and 2), but we believe this is a better approach than requiring, for example, that the mean MDM be identical between different inversions.

Table 1 Overview of all inversions used in this study. The column MDM indicates whether the flow-dependent MDM (FD) or the constant (C) MDM is used. α is the scaling factor from Eq. (2). The column Emis. indicates whether the CH₄ emissions or SF₆ emissions are used. ER stands for error reduction.

Table 2 Continuation of Table 1 for inversions with the synthetic setup and hourly observations, as well as for inversions with real observations.

2.3 Idealized experiments

In our idealized setup, one of these 10 ensemble members was considered to represent the true meteorology and used to generate the pseudo-observations of CH₄. All inversions were then performed using a different ensemble member in order to mimic the fact that the simulated meteorology (and the corresponding transport of CH₄) in general deviates from the true meteorology within the range of uncertainty in state-of-the-art meteorological analyses.

The ensemble spread (corresponding to the standard deviation) of the ${\mathrm{CH}}_{\mathrm{4}}^{\text{tot}}$ tracers was sampled at each station location and used to determine the temporally varying flow-dependent model–data mismatch (MDM), replacing the static mdm used in our previous inversion study (Steiner et al., 2024). In contrast to our previous study, we optimized only the emissions in the idealized setup. We did not optimize the background concentrations because the differences in background CH₄ concentrations introduced by deviating from true meteorology and using perturbed background ${\mathrm{CH}}_{\mathrm{4}}^{\text{bg}}$ concentrations in the driving data are part of the artificially created transport error and contribute to the ensemble spread that determines the MDM. However, if there were systematic biases in background CH₄ in an application with real data, it would still be necessary to optimize background concentrations together with the emissions. The performance of the inversions with the flow-dependent MDM was compared with inversions using the static MDM as used in Steiner et al. (2024). To ensure a fair comparison, the MDM was scaled in each inversion so that the innovation chi-squared value in each inversion was 1.

In the reference setup (inversions fc01 and fc02), we utilized an a priori variance of 0.07 (unitless) at the grid cell level, which corresponds to an emission uncertainty of 26 %. With such a low emission signal it is easier to demonstrate the benefits of our new approach, as the relative importance of transport uncertainties is larger. Since the same variance of 0.07 was used for generating the pseudo-observations (see Sect. 2.5), the reference setup corresponds to the optimal situation where the assumed a priori uncertainties are in perfect agreement with the true emission errors.

Starting from the reference setup, two sets of sensitivity experiments were conducted. In the first set, different combinations of a priori variances ranging from 0.007 to 0.7 were tested in order to analyze the impact of a priori assumptions deviating from the true emission error (inversions uc01 to uc16 and us01 to us16). In the second set, the true variances were varied between 0.07 and 1.0 in order to assess the potential of improving the inversion results using our proposed method for different ratios between flux uncertainties and transport uncertainties. All inversions were performed for two different emission fields. One emission field (used in inversions fc01 to fc16 and uc01 to uc16) is consistent Steiner et al. (2024) and comprises the sum of anthropogenic emissions from EDGARv6 (Crippa et al., 2021) together with various natural sources. These natural emissions include peatlands and mineral soils from JSBACH-HIMMELI version 2 (Jena Scheme for Biosphere–Atmosphere Coupling in Hamburg coupled to HelsinkI Model for Methane build-up and emission for peatlands; Raivonen et al., 2017; Reick et al., 2013); inland water emissions from the Global Carbon Project (GCP) CH₄ dataset (Saunois et al., 2020); and additional sources including termites (Saunois et al., 2020), ocean emissions (Weber et al., 2019), biofuel and biomass burning from the Global Fire Emission Database 4.1 s (GFED; van der Werf et al., 2017), and geological emissions (Etiope et al., 2019) scaled to a global total of 15 Tg. The second emission field (in the inversions fs01 to fs16 and us01 to us16) contains the same total amount but is spatially distributed according to the SF₆ emission field for the categories NFE and PRU in EDGARv7 (Crippa et al., 2022). The latter shows a more irregular emission distribution than CH₄ and, hence, serves as an additional test case with expectedly larger gradients. Figure 1 shows the spatial distribution of the two emission fields in our model domain.

Figure 1 Spatial distribution of CH₄ emissions (a), as well as CH₄ redistributed to follow EDGAR SF₆ emissions (b), remapped onto the simulation grid. The emissions are representative of the period of 2–11 July 2019.

In another experiment (inversions fc01 to fc12 and fs01 to fs12), we compared the aforementioned inversions where we assimilated only afternoon averages (nighttime averages for mountain stations) with inversions where we assimilated hourly observations filtered based on the ensemble spread of CH₄ concentrations. The filter excluded observations during times when the ensemble spread exceeded 5 ppb. This resulted in the exclusion of ca. 22 % of the 19 530 observations. Figure 2 shows the fraction of excluded observations as a function of the hour of the day. For these inversions, however, the assumption of temporally uncorrelated errors was no longer valid and was accounted for by introducing off-diagonal elements in the MDM error covariance matrix. The temporal error correlations were computed from the ensemble spread. For each of the hourly observations, the error correlation with the observations at the same station in the next 36 h was considered.

Figure 2 The fraction of pseudo-observations that were excluded from all available observations as a function of the hour of the day (UTC) by applying a threshold of 5 ppb for the ensemble spread in the ensemble of CH₄ concentrations.

2.4 Experiment with real-data

The setup for the application with real data (inversions real_c and real_f) closely followed the setup for the idealized experiment with the assimilation of daily afternoon or nighttime means (e.g., as in fc01). We used the same setup for the ICON simulations with a grid over Europe and the same state vector in the inversions. The transport uncertainties for the MDM were also derived from the same ensemble simulation. However, in order to stay closer to the real meteorology and background CH₄ concentrations, the forward simulations of the inversion were driven at the domain boundaries by ERA5 reanalysis data for meteorology (Hersbach et al., 2020) and by the CAMSv22r2 product (available via https://ads.atmosphere.copernicus.eu/, last access: 18 April 2024) for background CH₄ concentrations. Using the CAMSv22r2 product was necessary, as the CH₄ mole fractions from the experimental ensemble simulation had biases that were too large.

2.5 Pseudo-observations

Pseudo-observations, following the methodology outlined in Steiner et al. (2024), are generated with a forward simulation of ICON-ART, wherein the ${\mathrm{CH}}_{\mathrm{4}}^{\text{tot}}$ tracer concentrations are sampled at the station locations. In this simulation, the emission field is perturbed using a set of true scaling factors, which we aim to reproduce as accurately as possible through the inversions starting from unscaled emissions. The true scaling factors are a field of spatially correlated random perturbations with a correlation length of 200 km. In our standard configuration, a variance of 0.07 is applied to generate this perturbation field. We systematically explored different configurations, varying the true and a priori variances, with true variances ranging from 0.007 to 1.0. In addition, to mimic measurement noise, a 2 ppb noise was introduced. Pseudo-observations were generated at the stations available in the dataset of the European Obspack 2022-1 L2 release (ICOS RI et al., 2022) for the year 2019.

2.6 Real observations

In the application with real observations, the Obspack dataset of quasi-continuous in situ observations from 28 stations was used. Most of the stations are members of the atmosphere network of the Integrated Carbon Observation System (ICOS; Heiskanen et al., 2022). As in our previous study, we distinguish mountain sites from sites in flat terrain. Stations where the model topography was more than 200 m lower than the actual topography (due to coarse grid representation) were classified as mountain stations. For sites in flat terrain, only daytime (11:00 to 16:00 LT, local time) mean values were assimilated, as is usually done in atmospheric inverse modeling to avoid difficulties in representing shallow nocturnal boundary layers. In contrast, only nighttime mean values between 23:00 and 06:00 LT were assimilated for mountain sites, as these are least influenced by pollution from daytime upslope valley winds, which are difficult to represent in a coarse-resolution model. The height at which the model output was sampled was different for mountain sites than for sites in flat terrain to account for the fact that the smooth model topography typically underestimates the real altitude of mountain sites. For sites in flat terrain, the (relative) height of the observation above ground was preserved, whereas for mountain sites a height in between the height relative to the model topography and the absolute height of the station was chosen. Mountain stations are indicated in Fig. 8 with triangles and lowland stations with circles. For stations located on a hill but still within the daytime boundary layer (e.g., Beromünster), only measurements in the afternoon were used as for stations in flat terrain, but the vertical sampling of model fields was done in the same way as for mountain stations in order to maintain a realistic relative distance from near-surface emissions.

2.7 Observation error

For each pseudo-observation, we calculate the ensemble spread of ${\mathrm{CH}}_{\mathrm{4}}^{\text{tot}}$ in the ensemble. This ensemble spread is then incorporated into the flow-dependent model–data mismatch (MDM) and scaled to achieve an optimal innovation chi-squared value (see Sect. 2.2). Specifically, the flow-dependent MDM is computed as

$$\begin{array}{}\text{(2)}& \text{MDM}=\mathit{\alpha }\cdot \text{SD}\left({\mathrm{CH}}_{\mathrm{4}}^{\text{tot}}\right)+\mathrm{2}\left(\mathrm{ppb}\right),\end{array}$$

where the factor α varies depending on the inversion (to keep an innovation chi-squared value of 1) but remains constant across all stations and time steps within one inversion. The additional term of 2 ppb accounts for the 2 ppb noise introduced into the pseudo-observations (see Sect. 2.5).

Inversions with the flow-dependent MDM are compared to inversions with the static MDM, which was implemented following a similar principle as that outlined in Steiner et al. (2024). The static MDM varies between stations but remains constant over time. In this study, we calculated the average of the flow-dependent MDMs for the observations assimilated at each station, adjusting them with a factor to maintain an innovation chi-squared value of 1.

The MDMs for application to real observations were created in the same manner, with the only difference being that the factor x was chosen to be station-dependent, which allows us to achieve a chi-squared value of 1 for each station separately. This became necessary because, unlike in the synthetic setup, some regions exhibited significantly larger biases in the background concentrations than other regions. With this adjustment, stations where a large bias occurred had a lower weight than stations with good a priori agreement.

To construct the R matrix for inversions assimilating hourly observations, we calculated, for each station and hour of the day, the mean temporal correlation with observations for the next 36 h during the inversion period. We then fitted a function to these data for each station and hour of the day. The function is a combination of exponential decay and a Gaussian distribution:

$$\begin{array}{}\text{(3)}& \exp {\displaystyle \frac{-\mathrm{\Delta }t}{a}}+b\exp {\displaystyle \frac{-(\mathrm{\Delta }t+c{)}^{\mathrm{2}}}{\mathrm{2}{d}^{\mathrm{2}}}},\end{array}$$

where the parameters a, b, c and d are fitted. We chose this fit because it is able to represent the decay in the first hours and the correlated errors between two nights at lowland stations (see Fig. 7). The mean value of the decay time a over all stations and hours is 7.8 h, with large variability indicated by a standard deviation of 4.0 h. The inversion of the R matrix, which only considers the covariances fitted to the next 36 h of observations, is very unstable due to poor conditioning. As a result, unrealistic results are produced when inverting HP^bH+R. To address this issue, we conditioned the R matrix by multiplying the fitted covariances by a factor that exponentially decreases with time ($\exp (-\mathrm{\Delta }t/\mathrm{24}$ h)), as proposed also by Ghosh et al. (2021).

3 Results and discussion

3.1 Characteristics of the transport error

We start the analysis by illustrating examples of the ratio between flux and transport uncertainties, which can be interpreted as the signal-to-error ratio of the flux signal: fluxes can only be retrieved reliably if this ratio is larger than one. The ratio was computed as the ratio between the spread in ${\mathrm{CH}}_{\mathrm{4}}^{\text{emis}}$ concentrations of the flux ensemble (used in the inversions) to the spread in ${\mathrm{CH}}_{\mathrm{4}}^{\text{tot}}$ concentrations of the meteorological ensemble. The time series of this ratio for Cabauw (lowland site) and Monte Cimone (mountain site) at two different magnitudes of flux signals, one with a variance of 0.07 and one with a variance of 1.0, are depicted in Fig. 3. This illustration shows that even in the scenario with a small flux error variance of 0.07, the signal of flux uncertainties is often still stronger than the signal of the transport error (ratio above 1). The low values at the beginning of the time series are due to a spinup effect: while the emission signal is still extremely small (${\mathrm{CH}}_{\mathrm{4}}^{\text{emis}}$ has not yet reached the stations), the spread in ${\mathrm{CH}}_{\mathrm{4}}^{\text{bg}}$ is already fully developed due to the perturbed initial and boundary conditions (ICBC).

Figure 3 Time series of the ratio between the emission signal and transport uncertainty at the stations Cabauw (a) and Monte Cimone (b) for two different levels of emission signal.

A snapshot of the ensemble spread of CH₄ at the lowest model levels at an arbitrary time step (Fig. 4) highlights the spatial variability in the transport error. The spread in the tracer of emitted ${\mathrm{CH}}_{\mathrm{4}}^{\text{emis}}$ reveals hotspots of large uncertainties, particularly downwind of strong CH₄ sources, but also elongated features of high uncertainty likely associated with frontal zones. This underscores the influence that atmospheric flow conditions have on the structure of the transport error. The spread in the background tracer ${\mathrm{CH}}_{\mathrm{4}}^{\text{bg}}$ is of similar magnitude but is much smoother.

Figure 4 Map of the ensemble spread of ${\mathrm{CH}}_{\mathrm{4}}^{\text{emis}}$ (a) and ${\mathrm{CH}}_{\mathrm{4}}^{\text{bg}}$ (b) in the lowest model level on 27 July 2019, 15:00 UTC. The arrows show the ensemble mean wind in the lowest model level.

Further insight into the structure of the transport uncertainty is obtained by plotting time series of vertical profiles of the ensemble spread. Figure 5 shows such a time series for the Dutch station Cabauw. It shows distinct periodic increases in the ensemble spread near the surface during the nights. These increases reflect the uncertainties associated with the simulation of boundary layer processes and their impact on boundary layer heights, which has a particularly strong impact on tracer concentrations in shallow nocturnal boundary layers. The periodic increases in uncertainty are occasionally superimposed by larger-scale events, such as the one from 9 to 12 July. The higher uncertainties at night support the common practice in GHG inversions to use only daytime observations, such as afternoon averages. At the same time, it is conceivable that the use of a meteorological ensemble also provides an opportunity to filter observations based on transport uncertainty rather than time of day and thereby use the information provided by the observations more efficiently.

Figure 5 Time series of vertical profiles of the standard deviation of ${\mathrm{CH}}_{\mathrm{4}}^{\text{bg}}$ (a, b) and ${\mathrm{CH}}_{\mathrm{4}}^{\text{emis}}$ (c, d) at Cabauw. The dashed lines in panels (a) and (c) indicate the upper boundaries of the plots in panels (b) and (d). The black line in the right-column plots indicates the inlet height at Cabauw.

The mean diurnal ensemble spread of CH₄ concentrations is depicted in Fig. 6, distinguishing between lowland and mountain stations. Each figure presents two scenarios: one where the ensemble is generated with perturbed initial and boundary conditions only (IC/BC, gray) and one where in addition to IC/BC the model physics is also perturbed (IC/BC and STTP, red). The spread attributed solely to perturbed IC/BC accounts for approximately 50 % of the spread in the two tracers when both perturbed model physics and IC/BC are considered. The higher spread in the nocturnal boundary layer at lowland sites and the peak during the early morning hours are caused approximately equally by the perturbed model physics and by the perturbed IC/BC conditions. At mountain stations, the spread is nearly constant over the day. This indicates that the ensemble is not fully capable of estimating the uncertainty at mountain sites, since a fundamental problem, the misrepresentation of thermally induced flow in the afternoon, is inherent to all ensemble members. This supports the common practice of assimilating only nighttime observations at mountain stations.

Figure 6 Diurnal profiles of the standard deviation of ${\mathrm{CH}}_{\mathrm{4}}^{\text{bg}}$ (a, b) and ${\mathrm{CH}}_{\mathrm{4}}^{\text{emis}}$ (b, c) concentrations in the ensemble for an ensemble simulation with perturbed IC/BC only (gray) and an ensemble simulation with perturbed IC/BC and perturbed model physics (red).

The temporal correlations illustrated in Fig. 7 (in analogy to Fig. 4 in Lauvaux et al., 2009) provide information on the temporal structure of the transport error. Two different periods of the day are examined separately (00:00 to 08:00 UTC and 12:00 to 20:00 UTC) to emphasize the differences between daytime and nighttime conditions. Each line shows the error correlation (y axis) of an hourly mean observation with the observations of the next 36 h (x axis), with separate analyses conducted for the lowland and mountain stations. At lowland stations, the nocturnal values exhibit significant error correlations with subsequent night-time observations but lack correlations with daytime values. Similarly to Lauvaux et al. (2009), the correlations with observations from the following night show that on the one hand, the nocturnal error structures are determined by static parameters that cause similar errors in different nights, but, on the other hand, the system is also sensitive to disturbances, as the correlations remain below 0.5. It can also be seen that the decline in correlations for the 08:00 UTC observation occurs earlier than for the midnight observation, which is a result of the earlier breakdown of the nocturnal boundary layer for the 08:00 UTC observation. For mountain sites, this pattern of recurring error correlations is much less pronounced. The daytime observations at both stations exhibit exponentially decaying correlations without a subsequent increase in the afternoon of the following day, contrasting with the observed nocturnal correlations.

Figure 7 Mean temporal correlations at Cabauw (a, b) and Monte Cimone (c, d) for observations at 00:00 to 08:00 UTC (a, c) and 12:00 to 20:00 UTC (b, d) with observations for the next 36 h. The error correlations over time (x axis) are indicated by one line for each observation.

We finally examined the correlations between the CH₄ transport uncertainty in our ensemble simulation and the ensemble spread of wind speed and direction, which could also be obtained from a reanalysis product without the need for an ensemble forward simulation. However, our CH₄ transport errors show only weak correlations with the ensemble spread of wind speed (0.15 on average, ranging from −0.19 to 0.47 at individual stations) and direction (0.10 on average, ranging from 0.14 to 0.42 at individual stations). This result supports the conclusions of Miller et al. (2015), who discussed and tested the incorporation of transport uncertainty into the MDM without the need for expensive ensemble simulations. They computed correlations between monthly CO₂ biases in atmospheric transport (relative to CO₂ boundary layer enhancements) and individual meteorological variables in a global ensemble simulation. The strongest correlations were found with inverse temperature over terrestrial regions (0.45) and with zonal winds over the oceans (0.29). However, many errors could not be explained by a single explanatory variable.

3.2 Flow-dependent observation error in an idealized setup

Figure 8 shows the true scaling factors alongside the optimized factors obtained from inversions using both the static (fc02) and flow-dependent MDM (fc01). Both inversions seem to be similarly successful at reproducing the large-scale patterns of the true state, especially in central Europe where emission fluxes and observation density are high. To better compare the quality of the results, Fig. 9 illustrates the improvement achieved by the flow-dependent MDM compared to the static MDM, both in terms of scaling factors (panel a) and emission fluxes (panel b). The predominance of green colors suggests a significant overall improvement when using the implementation of the new MDM. Since we perform these inversions in a synthetic setup where the ground truth is known, we can compute the flux error reduction precisely. Summed over the entire domain, this reduction amounts to 20.12 kg s⁻¹ (or 8.7 % of the a priori total error) with the static MDM and 36.57 kg s⁻¹ (or 15.8 % of the a priori total error) with the flow-dependent MDM, which corresponds to an improvement by 82 %.

Figure 8 True (a) and a posteriori scaling factors for the inversions with the static MDM (b) and the flow-dependent MDM (c).

Figure 9 Improvement in the error reduction for inversions using the flow-dependent MDM vs. the static MDM in terms of scaling factors (a) and emission flux (b). Green color indicates that the inversion with the flow-dependent MDM performs better, while red color indicates that the inversion with the static MDM performs better.

Figure 10 illustrates the relationship between error reduction and the ratio of the true to the a priori variance for the inversions with flow-dependent and static MDMs. In this plot, the true variance remains constant at 0.07 while the a priori variance varies, representing different levels of freedom to adjust the state in the inversion process. The analysis includes both CH₄ (uc01 to uc16) and SF₆ (us01 to us16) emission patterns. In the left half of the figure, the a priori variance is larger than the true variance. As a result, the system has too much freedom, resulting in strong adjustments of the state. This tendency to overfit leads to poor performance, particularly evident with the static MDM, where the error may become even larger than the a priori error due to over-fitting to biased observations. The right half of the figure shows situations where the a priori uncertainty is too low and the cost of the state adjustment is correspondingly high. In contrast to the previous situation with too high uncertainty, the degradation in performance is comparatively slow, and, in some instances, there may even be a slight performance gain, notably with the static MDM. This phenomenon can be attributed to the fact that decreasing the a priori uncertainty minimizes updates in all regions, due to higher associated costs in the cost function. As a result, regions initially subject to incorrect updates remain closer to the a priori state. This partially counteracts the performance degradation resulting from reduced updates in regions that initially perform well. This effect becomes apparent as we used (correlated) random perturbation factors for the ground truth, which are normally distributed around 1. In this case, the solution frequently benefits from maintaining proximity to the a priori state. However, in scenarios featuring substantial biases within the a priori on a larger scale, this proximity would likely result in a more rapid decline in performance, as depicted on the right side of the figure.

Figure 10 Error reduction (in kg s⁻¹) for inversions using the flow-dependent MDM (circles) and the static MDM (crosses) as a function of the ratio of true variance to a priori variance. Results are shown for inversions with CH₄ emissions (a) and SF₆ emissions (b).

In the reference setup with a true variance of 0.07, the transport error is relatively large compared to the emission signal. While this makes optimization challenging for the inversion system overall, it presents greater potential for improvement with the flow-dependent MDM. To evaluate the improvement with the flow-dependent MDM across various magnitudes of the emission signal, we illustrate the relationship between the relative reduction in errors (in percentage) and the true variance (${\mathit{\sigma }}_{\text{true}}^{\mathrm{2}}$) in Fig. 11. In all these inversions (fc01 to fc12 and fs01 to fs12), a perfect assumption is made for the a priori variance, resulting in ${\mathit{\sigma }}_{\text{prior}}^{\mathrm{2}}$ being equal to ${\mathit{\sigma }}_{\text{true}}^{\mathrm{2}}$. Moving the right side of the figure towards larger ${\mathit{\sigma }}_{\text{true}}^{\mathrm{2}}$ values, the ratio of the emission signal to transport error increases, making the system a simpler problem to optimize. This is reflected in the larger error reductions. However, as the flux signal increases, the distance between the lines representing the static and flow-dependent MDM reduces; i.e., the benefit of applying a flow-dependent MDM becomes smaller. The error reduction in the inversions is generally better for the SF₆ emission pattern compared to the CH₄ emissions, and the flow-dependent MDM has a larger effect. This is to be expected, as the transport error translates into a larger tracer concentration error when the emission pattern is more heterogeneous, as in the case of SF₆ emissions.

Figure 11 Relative error reduction for inversions using the flow-dependent MDM (circles) and the static MDM (crosses) as a function of the true variance. All inversions use the same a priori covariance as the true variance. Results are shown for inversions with CH₄ emissions (a) and SF₆ emissions (b).

3.3 Hourly observations with correlated errors in an idealized setup

To evaluate the effectiveness of assimilating hourly values vs. only daily afternoon mean values (nighttime for mountain sites), we performed five different inversions, each using a different R matrix. All five inversions are performed with the CH₄ emission field with the true as well as the a priori variance being 0.07. The results are summarized in Fig. 12, which shows the total error reduction across the domain for these inversions (positive y axis), as well as the reduction in the a posteriori uncertainty compared to the a priori uncertainty (negative y axis). The first two boxes show the error reduction in the two inversions assimilating daily afternoon and night averages, with both constant and flow-dependent MDM (R_a and R_b, respectively). These two inversions correspond to the two points on the far left in Fig. 11. The next box in Fig. 12 represents an inversion with hourly instead of daily observations with a diagonal R matrix with constant values for each station, corresponding to the static MDM (R_c). The next box represents an inversion also with a diagonal R matrix but with time-varying elements sampled from the meteorological ensemble, which corresponds to the flow-dependent MDM (R_d). The last box represents an inversion where temporal covariances are included in the R matrix as off-diagonal elements based on the sampled correlations between the hourly observations (R_e) as described in Sect. 2.7. Compared to the error reductions of 9.0 % and 15.6 % for the simulations R_a and R_b with daily observations, inversions assimilating hourly observations exhibit an improved error reduction of 15.8 % (R_c), 17.7 % (R_d) and 19.9 % (R_e). Thus, besides the overall improvement, the performance of inversions at assimilating hourly observations also increases when a flow-dependent MDM is used, and it further significantly improves if temporal covariances are considered. These results are consistent with those of Ghosh et al. (2021). In their synthetic study using a dense observational network in an urban area, they observed significant improvements in domain-total emission estimates for inversions using a diagonal R matrix constructed using the ensemble spread (equivalent to R_c). They found that non-diagonal R matrices that account for covariances (such as R_e) resulted in better estimation of the spatial emission structure but that this effect diminished when fewer stations were assimilated, which more closely resembles our widely spaced station setup. The negative y axis in Fig. 12 shows the reduction in uncertainty in the a posteriori P matrix compared to the a priori uncertainty. Assimilating hourly data results in a significantly larger reduction in the uncertainty in the P matrix as a greater number of observations is used. When comparing the two inversions that assimilate daily observations, it is evident that the one utilizing the flow-dependent MDM (R_b) shows a slightly smaller reduction in uncertainty despite its much better error reduction compared to the constant MDM (R_a). This highlights that uncertainty reduction does not necessarily correlate with inversion performance. Instead, it is a result of assumptions made about error characterizations and the amount of information provided by observations used in the inversion. The same observation also applies when comparing R_c, R_d and R_e. Here it can be argued that the state of the inversions that assimilate hourly data but lack temporal covariances (R_c and R_d) is adjusted too much to the observations, thus suffering from overfitting.

Figure 12 Boxplot of the domain-total error reduction (positive y axis) and uncertainty reduction in the a posteriori P matrix compared to the a priori P matrix (negative y axis) of the five different inversions (R_a–R_e).

3.4 Effect on real emission estimates

Figure 13 compares the results for the inversions (real_f and real_c) with real data for the month of July 2019. The maps of the increments generally show a very similar pattern, with differences between the two a posteriori emissions shown in Fig. 13c. The significant downward correction over Italy obtained with the static MDM is attenuated with the flow-dependent MDM, as is the case for the Moldova/Romania region. Similarly, the upward correction over southern England is damped. In contrast, the strong upward correction over the Benelux countries is further enhanced. While the spatial patterns of adjustments differ significantly, domain-total emissions are very similar: the flow-dependent MDM results in total a posteriori emissions of 988 kg s⁻¹, while the static MDM results in emissions of 981 kg s⁻¹ (with an a priori of 1150 kg s⁻¹ in both cases). In Fig. 13d we present these differences in relation to the a posteriori uncertainty in the inversion using the static MDM. It is evident that there are only a few regions where the differences exceed the a posteriori uncertainty, but in many regions this ratio is close to 1. Since the true emissions are unknown in this case, it is impossible to tell which of the two results is closer to reality. However, based on the results from the synthetic experiments, the results obtained with the flow-dependent MDM are to be preferred. The approach using a static MDM is more likely to assign too much weight to an observation collected during a meteorologically uncertain situation and, conversely, too little weight during a situation when the meteorology is predicted well.

Figure 13 A posteriori increments for the inversions with the static MDM (a) and the flow-dependent MDM (b). The difference between the a posteriori emissions for inversions with the flow-dependent MDM and with the static MDM is shown in (c), while the same difference divided by the a posteriori uncertainty is shown in (d).

4 Conclusions

This paper presents a comprehensive examination of flow-dependent transport uncertainties in GHG inversions. Leveraging meteorological ensemble simulations, we investigate the influence of realistic, temporally varying transport uncertainties in inversion results across various setups and flux strengths and compare it to the more traditional static approximation of the transport uncertainty.

4.1 Error structure characteristics

The spatial structure of the transport uncertainty exhibits highly variable patterns, especially when considering the tracer of emitted CH₄. This underscores that static MDM is a poor approximation of the real transport uncertainty. In contrast, the uncertainty in the background signal shows a larger-scale, more homogeneous structure. Aside from isolated weather situations, transport uncertainties are typically greatest during the night, reflecting the challenges models face in simulating low nocturnal boundary layer heights. This reaffirms the prevalent practice in current inversion systems of disregarding nocturnal observations or even assimilating only afternoon values at lowland stations. However, it is also apparent that such an ensemble simulation provides the opportunity to filter observations based on the corresponding uncertainty in the model rather than the time of day and even enables the assimilation of more observations. Similar to Miller et al. (2015), we could not find a clear correlation between transport uncertainties and wind speed that would support the approach of Bergamaschi et al. (2022) of assigning larger uncertainties to observations under low-wind conditions. However, with our synthetic setup, it was not possible to test their plausible hypothesis that local sources not represented by the model due to insufficient resolution have the largest influence when wind speeds are low.

4.2 Flow-dependent MDM in inversions

We compare the inversion results with a flow-dependent MDM and the standard static MDM in an idealized setup with synthetically produced observations. Both inversions show the best performance in central Europe, where emission fluxes and observation density are large. However, inversions with the flow-dependent MDM achieve a larger overall improvement. Depending on the flux signal and emission pattern, relative improvements from 9 % (for the largest flux signal of CH₄ emissions) to 81 % and 82 % (for the smallest flux signal of SF₆ and CH₄ emissions, respectively) are achieved. However, it is crucial to note that while our study achieved large improvements with the new error description, these advancements were observed within an idealized, synthetic setup where the (artificial) transport error is internally consistent with the uncertainties derived from the ensemble spread. Furthermore, we present an analysis where we depart from the assumption of perfect a priori uncertainty in inversions and highlight the importance of making realistic assumptions about a priori uncertainty as well. In particular, overestimating the a priori uncertainty quickly leads to overfitting to the biased observations. Conversely, being too conservative in the assumptions of a priori uncertainties results in a less pronounced decrease in performance, at least as long as the a priori assumptions do not have large-scale biases.

4.3 Assimilation of hourly observations

In this study, we also assessed the effectiveness of inversions assimilating hourly vs. daily observations. For hourly observations it was necessary to account for temporal correlations in the transport error and hence to include off-diagonal elements in the R matrix. Our analysis showed that these correlations typically exhibit an exponential decay with time, with nighttime observations showing more-persistent correlations within the same night. However, we observed a notable exception at lowland stations, where correlations increased again during the following night, peaking at values of 0.5. To incorporate these correlations into our R matrix, we fit a function to the sampled correlations. This function had to be damped with an exponentially decaying factor to facilitate robust results for the inverse of the otherwise ill-conditioned matrix HP^bH+R. The results demonstrate that hourly data assimilation leads to superior performance compared to daily assimilation of observations. Within the inversions with hourly observations, the performance improved when a flow-dependent MDM was used instead of a constant MDM, and it improved even more when temporal covariances were also considered. The results also indicate that uncertainty reduction in the P matrix does not necessarily correlate with inversion performance. It is rather a result of assumptions about error characterizations and the amount of information provided by observations used in the inversion.

4.4 Inversion with real observations and its limitations

Examining the inversion results together with real observations for July 2019, we find that in certain regions, such as Italy and Moldova/Romania or southern England, the flow-dependent MDM attenuates either the downward or the upward correction, while in other regions, such as the Benelux countries or Switzerland, the upward correction is amplified. However, it is noteworthy that these differences, although often comparable in magnitude to the uncertainties, rarely reach significance relative to the a posteriori uncertainty. We have applied the flow-dependent MDM to real observations assuming that the results are improved compared to the static MDM in a similar way as in the idealized setup. However, while the ensemble spread in the idealized setup accurately captures all transport uncertainties, this is not guaranteed in a setup with real data. Further analysis would be needed to verify whether the ensemble spread provided by the ECMWF EDA NWP system adequately reflects the differences between the measurements and the simulations. Furthermore, in the idealized setup only random uncertainties are accounted for, whereas in reality there may also be systematic transport errors. Systematic errors could result, for example, from a misrepresentation of vertical mixing in stable nocturnal boundary layers, which is known to be particularly difficult to simulate and, at the same time, to have a large impact on near-surface concentrations. While the assimilation of hourly data provided improved results in the idealized experiment, the preliminary result for real measurements revealed unexpectedly large differences from the results obtained with the assimilation of daytime data only. These real-world tests with hourly data raised concerns about the inclusion of nighttime observations in particular, which are very challenging for the model to represent correctly. Accurate inversion depends on the ensemble spread reliably capturing uncertainties and the model being free of significant nighttime biases. Given the need for further investigation of these issues, we have limited our demonstration to the afternoon/nighttime averages to ensure more-robust conclusions.

In conclusion, our findings demonstrate the advantages of integrating temporally varying, flow-dependent atmospheric transport uncertainties into inversions to enhance the accuracy of GHG flux estimations. Incorporating these uncertainties yields more-accurate estimates of GHG emissions, with significant improvements across a wide range of setups.