Last modified on Jan 08, 2025 14:55

Change log

Contributor(s)

1. Introduction

The CAMS flux inversions of CH₄ are based on the TM5-4DVar inverse modelling system (Bergamaschi et al., 2010, 2013; Meirink et al., 2008). The inversion system estimates monthly global fields of CH₄ surface fluxes that, when provided as input to the TM5 chemistry-transport model, provide the best possible match between simulated and observed CH₄ concentrations. The various components and configuration of this system are described below. Notations follow Table 1 from Rayner, Michalak, & Chevallier (2019).

1.1. State vector

The inversion system for CH₄ fluxes estimates the optimal value of a state vector x that consists of emission deviation factors, and eventually other parameters such as initial concentrations. An overview of the possible elements of the state is given in Table 1; the associated properties are described in following subsections.

1.2. Emission deviation values

The emission deviations are the main elements of the state vector: this is the target output of the whole system. The deviations are monthly time series of global 2D fields with shape (longitude,latitude). In the current configuration, deviation factors are included for 4 methane source category groups:

• wetlands;
• rice fields;
• biomass-burning;
• other sources (mainly anthropogenic).

An element

for grid cell (i,j) and time t quantifies a deviation from the background emissions

The actual emission

is computed from:

This semi-exponential formulation is needed to ensure that actual emissions are always positive, since

could also become negative to enforce a decrease of emissions to obtain a decrease of concentrations.

1.3. Initial concentrations

An optional second part of the state vector is a 3D field with initial concentrations. For latest inversions this has not been used, as the initial concentrations are taken from an inversion already and therefore not require further optimization.

1.4. Background state

The background state x^b is the initial estimate. The 'x^b' column in Table 1 lists the chosen values.
For the emission deviation part, the background values are simply zero, which implies no deviation from the emission inventories.
A background or first guess for the initial concentrations could be taken from an available model run, preferably the result of a previous 4D-Var inversion, such that concentrations are already in good agreement with observations. However, as described above, for latest inversions initial concentrations are not part of the state.

Table 1: Elements of the 4D-Var state and associated properties.

1.5. Cost function

In the inversion system, the optimal state x is defined as the state where the following cost-function reaches a minimum value:

This cost-function assigns a penalty when x differs from a background state x^b, and when the simulation H(x) differs from the observations y. The observation operator H(x) consists of the chemistry-transport model that simulates concentrations given all kinds of input data (meteorological fields etc.) and the emissions from the state estimate x, and applies a sampling operator to simulate the observations. The background covariance matrix B and the observation representation covariance matrix R define the relative weight of the two penalties in the total sum. The covariance matrices are described in detail in section 1.9 and in the observation descriptions.
The implementation of the elements of the cost function and the minimization procedure are described in the following sections.

1.6. Minimization algorithm

The state x for which the cost function in Eq. (2) reaches a minimum is obtained using the M1QN3 algorithm (Gilbert & Lemaréchal, 1989). The minimum is searched with an iterative procedure that provides in every iteration step a more optimal state. The input for each step is formed by:

• the current state estimate x_i;
• the evaluated costs J(x_i) for this estimate;
• the gradient ∇_xJ(x_i) of the cost function with respect to the elements of the state, evaluated in the current estimate.

The result of an iteration step is a new estimate x_i+1 of the optimal state. The iteration process is stopped when convergence is reached.

1.7. Cost-function gradient

The minimization procedure requires evaluation of the gradient of the cost function towards the elements of the state. The gradient can be computed from:

In here, H^T represents the adjoint of the tangent linear observation operator:

This adjoint operator mainly consists of the adjoint of the chemistry-transport model towards the state variables. The input for the adjoint observation operator is the departure vector:

which consists of the differences between the simulations H(x) and the observations y, weighted with the assumed observation representation error covariance R.
The inverse B^-1 of the background covariance is not explicitly computed; instead, a pre-conditioner is applied as described below.

1.8. Pre-conditioner

The pre-conditioner is a state transformation that is used to avoid computations with the full covariance matrix B and its inverse. The pre-conditioned state is defined as:

with the reverse transformation:

Formulated in terms of the new state vector w, the minimization procedure should be applied on the cost function:

with gradient:

The iteration loop of the optimization starts with w=o, which is the pre-conditioned equivalent of x=x^b.

1.9. Background covariance

The background covariance B statistically describes how the true state x might deviate from the background value x^b. In this application, it describes the uncertainty in emissions (and eventually the initial concentrations). The full B matrix is not computed but is only present in the algorithm as a parameterization; the operations with B (or the inverse B^-1, or square root B^-1/2 of the inverse) are implemented using these parameterizations.
The uncertainties in the emission deviation and the other parameters in the state (for example initial concentrations) are assumed to be uncorrelated. The B matrix has therefore a block-diagonal form:

where in this example B_e is the background covariance for the emission deviations, and B_c is the background covariance for the initial concentrations.

1.9.1. Background covariance for emission deviations

The emission deviations are assumed to be uncorrelated between the four categories; their covariance is therefore block-diagonal:

Each of the four diagonal blocks describes the covariance in the monthly time series of 2D emission deviation fields. This covariance is parameterized as the product of a standard deviation and a correlation:

where the S matrix is a diagonal matrix formed from a standard deviation field, and the correlation matrix C holds the spatial and temporal correlations. The assumed standard deviations are listed in the 'sigma' column of Table 1. A standard deviation of 1.0 (100%) is assumed, except for category 'other' for which a standard deviation of 0.50 (50%) is assumed.
The correlations within the time series of emission deviation fields are assumed to be separable in a spatial (horizontal) and temporal component, and can therefore be written as a Kronecker product:

The horizontal correlations are parameterized as being homogeneous (the same everywhere) and isotropic (independent of direction). The correlation between the emission deviations in two grid cells therefore only depends on the distance between the grid cells; here a Gaussian decay is assumed:

where d_i,j is the distance between the centers of cell i and j in km, and L is the correlation length scale. Table 1 shows the chosen length scales, which are in this application 500 km for each of the emission categories.
Temporal correlations are omitted (C_t=I) for the wetlands, rice, and biomass burning categories. Thus, for these categories, large fluctuations between the months are allowed, which is in line with the strongly fluctuating seasonal pattern in the corresponding emissions. For category 'other' (mainly anthropogenic), the deviations are assumed to be smoother in time, and therefore a temporal correlation is prescribed with an exponential decay that depends on the distance in months:

where m and n are the month numbers, and τ is the temporal correlation scale of 9.5 months (Table 1).

1.9.2. Background covariance for initial concentrations

The background covariance for the initial concentrations (if used) is assumed to be separable in a horizontal and vertical component. It is written as a Kronecker product between a horizontal correlation and a vertical covariance:

Similar as for the emissions, the horizontal correlations are assumed to be homogeneous and isotropic, and are parameterized using a Gaussian decay: 

where d_i,j is the distance between the centers of cell i and j, and L=500 km is the correlation length scale (Table 1).
The vertical covariance is shown in Figure 1. This covariance was constructed in an earlier study (Meirink et al., 2006) by sampling a covariance from the difference between two model runs driven by different meteorological inputs (NMC method). Note that this vertical covariance is defined on the 25 vertical layers used for the coarse resolution version of the TM5 model; since for this release a model version on 34 layers is used, the vertical covariance has been projected on the new layer definition.

Figure 1: Vertical covariance prescribed for the background covariance of the initial concentrations. Taken from Fig 1 in (Meirink, Eskes, and Goede 2006).

1.9.3. Eigenvalue decomposition

The operations involving the background covariance are actually not using the matrix B itself, but only the square root B^1/2 and its transpose B^T/2. These are implemented using an eigenvalue decomposition of the correlation matrices; in a general notation:

where Λ is the diagonal matrix with (strictly positive) eigenvalues. In terms of these matrices, the transformation from pre-conditioned to normal state then becomes:

and the pre-conditioned gradient:

The expression of the covariances in block-diagonal matrices and Kronecker products remains for the eigenvalue decomposition. It is therefore sufficient to decompose only the individual horizontal and temporal correlation matrices, and the vertical covariance. 

1.10. Observations

The observation vector y collects all measured values that should be used to optimize the state. The observations are collected for the entire time window of the 4D-Var optimization.
In the current application, the observations are ground-based observations from the NOAA network, and satellite column-mean mixing ratios from satellites.

1.11. Observation simulations

The observation operator H(x) simulates the observations in y given a state vector x.
The core element of H(x) is the chemistry transport model TM5-MP, which is here used to simulate methane concentrations given an estimate of the actual emissions, initial concentrations, and other input data such as meteorology and methane sink rates. The properties and configuration of the model are described in more detail in 4.
The other element of H(x) is the algorithm that computes an equivalent of an observation from the model fields. Typically, this consists of spatial interpolations from the model grid to the observation locations, aggregation of concentrations into total columns, and averaging in time.
The observation representation error R describes the uncertainty in the observation representation. It should quantify how much the simulations in H(x) are expected to differ from the observations y. Differences are typically due to instrumental errors and model representation errors that arise from the limited horizontal and vertical resolutions. The exact formulation of the observation representation errors is left for the description of the observation data.

2. EDGAR v8.0 emission processing

The EDGAR v8.0 inventory¹ for Greenhouse Gasses (Crippa et al., 2024) covers the period 1970-2022 when production of the v23r2 release started (mid 2024). In the course of the project, an update named EDGAR_2024_GHG became available, which also includes 2023.

From this inventory, the source category 'AGS (Agricultural Soils)' is used for super-category 'rice' in the inversions, while the other inventory categories are used in the 'other' super-category (which has mainly the anthropogenic emission).

2.1. Monthly time profiles

In the EDGAR v8.0 inventory, monthly files were available for the years 2000-2022. Initially some inconsistencies on the AGS (Agricultural Soils) category were found, which were corrected by the EDGAR team. These corrections have also been included in the EDGAR_2024_GHG release.
The time series of monthly files was extended to the full inversion period 1983-2024 (including spin down) by using the monthly profile for the year 2000 for the earlier years, and the profile of 2022 for recent years.

2.2. Continuation for recent years

The time series of EDGAR emissions has been extended to 2023-2024 following the same method as was used for previous releases, using latest available information. As basis, yearly total production and/or usage numbers for fossil fuels have been obtained from statistics collected by British Petroleum (BP), as well as yearly numbers of agricultural production obtained from statistics collected by the United Nations Food and Agriculture Organization (FAO). Depending on the source category, the relative change in one of these statistics has been used as proxy for the growth of emissions with respect to 2021. For most recent years, when these statistics are not yet available, a linear trend is estimated from the previously collected data.

The continuation of emissions is illustrated in Figure 2 for two different source categories. For the ENF (Enteric Fermation) category, the proxy derived from cattle number shows a steady increase in time, and this is therefore continued in the emission estimate for recent years. For the PRO (fuel exploitation) category however, a clear irregularity is seen in the in the fossil fuel production proxies for 2020. This reflects the decrease in production related to the outbreak of global pandemic in that year, which led at least temporary to lower production anticipating on lower consumption. For 2022-2023, the global emission estimate is however assumed to be increased with respect to 2021 based on the proxies.

Figure 2: Examples of continuation of EDGAR emission timeseries after 2021 for Enteric Fermentation (top) and Fuel Exploitation (bottom). The yearly total EDGAR emissions (blue) are continued using a proxy (red) that is derived from FAO or BP statistics. When these statistics are available the relative change with respect to 2021 is used (solid red markers); for most recent years, a value is estimated using linear extrapolation (open red markers).

3. Atmospheric loss rates

The destruction of CH₄ by other molecules in the atmosphere is the most important sink for methane. The only other substantial removal process is uptake by soils, but this amount is small compared to the chemical loss. Therefore, the balance between emission of CH₄ and destruction by atmospheric sinks determines roughly the overall concentration of methane. To be able to estimate methane fluxes with inverse modelling, it is necessary to estimate the atmospheric destruction.

3.1. Methane loss reactions

The most import loss reaction is that with OH, which takes place in all layers of the atmosphere with a temperature (T) dependent reaction rate (k):

CH₄ + OH → ... , k = 2.45 ∙ 10^-12 e^-1775/T

In the stratosphere also the following loss reactions are important:

CH₄ + O(¹D) → ... , k = 1.5625 ∙ 10^-10
CH₄ + Cl^- → … , k = 7.3 ∙ 10^-12 e^-1280/T

To simulate the chemical loss of methane it is therefore necessary to have estimates of the concentrations of OH in the troposphere, and OH, O(¹D) and Cl^- in the stratosphere. This chapter describes how timeseries of 3D concentrations have been constructed. 

3.2. Simulated time series of sinks concentrations

A timeseries of sinks concentrations has been defined based on model simulations performed within the scope of the CAMS projects. Since none of the available simulations covered the entire timeseries, a harmonized set was created by imposing trends on the available partly time series.
The model simulations that were used are:

• The CAMS Re-analysis (Inness et al., 2019), for the current inversion available from 2003 up to June 2021. This product does not contain any COVID lock-down emission scenarios for 2020 or 2021.
• Dedicated simulations with IFS/CB05/BASCOE (V Huijnen et al., 2016) performed for 1990-2009. For the stratospheric sinks, these runs provided a climatology.
• New simulations for 1998-1999 using an updated version of IFS/CB05/BASCOE, which had lower CO production by bio-mass burning. This increased the tropospheric OH concentrations over regions covered by tropical forest, since those concentrations are believed to be too low in the previous simulations. The spatial patterns of tropospheric OH are obtained from these runs.
• Simulations for 2020 with IFS/CB05 for the "CONFORM" project, which compared the impact of emission changes following a global COVID19 lock-down scenario with reference simulations. An evaluation of the IFS model simulations (Vincent Huijnen et al., 2021) showed that the overall tropospheric OH concentration decreased which could be attributed to the lower emissions of pollutants. However, at places with high NO_x emissions the OH surface concentrations actually increased as a result of the shifted ratio between NOx and volatile organic compounds such as formaldehyde, since emissions of the later were not reduced in the scenarios.

3.3. Harmonized time series of sinks

A harmonized time series of monthly tropospheric OH has been obtained with the following steps.

1. Monthly spatial patterns are initially taken from the 1998-1999 IFS/CB05/BASCOE simulations. Each month is therefore an average from two years.
2. The monthly spatial patterns are combined with a yearly time series of averaged tropospheric concentrations to impose inter-annual variation. This timeseries is derived using the following steps.
   1. For 2003-2020, initially the time series of yearly averaged tropospheric OH concentrations in the CAMS Reanalysis is used. Note that these simulations did not include COVID lock-down scenarios for the emissions.
   2. For 2021-2023, the values for 2020 are re-used.
   3. A linear trend for 1990-2002 has been derived from the first (long) series of IFS/CB05/BASCOE simulations. This trend was then adjusted to match with the trend over 2003-2009 in the CAMS Reanalysis. The inter-annual variation within the 1990-2002 period was kept after adjusting the linear trend.
   4. The linear trend used for the 1990's was extrapolated backwards in time to 1979; see also section 3.5.
3. For 2020 the spatial patterns are scaled following the (gridded) ratio between tropospheric OH from the IFS run with COVID lock-down emission scenario and the reference run; see section 3.6.
4. The total time series 1979-2022 was increased with 10% using a calibration factor obtained from simulations with a methyl-chloroform tracer, as will be described in section 3.4.

For the stratospheric concentrations of OH, O(¹D), and Cl^-, a climatology has been created based on the 1990-2009 series of IFS/CB05/BASCOE simulations.

3.4. Calibration using methyl-chloroform simulations

To validate the new timeseries of sink fields, a simulation of methyl-chloroform (MCF) concentrations has been performed. After the ban on production of MCF in the 1990's following the Montreal Protocol, the concentrations of MCF decayed almost exclusively due to destruction by OH. Therefore, simulations of MCF could be used to investigate whether estimated OH concentrations roughly lead to the observed decay of MCF.
During the preparation of the previous v20r1 release of the inversion, the TM5-MP model was used to simulate MCF over 1998-2019 in order to validate the tropospheric OH concentrations. For the current v21r1 release, neither the model nor the OH concentrations changed for these years, and therefore it was not necessary to repeat the experiment, and the results could be re-used.
The model was configured for MCF simulations following (Naus et al., 2020). The same MCF emissions were used (which are rather small due to the ban on production), as well as the deposition velocities to the ocean, and the destruction rates due to photolysis by ultra-violet light in the stratosphere. The simulation over 1998-2019 started from optimized initial concentrations that were also provided with the input fields. The OH concentrations (tropospheric and stratospheric) created for the methane inversions were used to simulate the chemical degradation.
Figure 3 shows a time series of monthly MCF observations and simulations, averaged over all available surface station data. The observations show a clear decay in time, which is simulated well by a simulation using the created OH sinks (red line). The concentrations are slightly too high, and therefore the OH concentrations were increased using multiplication factor applied to the entire series. The best match (using a single calibration factor) was found for an increase by roughly 10%. The experiments suggested that further tuning could be beneficial, with an even stronger increase during the first years, and a smaller increase later on. This exercise is left for future inversions, when also new simulations of OH and other sinks will be available.

Figure 3: Timeseries of monthly average MCF observations and simulations using time series of OH concentrations.

3.5. Backwards extension to 1983

The tropospheric OH concentrations obtained from the CAMS simulations shows a linear increasing trend for the 1990's. This linear trend has been extrapolated backwards in time to 1983 to obtain sink fields for the entire inversion window. The assumption that a linear trend for 1983-2000 could be used is based on simulations with global atmospheric chemistry models described in (D. S. Stevenson et al., 2020), from which an illustration has been provided here as Figure 4.

Figure 4: Timeseries of global annual mean tropospheric OH taken from Figure 2.(a) of (D. S. Stevenson et al., 2020). The red ellipses mark the years 1980-2000 which are of interest for this inversion study. Although inter-annual variation is resent, an overall liner increasing trend is visible for these years.

3.6. Impact of COVID-19 lockdowns in 2020

The effect of the COVID-19 lockdowns in 2020 is incorporated in the (tropospheric) OH sink field. The lockdown correction is based on the CAMS report (Huijnen et al., 2021), which investigates the changes in primary pollutants and secondary compounds in response to modified emissions during the COVID-19 pandemic using IFS/CB05 simulations. The simulations show an overall lower OH concentration in the troposphere due to reduced NO_x emissions.
In order to produce a consistent and continuous OH field compared to the years before the pandemic, the OH fields available from the IFS/CB05 simulations are not used as input for the inversion directly. Instead, the original sink field obtained for 2020 has been scaled with the relative changes in the OH field due to the COVID-19 pandemic following the IFS/CB05 simulations. The scaling is done on a monthly resolution and layer by layer and leads to a decrease in the global average of the OH sink field for 2020. For the year 2021, the same sink field as in 2019 is used without further scaling or imposing trends.
Figure 5 shows the relative changes in surface concentrations between the original OH field and the new version, which includes the scaling for the lockdown effect. Overall OH concentrations could decrease up to 40%, in particular over South America, North America, southern Europe and Africa. However, due to interaction between ozone, solar radiation and formaldehyde, local increased levels of OH are possible, for example over Great Brittan and parts of eastern China.
The monthly scaling of the timeseries of tropospheric OH is further illustrated in Figure 6. The original time series the years shown (blue) was based on IFS/CB05 simulations from the CAMS Reanalysis (Inness et al., 2019) and application of the calibration factor derived from methyl-chloroform simulations. The new timeseries is based on this previous time series but for the year 2020 it has been scaled with the relative changes of the OH field as modelled by (Huijnen et al., 2021). Compared to the reference series a drop in OH concentrations of 0.002 ppt is introduced in March 2020, which corresponds to a decrease of 5%. This difference between the time series vanishes throughout the year.

Figure 5: Relative difference between original OH concentration at the surface (v20) and the values obtained from the simulations including the effect of the COVID-19 pandemic (v21).  

Figure 6: Timeseries of monthly averaged global tropospheric OH concentrations for the years 2019 and 2020. Average concentrations for COVID-19 scenario (v21) are shown in cyan and original concentration (v20) are shown in blue. 

4. Simulation model

The global atmospheric transport model used in the inversions is TM5 (V. Huijnen et al., 2010; Krol et al., 2005). For the production of the v20r1 inversions, the model was updated to "TM5-MP", the "massive parallel" version. The TM5-MP version is able to employ parallel computing more efficiently, but does not significantly change the simulations.
The components of the model relevant for the production chain are described in the following sections.

4.1. Horizontal resolution

For the current release, the model is configured on a global 1^o×1^o horizontal grid. Near the poles, the number of cells is reduced to avoid presence of small grid cells and associated requirement of small time steps.

4.2. Vertical layers

Simulations are performed at 34 vertical layers, defined as a coarsening of ECMWF's 137-layer definition used for the ERA-5 meteorological input. Figure 7 shows the distribution of the original layers and the coarsening. Especially in the stratosphere and the lower troposphere the number of layers is reduced.

Figure 7: Vertical layers in meteorological input (ERA-5,L137) and TM5 34 layer coarsening.

4.3. Meteorological input

The meteorological input for the model is taken from ECMWF's ERA5 re-analysis (Hersbach et al., 2020). The original data has a resolution of ~30 km, which is processed into TM5 input fields on 1^ox1^o degrees. The original 137 layers are coarsened to the model layers as described in section 4.2. The temporal resolution of the data processed for TM5 is hourly for the surface fields, and 3-hourly for 3D fields; at run time, linear interpolation in time is used to obtain the values at the model time step. The configuration of the meteorological input is summarized in Table 2.

Table 2: Configuration of meteorological input.

4.4. Emissions of CH₄

The emissions of CH₄ are read as gridded timeseries. The emissions selected for this application are described in more detail in the Evaluation and Quality Control document of the release. The processing of the EDGAR anthropogenic emissions is outlined in detail in Chapter 2 of this document.

4.5. Chemical loss of CH₄

Methane is chemically destroyed in the atmosphere. Apart from a small soil sink, this is the only mechanism that removes methane. The most import loss reactions are that with OH in the troposphere, and with OH, O(¹D) and Cl^- in the stratosphere. 3 describes how time series of these concentrations have been created for the current inversion.
To simulate the chemical destruction of methane, the TM5-MP model configured for this study uses monthly mean concentrations of sink tracers. For each tracer, the loss rate (1/s) is computed using the reaction rates listed in section 3.1 in combination with the monthly averaged temperature. The rates are assigned to either the troposphere or the stratosphere using the tropopause parameterization of (Lawrence et al., 2001).

4.6. Interpretation of mole mixing ratios

To compare the model simulations with observations it is necessary to know how to interpret simulated mixing ratios with respect to dry or wet air. To ensure mass conservation, TM5 uses a constant global total air mass. This should be interpreted as a dry air mass, since the actual air mass is not constant but fluctuates with the water load.
Mole mixing ratios of methane could therefore be interpreted as dry air mole fractions. As also NOAA methane observations are reported as dry air mole fractions (gml.noaa.gov/ccgg/trends_ch4), the mole fractions from the model and the observations can be compared directly without conversions based on humidity.

5. Physical Parallelization for long inversion windows

To limit processing time, an algorithm named "Physical Parallelization" (PP) is used to process the long inversion window of more than 40 years in parallel chunks.

5.1. Algorithm

The PP algorithm is described in detail in (Pandey et al., 2022). The approach is to solve the 4D-var optimization over the entire multi-decadal window, but to perform model simulations in parallel over shorter windows. For each of the shorter windows, an approximation of the impact of earlier emissions is taken into account using a series of box models (Figure 9). The simulations from the box models are combined with a sensitivity run and added to the full model simulations as a bias correction.

Figure 9: Simulation of long time windows using "Physical Parallelization".

The actual algorithm is illustrated in Figure 10. For latest release v23r2, it has been implemented using the following steps.

1. The entire time window is defined to be the target time window 1983-2023, with in addition a 5 months spin-down period in order to use information from recent observations to the optimize emissions at the end of the target window. A spin-up period is not considered as necessary since the window is already initialized from optimized concentrations.
2. The 41-year period 1979-2020 is divided in 13 blocks of 3 years, and a 14^th block of 2 years. Except for the first block, each block is preceded by 6 extra months to serve as overlap between blocks. The overlap is used to compute sensitives of observation-minus-simulation differences in the following months towards emission changes in the overlap period. The 14^th block also contains the 5-month spin-down period.
3. An initial concentration field co for the start of the entire window is obtained from the v22r2 release of the inversion that covered 1979-1922.
4. A standard model simulation is performed for the entire window, starting from c₀, and using a priori emissions x^a. This model simulation should be interpreted a single serial run, but is implemented as a sequence of simulations over the blocks, from the second block onwards initialized from concentrations saved from the previous block.
5. For each block k a sensitivity simulation is performed with initial concentrations of 1 ppb in the entire atmosphere, and zero emissions. These runs could be performed in parallel over the blocks. The output consist of simulations of the observations that will be used for the inversion: h_k=H_k(c=1,x=0)
6. The 4D-var optimization cycle is performed over the entire window. Each iteration provides an update xⁱ of the emissions using two model simulations:
   1. A forward simulation over the entire window is gathered from simulations performed in parallel over the blocks:
      1. The simulation for a block k is initialized from an a priori concentration field cka taken from the model simulation (step 4), and uses the latest emission estimate xⁱ. The output consists of simulated observations:
         $$m^i_k = \mathbf{H}_k ( \mathbf{c}^a_k, x^i)$$
      2. Since the model simulation starts from the a priori concentration field cka, these simulations do not take into account the change in estimated emissions for previous blocks. A correction is computed based on the difference between latest emission estimate and a priori emissions, and the output of the sensitivity runs (see section 6.2 for full description of all the entities):
         $$\mathbf{n}_k^i = \mathbf{h}_k \sum_{l=1}^{k-1} S_{l,k} \mathbf{F} \left(x^i_l - x^a_l\right)$$
      3. The simulated observations for the entire time window are gathered from the simulations per block and the corresponding corrections:
         $$\mathbf{m}^i = [\dots, \mathbf{m}_k^i + \mathbf{n}_k^i, \dots] $$
   2. An adjoint simulation over the entire window is gathered from simulations performed in parallel over the blocks:
      1. An adjoint forcing is computed per block using the simulated observations and corrections, the observations y_k, and observation representation errors R_k that in this application also depend on information from the model simulation of step 4:
         $$f^i_k = \mathbf{R}_k^{-1} (m^i_k + n^i_k - y) $$
      2. A first contribution to the gradient of the 4D-var cost-function towards changes in the emissions is computed by the adjoint simulation over the block:
         $$g^i_k = \mathbf{H}_k^T f^i_k$$
      3. Additional contributions are computed using the adjoint of the correction formula, which account for the sensitivity of the corrections in future blocks for emissions in block k:
         $$\delta g_k^i = \mathbf{F}^\top \sum_{i=1}^{k-1} S_{l,k} \left(\mathbf{h}_k^\top \mathbf{f}_k^i\right)$$
      4. A gradient for the entire window is collected by adding contributions from all blocks:
         $$g^i = \sum_{k=1}^{K} g^i_k$$

   3. A new estimate of the emissions is computed using the gathered simulations and gradients.
7. Finally, an a posterior simulation after n iterations is performed in serial, starting from the initial concentration co, and using the latest optimized emissions xⁿ.

Figure 10: Illustration of PP algorithm, taken from "Figure 1" in (Pandey et al., 2022).

5.2. Computation of simulation corrections

An important element of the PP algorithm is the computation of corrections to the observation simulations for emission changes in previous blocks, as described under step 6.a.ii:

In the current inversion, the emission fields x have a monthly temporal resolution. This resolution is therefore also used by the concentration operator F and the sink operator S_l,k, but to keep notations simple and in line with the original algorithm description, the descriptions below are done in terms of blocks.

5.3. Global emission-to-concentrations operator

The emission-to-concentration operator F is a box model that simulates the global concentration change δx in ppb for an emission change , if this is distributed equally over the entire atmosphere. This is a crude conversion from emissions to immediate well mixed concentrations, while in reality it takes months to years before an emission of any trace gas could be called well mixed over the atmosphere.

For a scalar δx in kg CH₄, the simulation is:

where:

is the number of moles in 1 kg emitted CH₄, with M_CH4 the molar mass [kg/mole], and where:

is the number of moles air in the atmosphere based on the average global surface pressure p_s,g=985 hPa, the global area A_g computed using an Earth radius of 6371 km, the gravity constant g, and the molar mass of air M_air in kg/mole.

For the grided emission fields, the total emission change is calculated and the same factor F is applied.

5.4. Sink operator

Over time, the well-mixed background concentration that is left from the F operator is chemically destroyed by atmospheric sinks, mainly OH in the troposphere. The operator S_k,l describes the fraction of methane that is left at the start of block l after emission and global uniform mixing in block k. In (Pandey et al., 2022) a sink rate was computed assuming an exponential decay with a constant atmospheric lifetime τ of 9 years:

As described in Section 4.5, the TM5-MP model simulates the sinks using monthly timeseries of concentration fields, and these are therefore not constant. In this application, the output of the sensitivity simulations is used to compute an effective lifetime as present in the model. This is visible in the simulated surface concentrations h_k at observations sites as computed in step 5 of the algorithm. Figure 11 shows the monthly averages over all these simulations for each of the 11 blocks. Starting from the uniform 1 ppb initialization, the concentrations decrease in time. An average life time for each block is computed by fitting an exponential decay function with the data, and resampling this at monthly resolution; the result is shown in Figure 12. The lifetimes computed in this way are about 11.3 years in 1979, and decrease to about 9.5 years around 2008 and thereafter. This is in line with the OH concentrations used for the atmospheric sinks, that are currently estimated to be increasing in the 1980's and 1990's than later on, leading to a longer life at the start of the time window.

Figure 11: Monthly averaged concentrations at surface station locations as simulated per PP block in the sensitivity simulations.

Figure 12: Effective CH₄ lifetime in TM5-MP model simulations as computed from sensitivity runs.

5.5. Example timeseries

As illustration of the PP method, timeseries of observed and simulated CH₄ at station Mauna Loa are shown in Figure 13 for the first 4 blocks as were used in the v20r1 release (including a 6-month spin-up that was omitted in later releases). The model simulations (black) form a continuous time series, starting from a posteriori concentrations that match the observations at the begin of the spin-up period. The simulated CH₄ growth rate seems to be too high, since the model is more and more over-estimating the observations. The inversion tries to correct for this by adjusting the emissions; a net decrease of emissions is expected, although locally also increments could be present. A new emission estimate is obtained in every iteration of the 4D-var cycle. The cyan lines show the forward simulations mki for iteration i=3; it is based on the same initial concentration as used for the model simulation, but uses the latest estimate of the emission; as expected, the net emission reduction leads to lower simulations. The blue lines are the simulated concentrations after correction for previous emission changes, thus mki+nki for i=3. For the first block, the correction is zero, and the corrected simulations are equal to the forward simulation. The results show that already after 3 iterations the system is able to provide simulations close to the observations, also for the blocks at the end of the time window (not shown here).

Figure 13: Observed and simulated concentrations at site Mauna Loa during the first 4 blocks of the PP algorithm. Black line is the standard model simulation; cyan is simulation in iteration 3, blue line is the same after correction for previous emission changes.

Support

All queries about the TM5-MP/4D-VAR system or its products can be addessed through https://confluence.ecmwf.int/site/support.

References

Bergamaschi, P., Houweling, S., Segers, A., Krol, M., Frankenberg, C., Scheepmaker, R. A., Dlugokencky, E., Wofsy, S. C., Kort, E. A., Sweeney, C., Schuck, T., Brenninkmeijer, C., Chen, H., Beck, V., & Gerbig, C. (2013). Atmospheric CH 4 in the first decade of the 21st century: Inverse modeling analysis using SCIAMACHY satellite retrievals and NOAA surface measurements. Journal of Geophysical Research: Atmospheres, 118(13), 7350–7369. https://doi.org/10.1002/jgrd.50480

Bergamaschi, P., Krol, M., Meirink, J. F., Dentener, F., Segers, A., van Aardenne, J., Monni, S., Vermeulen, A. T., Schmidt, M., Ramonet, M., Yver, C., Meinhardt, F., Nisbet, E. G., Fisher, R. E., O’Doherty, S., & Dlugokencky, E. J. (2010). Inverse modeling of European CH4 emissions 2001-2006. Journal of Geophysical Research, 115(D22), D22309. https://doi.org/10.1029/2010JD014180

Crippa, M., Guizzardi, D., Pagani, F., Schiavina, M., Melchiorri, M., Pisoni, E., Graziosi, F., Muntean, M., Maes, J., Dijkstra, L., Van Damme, M., Clarisse, L., and Coheur, P.: Insights into the spatial distribution of global, national, and subnational greenhouse gas emissions in the Emissions Database for Global Atmospheric Research (EDGAR v8.0), Earth Syst. Sci. Data, 16, 2811–2830, doi:10.5194/essd-16-2811-2024, 2024.

Gilbert, J. C., & Lemaréchal, C. (1989). Some numerical experiments with variable-storage quasi-Newton algorithms. Mathematical Programming, 45(1), 407–435. https://doi.org/10.1007/BF01589113

Hersbach, H., Bell, B., Berrisford, P., Hirahara, S., Horányi, A., Muñoz-Sabater, J., Nicolas, J., Peubey, C., Radu, R., Schepers, D., Simmons, A., Soci, C., Abdalla, S., Abellan, X., Balsamo, G., Bechtold, P., Biavati, G., Bidlot, J., Bonavita, M., … Thépaut, J.-N. (2020). The ERA5 global reanalysis. Quarterly Journal of the Royal Meteorological Society, 146(730), 1999–2049. https://doi.org/https://doi.org/10.1002/qj.3803

Huijnen, V., Williams, J., van Weele, M., van Noije, T., Krol, M., Dentener, F., Segers, A., Houweling, S., Peters, W., de Laat, J., Boersma, F., Bergamaschi, P., van Velthoven, P., Le Sager, P., Eskes, H., Alkemade, F., Scheele, R., Nédélec, P., & Pätz, H.-W. (2010). The global chemistry transport model TM5: description and evaluation of the tropospheric chemistry version 3.0. Geoscientific Model Development, 3(2), 445–473. https://doi.org/10.5194/gmd-3-445-2010

Huijnen, V, Flemming, J., Chabrillat, S., Errera, Q., Christophe, Y., Blechschmidt, A.-M., Richter, A., & Eskes, H. (2016). C-IFS-CB05-BASCOE: stratospheric chemistry in the Integrated Forecasting System of ECMWF. Geosci. Model Dev., 9(9), 3071–3091. https://doi.org/10.5194/gmd-9-3071-2016

Huijnen, Vincent, Bouarar, I., Brasseur, G., Remy, S., Schreurs, T., & Williams, J. (2021). Global trace gas changes in the IFS mini-ensemble in response to changes in primary emissions during the COVID-19 pandemic. https://doi.org/10.24380/4g5p-gj91

Inness, A., Ades, M., Agustí-Panareda, A., Barré, J., Benedictow, A., Blechschmidt, A.-M., Dominguez, J. J., Engelen, R., Eskes, H., Flemming, J., Huijnen, V., Jones, L., Kipling, Z., Massart, S., Parrington, M., Peuch, V.-H., Razinger, M., Remy, S., Schulz, M., & Suttie, M. (2019). The CAMS reanalysis of atmospheric composition. Atmos. Chem. Phys., 19(6), 3515–3556. https://doi.org/10.5194/acp-19-3515-2019

Krol, M., Houweling, S., Bregman, B., van den Broek, M., Segers, A., van Velthoven, P., Peters, W., Dentener, F., & Bergamaschi, P. (2005). The two-way nested global chemistry-transport zoom model TM5: algorithm and applications. Atmospheric Chemistry and Physics, 5, 417–432.

Lawrence, M. G., Jöckel, P., & von Kuhlmann, R. (2001). What does the global mean OH concentration tell us? Atmos. Chem. Phys., 1(1), 37–49. https://doi.org/10.5194/acp-1-37-2001

Meirink, J. F., Bergamaschi, P., & Krol, M. C. (2008). Four-dimensional variational data assimilation for inverse modelling of atmospheric methane emissions: method and comparison with synthesis inversion. Atmospheric Chemistry and Physics, 8(21), 6341–6353. https://doi.org/10.5194/acp-8-6341-2008

Meirink, J. F., Eskes, H. J., & Goede, A. P. H. (2006). Sensitivity analysis of methane emissions derived from SCIAMACHY observations through inverse modelling. Atmospheric Chemistry and Physics, 6(5), 1275–1292. https://doi.org/10.5194/acp-6-1275-2006

Naus, S., Montzka, S. A., Patra, P. K., & Krol, M. C. (2020). A 3D-model inversion of methyl chloroform to constrain the atmospheric oxidative capacity. Atmos. Chem. Phys. Discuss., 2020, 1–23. https://doi.org/10.5194/acp-2020-624

Pandey, S., Houweling, S., & Segers, A. (2022). Order of magnitude wall time improvement of variational methane inversions by physical parallelization: a demonstration using TM5-4DVAR. Geosci. Model Dev., 15(11), 4555–4567. https://doi.org/10.5194/gmd-15-4555-2022

Rayner, P., Michalak, A. M., & Chevallier, F. (2016). Fundamentals of Data Assimilation. Geoscientific Model Development Discussions, 1–21. https://doi.org/10.5194/gmd-2016-148

Stevenson, D. S., Zhao, A., Naik, V., O’Connor, F. M., Tilmes, S., Zeng, G., Murray, L. T., Collins, W. J., Griffiths, P. T., Shim, S., Horowitz, L. W., Sentman, L. T., & Emmons, L. (2020). Trends in global tropospheric hydroxyl radical and methane lifetime since 1850 from AerChemMIP. Atmos. Chem. Phys., 20(21), 12905–12920. https://doi.org/10.5194/acp-20-12905-2020

6. Related articles

• Page:
  FTP access to CAMS global data
• Page:
  FTP access to CAMS GFAS data
• Page:
  How to reference or acknowledge data from CAMS (Copernicus Atmosphere Monitoring Service) in a publication
• Page:
  CAMS Regional: European air quality reanalyses data documentation
• Page:
  CAMS Regional: European air quality analysis and forecast data documentation