Articles | Volume 14, issue 5
Model description paper
26 May 2021
Model description paper |  | 26 May 2021

Variational regional inverse modeling of reactive species emissions with PYVAR-CHIMERE-v2019

Audrey Fortems-Cheiney, Isabelle Pison, Grégoire Broquet, Gaëlle Dufour, Antoine Berchet, Elise Potier, Adriana Coman, Guillaume Siour, and Lorenzo Costantino

Up-to-date and accurate emission inventories for air pollutants are essential for understanding their role in the formation of tropospheric ozone and particulate matter at various temporal scales, for anticipating pollution peaks and for identifying the key drivers that could help mitigate their concentrations. This paper describes the Bayesian variational inverse system PYVAR-CHIMERE, which is now adapted to the inversion of reactive species. Complementarily with bottom-up inventories, this system aims at updating and improving the knowledge on the high spatiotemporal variability of emissions of air pollutants and their precursors. The system is designed to use any type of observations, such as satellite observations or surface station measurements. The potential of PYVAR-CHIMERE is illustrated with inversions of both carbon monoxide (CO) and nitrogen oxides (NOx) emissions in Europe, using the MOPITT and OMI satellite observations, respectively. In these cases, local increments on CO emissions can reach more than +50 %, with increases located mainly over central and eastern Europe, except in the south of Poland, and decreases located over Spain and Portugal. The illustrative cases for NOx emissions also lead to large local increments (> 50 %), for example over industrial areas (e.g., over the Po Valley) and over the Netherlands. The good behavior of the inversion is shown through statistics on the concentrations: the mean bias, RMSE, standard deviation, and correlation between the simulated and observed concentrations. For CO, the mean bias is reduced by about 27 % when using the posterior emissions, the RMSE and the standard deviation are reduced by about 50 %, and the correlation is strongly improved (0.74 when using the posterior emissions against 0.02); for NOx, the mean bias is reduced by about 24 % and the RMSE and the standard deviation are reduced by about 7 %, but the correlation is not improved. We reported strong non-linear relationships between NOx emissions and satellite NO2 columns, now requiring a fully comprehensive scientific study.

1 Introduction

The degradation of air quality is a worldwide environmental problem: 91 % of the world's population have breathed polluted air in 2016 according to the World Health Organization (WHO), resulting in 4.2 million premature deaths every year (WHO, 2016). The recent study of Lelieveld et al. (2019) even suggests that the health impacts attributable to outdoor air pollution are substantially higher than previously assumed (with 790 000 premature deaths in the 28 countries of the European Union against the previously estimated 500 000; EEA, 2018). The main regulated primary (i.e., directly emitted in the atmosphere) anthropogenic air pollutants are carbon monoxide (CO), nitrogen oxides (NOx= NO + NO2), sulfur dioxide (SO2), ammonia (NH3), volatile organic compounds (VOCs) and primary particles. These primary air pollutants are precursors of secondary (i.e., produced in the atmosphere through chemical reactions) pollutants such as ozone (O3) and particulate matter (PM), which are also threatening to both human health and ecosystems. Monitoring concentrations and quantifying emissions are still challenging and limit our capability to forecast air quality to warn population and to assess (i) the exposure of population to air pollution and (ii) the efficiency of mitigation policies.

Bottom-up (BU) inventories are built in the framework of air quality policies such as The Convention on Long-Range Transboundary Air Pollution (LRTAP,, last access: March 2019) for air pollutants. Based on national annual inventories, research institutes compile gridded global or regional monthly inventories (mainly for the US, Europe and China) with a high spatial resolution (currently regional- or city-scale inventories are typically finer than 0.1×0.1). These inventories are constructed by combining available (economic) statistics data from different detailed activity sectors with the most appropriate emission factors (defined as the average emission rate of a given species for a given source or process, relative to the unit of activity in a given administrative area). It is important to note that the activity data (often statistical data) have an inherent uncertainty and that their reliability may vary between countries or regions. In addition, the emission factors bear large uncertainties in their quantification (Kuenen et al., 2014; EMEP/EEA, 2016; Kurokawa et al., 2013). Moreover, these inventories are often provided at the annual or monthly scale with typical temporal profiles to build the weekly, daily and hourly variability of the emissions. The combination of uncertain activity data, emission factors and emission timing can be a large source of uncertainties, if not errors, for forecasting or analyzing air quality (Menut et al., 2012). Finally, since updating the inventories and gathering the required data for a given year is costly in time, manpower and money, only a few institutes have offered estimates of the gaseous pollutants for each year since 2011 (i.e, European Monitoring and Evaluation Programme EMEP updated until the year 2017, MEIC updated until the year 2017 to our knowledge). Nevertheless, using knowledge from inventories and air quality modeling, emissions have been mitigated. For example, from 2010 to nowadays, emissions in various countries have been modified and/or regional trends have been reversed downwards (e.g., the decrease in NOx emissions over China since 2011; de Foy et al., 2016), leading to significant changes in the atmospheric composition. Consequently, the knowledge of precise and updated budgets, together with seasonal, monthly, weekly and daily variations of gaseous pollutants driven, amongst other processes, by the emissions are essential for understanding their role in the formation of tropospheric ozone and PMs at various temporal scales, for anticipating pollution peaks and for identifying the key drivers that could help mitigate these concentrations.

In this context, complementary methods have been developed for estimating emissions using atmospheric observations. They operate in synergy between a chemistry-transport model (CTM) which links the emissions to the atmospheric concentrations, atmospheric observations of the species of interest and statistical inversion techniques. A number of studies using inverse modeling were first carried out for long-lived species such as greenhouses gases (GHGs) (e.g., carbon dioxide CO2 or methane CH4) at the global or continental scales (Hein et al., 1997; Bousquet et al., 1999), using surface measurements. Later, following the development of monitoring station networks, the progress of computing power and the use of inversion techniques more appropriate to non-linear problems, these methods were applied to shorter-lived molecules such as CO. For these various applications (e.g., for CO2, CH4, CO), the quantification of sources was solved at the resolution of large regions (Pétron et al., 2002). Finally, the growing availability and reliability of observations since the early 2000s (in situ surface data, remote sensing data such as satellite data) and the improvement of the global CTMs, computational capacities and inversion techniques have increased the achievable resolution of global inversions, up to the global transport model grid cells, i.e., typically with a spatial resolution of several hundreds of square kilometers (Stavrakou and Müller, 2006; Pison et al., 2009; Fortems-Cheiney et al., 2011; Hooghiemstra et al., 2012; Yin et al., 2015; Miyazaki et al., 2017; Zheng et al., 2019).

Today, the scientific and societal issues require an up-to-date quantification of pollutant emissions at a higher spatial resolution than the global one, which will lead to a wide use of regional inverse systems. However, although they are suited to reactive species such as CO and NOx, and their very large spatial and temporal variability, they have hardly been used to quantify pollutant emissions. Some studies inferred NOx (Pison et al., 2007; Tang et al., 2013) and VOC emissions (Koohkan et al., 2013) from surface measurements. Konovalov et al. (2006, 2008, 2010), Mijling et al. (2012, 2013), van der A et al. (2008), Lin et al. (2012) and Ding et al. (2017) have also shown that satellite observations are a suitable source of information to constrain NOx emissions. These regional inversions using satellite observations were often based on Kalman filter (KF) schemes (Mijling et al., 2012, 2013; van der A et al., 2008; Lin et al., 2012; Ding et al., 2017).

Variational inversion systems allow for the solving of high-dimensional problems, typically solving for the fluxes at high spatial and temporal resolution, which can be critical to fully exploit satellite images. Here, we present the Bayesian variational atmospheric inversion system PYVAR-CHIMERE for the monitoring of anthropogenic emissions of reactive species at the regional scale. It is based on the Bayesian variational assimilation code PYVAR (Chevallier et al., 2005) and on the regional state-of-the-art CTM CHIMERE (Menut et al., 2013; Mailler et al., 2017). CHIMERE is dedicated to the study of regional atmospheric pollution events (e.g., Ciarelli et al., 2019; Menut et al., 2020), included in the operational ensemble of the Copernicus Atmosphere Monitoring Service (CAMS) regional services. The main strengths of PYVAR-CHIMERE come from the strengths of CHIMERE and from its high modularity for the definition of the control vector. CHIMERE is indeed an extremely flexible code, in particular for the definition of the chemical scheme.

The PYVAR-CHIMERE system takes advantage of the previous developments for the quantification of fluxes of long-lived GHG species such as CO2 (Broquet et al., 2011) and CH4 (Pison et al., 2018) at the regional to the local scales, but it now solves for reactive species such as CO and NOx. It has also a better level of robustness, clarity, portability and modularity than these previous systems. Variational techniques require the adjoint of the model to compute the sensitivity of simulated atmospheric concentrations to corrections of the fluxes. CHIMERE is one of the few CTMs for which the adjoint has been coded. Global models include GEOS-CHEM (Henze et al., 2007), IMAGES (Stavrakou and Müller, 2006), TM5 (Krol et al., 2008), GELKA (Belikov et al., 2016) and LMDz (Chevallier et al., 2005; Pison et al., 2009); limited-area models include CMAQ (Hakami et al., 2007), EURAD-IM (Elbern et al., 2007), RAMS/CTM-4DVAR (Yumimoto and Uno, 2006) and WRF-CO2 4D-Var (Zheng et al., 2018).

The principle of variational atmospheric inversion and the configuration of PYVAR-CHIMERE are described in Sects. 2 and 3, respectively. Details about the forward, tangent-linear (TL) and adjoint codes of CHIMERE are also given. Then, the potential of PYVAR-CHIMERE is illustrated in Sect. 4 with the optimization of European CO and NOx emissions, constrained by observations from the Measurement of Pollution in the Troposphere (MOPITT) and from the Ozone Monitoring Instrument (OMI) satellite instruments, respectively.

2 Principle of Bayesian variational atmospheric inversion

In what follows, we use the notations and equations used in the inverse modeling community (Rayner et al., 2019). The Bayesian variational atmospheric inversion method adjusts a set of control parameters, including parameters related to the emissions whose estimate is the primary target of the inversion.

The prior information about the parameters x to be optimized during the inversion process is given by the vector xb. The parameters to be optimized can be surface fluxes but may also include initial or boundary conditions for example, as explained in Sect. 3.4. The adjustments are applied to prior values, usually taken, for the emissions, from pre-existing BU inventories. The principle is to minimize, on the one hand, the departures from the prior estimates of the control parameters, which are weighted by the uncertainties in these estimates (called hereafter “prior uncertainties”), and, on the other hand, the differences between simulated and observed concentrations, which are weighted by all other sources of uncertainties explaining these differences (hereafter referred to collectively as “observation errors”). In statistical terms, the inversion searches for the most probable estimate of the control parameters given their prior estimates, observations, CTM and associated uncertainties. The solution, which will be called posterior estimate, is found by the iterative minimization of a cost function J (Talagrand, 1997), defined as follows:

(1) J ( x ) = x - x b T B - 1 x - x b + H ( x ) - y T R - 1 H ( x ) - y .

H is the non-linear observation operator that projects the control vector x onto the observation space. In most of the variational atmospheric inversion cases (such as those described in Sect. 4), the observation operator includes the operations performed by the CTM in linking the emissions to the concentrations and any other transformation to compute the simulated equivalent of the observations such as an interpolation or an extraction and averaging of the simulated concentration fields (see Sect. 3.5). The observations in y could be surface measurements and/or remote sensing data such as satellite data. The prior uncertainties and the observation errors are assumed to be unbiased and to have a Gaussian distribution. Consequently, the prior uncertainties are characterized by their covariance matrix B and the observation errors are characterized by their covariance matrix R. By definition, the observation errors combine errors in both the data and the observation operator, in particular the following:

  • measurement errors and errors in the conversion of satellite measurement into concentration data,

  • errors from the CTM,

  • representativity errors due to the comparison between point measurements and gridded models or due to the representation of the fluxes as gridded maps at a given spatial resolution, and

  • aggregation errors associated with the optimization of emissions at a given spatial and/or temporal resolution (as specified in the control vector) that is different from (usually coarser than) that of the CTM (Wang et al., 2017).

Figure 1Simplified scheme of the iterative minimization in PYVAR-CHIMERE. PYVAR, CHIMERE and text sources are displayed in blue, in orange and in grey, respectively.


For inversions with observation and control vectors with a high dimension, the minimum of J cannot be found analytically due to computational limitations. It can be reached iteratively with a descent algorithm. In this case, the iterative minimization of J is based on a gradient method. J is calculated with the forward observation operator (including the CTM), and its gradient relative to the control parameters x is provided by the adjoint of the observation operator (including the adjoint of the CTM). The gradient is defined as follows:

(2) J ( x ) = B - 1 x - x b + H * R - 1 H ( x ) - y ,

where H* is the adjoint of the observation operator.

The high non-linearity of the chemistry for reactive species makes it difficult to use its TL code to approximate the actual observation operator, and, more generally, it makes the inversion problem highly non-linear. Therefore, in PYVAR-CHIMERE, we use the M1QN3 limited-memory quasi-Newton minimization algorithm (Gilbert and Lemaréchal, 1989), which relies on the actual CHIMERE non-linear model to compute J at each iteration of the minimization. As with most quasi-Newton methods, it requires an initial regularization of x, the vector to be optimized, for better efficiency. We adopt the most generally used regularization, made by minimizing in the space defined by the following:

(3) χ = B 1 2 x - x b ,

instead of the control space defined by x. Although more advanced regularizations can be chosen, the minimization with χ is preferred because it simplifies the equation to solve. In the χ space, Eq. (2) can be re-written as follows:

(4) J χ = χ + B 1 2 H * R - 1 H ( x ) - y .

The criterion for stopping the algorithm is based on a threshold set on the ratio between the final and initial gradient norms or on the maximum number of iterations to perform. As shown in Fig. 1, the minimization algorithm repeats the forward-adjoint cycle to get an estimate close to the optimal solution of the inversion problem for the control parameters. This approximation of the optimal estimate is found by satisfying the convergence criteria of the minimizer with a given reduction of the norm of the gradient of J. Nevertheless, due to the non-linearity of the problem, the minimization may reach a local minimum only, instead of the global minimum.

Finally, the calculation of the uncertainty in the estimate of emissions from the inversion, known as “posterior uncertainty”, is challenging in a variational inverse system (Rayner et al., 2019). Even though the posterior uncertainty can be explicitly written in various analytical forms, it requires the inversion of matrices that are too large to invert given the current computational resources in our variational approach. As a trade-off between computing resources and comprehensiveness, the analysis error may be evaluated by an approach based on a propagation of errors through sensitivity tests (e.g., as in Fortems-Cheiney et al., 2012). It can also be estimated through a Monte Carlo ensemble (Chevallier et al., 2007), implemented in PYVAR. Nevertheless, it should be noted that the cost of the Monte Carlo experiments used to derive these posterior uncertainties is huge.

Figure 2Simplified scheme of how PYVAR scripts are used to drive CHIMERE for an inversion using satellite observations. PYVAR, CHIMERE and text sources are displayed in blue, in orange and in grey, respectively. “AK” refers to averaging kernels as detailed in Sect. 3.5.


3 The PYVAR-CHIMERE configuration

3.1 PYVAR adapted to CHIMERE

The PYVAR-CHIMERE inverse modeling system is based on the Bayesian variational assimilation code PYVAR (Chevallier et al., 2005) and on a previous inversion system coupled to CHIMERE (Pison et al., 2007). PYVAR is an ensemble of Python scripts, which deals with preparing the vectors and the matrices for the inversion, drives the required Fortran codes of the transport model and computes the minimization of the cost function to solve the inversion. Previously used for global inversions with the LMDz model (e.g., Pison et al., 2009; Chevallier et al., 2010; Fortems-Cheiney et al., 2011; Yin et al., 2015; Locatelli et al., 2015; Zheng et al., 2019), PYVAR has been adapted to CHIMERE with an adjoint code without chemistry by Broquet et al. (2011). In order to couple PYVAR to the new state-of-the-art version of CHIMERE (see Sect. 3.2), to include chemistry and to increase its modularity, flexibility and clarity, the new system described here has been developed. It includes elements of the inversion system (coded in Fortran90) of Pison et al. (2007).

3.2 Development and parallelization of the adjoint and tangent-linear codes of CHIMERE

To compute the sensitivity of simulated atmospheric concentrations to corrections to the fluxes, the adjoint of CHIMERE has been developed. Originally, the sequential adjoint was coded (Menut et al., 2000; Menut, 2003; Pison et al., 2007). The adjoint has been coded by hand line by line, following the principles formulated by Talagrand (1997). It contains exactly the same processes as the CHIMERE forward model. The code has been parallelized, which required a redesigning of the entire code, associated with a full testing scheme (see Sect. 3.3). Furthermore, the TL code has been developed and validated (see Sect. 3.3). Changes have been implemented in the forward CHIMERE code embedded in PYVAR-CHIMERE to match requirements of the studies conducted with this system. These changes have been implemented in both the adjoint and the TL codes. Compared to the CHIMERE 2013 version (Menut et al., 2013), the most important of these changes are, regarding geometry, the possibility of polar domains and the use of the coordinates of the corners of the cells instead of only the centers, allowing the use of irregular grids. Regarding transport, the non-uniform Van Leer transport scheme on the horizontal has been implemented, which is consistent with the use of irregular grids. Finally, various switches have been added to keep the system consistent for GHG studies. For example, we can avoid going into the chemistry, deposition or wet deposition routines when the targeted species do not require them (e.g., no chemistry for methane or carbon dioxide at a regional scale).

PYVAR-CHIMERE is currently implemented with a full module of gaseous chemistry. As a compromise between the robustness of the method for reactive species, the time required to code the adjoint and the computational cost with a full chemical scheme, the aerosols modules of CHIMERE have not been included in the adjoint of CHIMERE yet and are therefore not available in PYVAR-CHIMERE. The development and maintenance of the adjoint means that the version used is necessarily one or two versions behind the distributed CHIMERE version (, last access: March 2019). It should also be noted that PYVAR-CHIMERE only infers anthropogenic emissions at this stage. The optimization of biogenic emissions, which are linearly interpolated at the sub-hourly scale in CHIMERE, is currently under development.

As an example, Fig. 2 presents a simplified scheme of how PYVAR scripts are used to drive this version of CHIMERE for forward simulations and inversions using satellite observations.

Table 1Examples for the definition of the control vector and for the construction of the B matrix, as illustrated in Sect. 4.

a Add, mult or scale. b Fluxes, initial conditions, lateral boundary conditions or top boundary conditions. c Fixed values (fx) or percentages (%).

Download Print Version | Download XLSX

3.3 Accuracy of tangent-linear and adjoint codes

Different procedures have been implemented to test the accuracy of the TL and adjoint codes. To test the linearity of the TL code, we compute a Taylor diagnostic. It consists of computing the TL code at x0 for given increments Δx, dHx0 (Δx), then the TL code at x0 for λ×Δx with λ an arbitrary small number, dHx0(λΔx). Theoretically, if the TL code is well coded, λdHx0(Δx)=dHx0(λΔx) by definition. In practice, the difference must be lower than 10 times the epsilon of the machine on which it is run.

The adjoint code is also tested, by verifying that dHx0(Δx)|dHx0(Δx)=Δx|H*(dHx0(Δx)) where H* stands for the adjoint at x. What is actually computed is the ratio of the difference between the two scalar products to the second one and the accuracy of the computation. The difference should be a few times greater than the epsilon of the machine on which it is run.

3.4 Definition of the control vector

The control vector is specified by the user in a text file. This file is formatted following Table 1. The parameters to be inverted may be fluxes and/or initial conditions and/or boundary concentration conditions, at the grid-cell resolution or for one region encompassing up to the whole domain. Several types of corrections can be applied; they are defined in the code as “add”, “mult” or “scale”. Both the corrections “add” and “mult” are applied to gridded control variables. For correction type “add”, the control variables are increments added to the corresponding components of the model inputs. For correction type “mult”, the control variables are scaling factors multiplying the corresponding components of the model inputs. The difference between the two options “add” and “mult” plays a role when inverting fluxes which can switch from positive to negative values (like CO2 natural fluxes). For type “scale”, the control variables are scaling factors applied to maps different from the maps of emissions used as prior input of the forward model: for example, activity maps can be used and scaled to get emissions; the obtained values are then added to the corresponding components of the model inputs. With these various types, it is possible to define the control variables as the budgets of emissions for different regions, types of activities and/or processes, which can thus be directly rescaled by the inversions, similarly to what is done in systems where the control vector is not gridded (Wang et al., 2018).

Different simple but efficient ways of building the error covariance matrix B are implemented in PYVAR-CHIMERE. The variances and correlations are defined independently. The variances are specified by the user through the specification of the values for the corresponding standard deviation (i.e., the diagonal matrix of standard deviations ; Table 1) which can be made in terms of fixed values (“fx” in the code) or percentages (%, “pc” in the code). For correction types “mult” and “scale”, as well as for correction type “add” with a fixed value, the value is directly used as the uncertainty in the corresponding components of the control vector. For correction type “add” with a percentage provided, maps of standard deviation of uncertainty are built by applying this percentage to the matching input fields (fluxes, initial conditions, boundary conditions). The user may also provide a script to build personalized maps of variances.

Potential correlations between uncertainties in different types of control variables (e.g., between fluxes and boundary conditions) and correlations between uncertainties in different species (e.g., between fluxes of CO and NOx) are not coded yet. Such correlations increase the observation constraint on the emissions in the inversion process by transferring information from one species to the other. The level (and sometimes the sign) and thus the impact on the inversion of such correlations highly depend on the study cases and are often debated due to the lack of precise characterization of the uncertainties in inventories of anthropogenic emissions of GHG and pollutants (Super et al., 2020). Only correlations for a given type of control variable and a given species are so far taken into account so that the B matrix is block diagonal. For a given type of control variable and a given species (in the illustration in Sect. 4.2.2: CO, NO or NO2 fluxes), spatial and temporal correlations can be defined using correlation lengths through time LT and space LS. Those lengths are used to model temporal and/or spatial auto-correlations using an exponentially decaying function: the correlation r between parameters at a given location but separated by duration d(xi,xj) or at a given time but distant by d(xi,xj) is given by

(5) r x i , x j = exp - d x i , x j L ,

where L is the corresponding correlation length. There is no correlation between uncertainties in land and ocean flux. Note that the spatial correlations are computed for each vertical level independently when dealing with control variables with vertical resolution (3D fields of fluxes when accounting for emission injection heights, or boundary/initial conditions). Vertical correlations in the uncertainties in such variables have not been coded yet. Apart from this, the system assumes that temporal correlations and spatial correlations depend on the time lag and distance but not on the specific time and location of the corresponding parameters. It also assumes that the correlation between uncertainties at different locations and different times can be derived from the product of the corresponding autocorrelation in time and space.

Each block of B can thus be decomposed based on Kronecker products:

(6) B = C T C S ,

where is the Kronecker product, and CT and CS are the temporal and spatial correlations, respectively. The calculations involving B1/2 (in Eqs. 3, 4) are simplified in PYVAR-CHIMERE using the eigendecomposition of CT and CS. Its square root can be calculated according to the following:

(7) C T 1 / 2 = V C T D C T 1 / 2 V C T T

(and similarly for CS), where VCT is the matrix with the eigenvectors as columns, and DCT is the diagonal matrix of eigenvalues of CT. It is possible to choose a threshold under which the eigenvalues are truncated when computing the spatial correlations in order to save computation time and memory, but not when computing the temporal correlations.

3.5 Equivalents of the observations

During forward simulations, the equivalents of the components of y (i.e, the equivalents of the individual data) are calculated by PYVAR-CHIMERE. It includes the CTM and an interpolation (see below the vertical interpolation from the model's grid to the satellite levels) or an extraction and averaging (e.g., extracting the grid cell matching the geographical coordinates of a surface station and averaging over 1 h). As a compromise between technical issues such as the time required for reading and writing files, the observation operator H that generates the equivalent of the observations by the model (i.e., H(x)) has been so far partly embedded in the code of CHIMERE. It makes it easier to use finer time intervals than available in the usual hourly outputs of CHIMERE to compute the required information (e.g., within the finer CTM physical time steps).

To make comparisons between simulations and satellite observations, the simulated vertical profiles are first interpolated on the satellite's levels (with a vertical interpolation on pressure levels) in CHIMERE. Then, the averaging kernels (AKs), when available, are applied to represent the vertical sensitivity of the satellite retrieval. Two types of formula, depending on the satellite observations used, have been detailed in PYVAR-CHIMERE for the use of AKs:

(8) c m = A K c m ( o )


(9) c m = x a + A K c m ( o ) - x a ,

where cm is the modeled column, AK contains the averaging kernels that can be provided in the form of a vector (e.g., OMI product) or matrix (e.g., MOPITT product), xa is the prior state vector (provided together with the AKs when relevant), and cm(o) is the vertical distribution of the original model partial columns interpolated to the pressure grid of the AKs.

3.6 Numerical language

The PYVAR code is in Python 2.7, the CHIMERE CTM is coded in Fortran90. The CTM requires several numerical tools, compilers and libraries. The PYVAR-CHIMERE system was developed and tested using the software versions as described in Table 2.

Table 2URL addresses for the development and the use of the PYVAR-CHIMERE system and its modules.

Download Print Version | Download XLSX

PYVAR-CHIMERE's computation time for one node of 10 CPUs is about 4 h for 1 d of inversion (with  10 iterations) for the European domain size of 101 (longitude) × 85 (latitude) × 17 (vertical levels) used in Sect. 4. As described in Menut et al. (2013) for CHIMERE, the model parallelization results from a Cartesian division of the main geographical domain into several sub-domains, each one being processed by a worker process. To configure the parallel sub-domains, the user has to specify two parameters in the model parameter file: the number of sub-domains for the zonal and meridian directions. The total number of CPUs used is therefore the product of these two numbers plus one for the master process. The optimal number of CPUs for the parallelization of the transport scheme depends on the size of the tiles and also of the technical characteristics of the machine, because of the time required to exchange halos.

4 Potential of PYVAR-CHIMERE for the inversion of CO and NOx emissions

The potential of the PYVAR-CHIMERE system to invert emissions of reactive species is illustrated with the inversion of CO and NOx anthropogenic emissions in Europe based on MOPITT CO data and OMI NO2 data, respectively. We have chosen to present an illustration of CO inversion over 7 d, the first week of March 2015. Considering the short lifetime of NOx of a few hours (Valin et al., 2013; Liu et al., 2016), we have chosen to present illustration of NOx inversion over 1 d, 19 February 2015. These particular periods have been chosen as they present a representative number of super-observations during winter, and as the emissions are high during that period. All the information required by the system to invert CO and NOx emissions is listed in Table 1.

4.1 Data and model description

4.1.1 Observations y

We use CO data from the MOPITT instrument (Deeter et al., 2019). MOPITT has been flown onboard the NASA EOS-Terra satellite, on a low sun-synchronous orbit that crosses the Equator at 10:30 and 22:30 LST. The spatial resolution of its observations is about 22×22 km2 at nadir. It has been operated nearly continuously since March 2000. MOPITT CO products are available in three variants: thermal-infrared TIR only, near-infrared NIR only and the multispectral TIR-NIR product, all containing total columns and retrieved profiles (expressed on a 10-level grid from the surface to 100 hPa). We choose to constrain CO emissions with the MOPITT surface product for our illustration. Among the different MOPITTv8 products, we choose to work with the multispectral MOPITTv8-NIR-TIR one, as it provides the highest number of observations, with a good evaluation against in situ data from NOAA stations (Deeter et al., 2019). The MOPITTv8-NIR-TIR surface concentrations are sub-sampled into “super-observations” in order to reduce the effect of errors that are correlated between neighboring observations: we selected the median of each subset of MOPITT data within each 0.5×0.5 grid cell and each physical time step (about 5–10 min). After this screening, 8437 “super-observations” remain in the 7 d inversion (from 10 667 raw observations). It is important to note that the potential of MOPITT to provide information at a high temporal resolution, up to the daily scale, is hampered by the cloud coverage (see the blanks in Fig. 5b).

The observational constraint on NO2 emissions comes from the OMI QA4ECV tropospheric columns (Muller et al., 2016; Boersma et al., 2016, 2017). The Ozone Monitoring Instrument (OMI), a near-UV–visible nadir solar backscatter spectrometer, was launched onboard EOS Aura in July 2004. It has been flown on a 705 km sun-synchronous orbit that crosses the Equator at 13:30 LT. Our data selection follows the criteria of the OMI QA4ECV data quality statement. As the spatial resolution of the OMI data is finer than this of the chosen CHIMERE model grid (13×24 km2 against 0.5×0.5, respectively), the OMI tropospheric columns are sub-sampled into “super-observations” (median of the OMI data within the 0.5×0.5 grid cell and each physical time step and its corresponding AKs).

Figure 3Mean CO surface concentrations from 1 to 7 March 2015 simulated by CHIMERE (a) with anthropogenic and biogenic emissions, and (b) without emissions, in ppbv, at the 0.5× 0.5 grid-cell resolution.

Figure 4Mean CO surface concentrations from 1 to 7 March 2015 simulated by CHIMERE using for initial and boundary conditions, (a) the climatological values from the LMDZ-INCA global model; (b) the climatological values from a MACC reanalysis, in ppbv; and (c) the relative differences between these two simulations , in %, at the 0.5× 0.5 grid-cell resolution.

4.1.2 CHIMERE setup

CHIMERE is run over a 0.5×0.5 regular grid (about 50×50 km2) and 17 vertical layers, from the surface to 200 hPa (about 12 km), with 8 layers within the first two kilometers. The domain includes 101 (longitude) × 85 (latitude) grid cells (31.5–74 N, 15.5 W–35 E; see Fig. 3). CHIMERE is driven by the European Centre for Medium-Range Weather Forecasts (ECMWF) meteorological forecast (Owens and Hewson, 2018). The chemical scheme used in PYVAR-CHIMERE is MELCHIOR-2, with more than 100 reactions (Lattuati, 1997; CHIMERE,, last 8 June 2017) including 24 for inorganic chemistry. The prior anthropogenic emissions for CO and NOx emissions are obtained from the TNO-GCHco-v1 inventory (Super et al., 2020), the last update of the TNO-MACCII inventory (Kuenen et al., 2014). This inventory is based on the EMEP Centre on Emission Inventories and Projections (CEIP) official country reporting for air pollutants done in 2017. It is an inventory at 6×6 km2 horizontal resolution. From the annual and national budgets, each sector is assigned to a specific proxy to quantify the spatial variability of the emissions within each country. Temporal profiles are also provided per gridded nomenclature for reporting (GNFR) sector code (variations due to the month, weekday and hour). Following the Generation of European Emission Data for Episodes (GENEMIS) recommendations (Kurtenbach et al., 2001; Aumont et al., 2003), NOx emissions are speciated as 90 % of NO, 9.2 % of NO2 and 0.8 % of nitrous acid (HONO). The TNO-GHGco-v1 inventory has been aggregated to the CHIMERE grid.

The prior anthropogenic emissions for VOCs are obtained from the EMEP inventory (Vestreng et al., 2005; EMEP/CEIP:, last access: March 2019). Biogenic emissions come from the Model of Emissions of Gases and Aerosols from nature (MEGAN) (Guenther et al., 2006). Different climatological values from the LMDZ-INCA global model (Szopa et al., 2008) or from a Monitoring Atmospheric Composition and Climate (MACC) reanalysis are used to prescribe concentrations at the lateral and top boundaries and the initial atmospheric composition in the domain. Full access to and more information about the MACC reanalysis data can be obtained through the MACC-II web site (, last access: March 2019). In order to ensure realistic fields of simulated CO and NO2 concentrations from the beginning of the inversion period, runs have been preceded with a 10 d spinup.

Figure 5Mean CO collocated surface concentrations from 1 to 7 March 2015 (a) simulated by CHIMERE using the prior TNO-GHGco-v1 emissions and the climatological values from the LMDZ-INCA global model for initial and boundary conditions, (b) observed by MOPITTv8-NIR-TIR and (c) simulated by CHIMERE using the posterior emissions, in ppbv, at the 0.5× 0.5 grid-cell resolution. Relative differences between MOPITT and (d) the prior CHIMERE simulation or (e) the posterior CHIMERE simulation, in %. Statistics for the comparison between simulations and observations are given in Table 4 for the area in the purple box.

4.1.3 CO sensitivity to emissions and to initial and boundary conditions

With its lifetime of about 2 months, CO could be strongly influenced by the initial and lateral boundary conditions prescribed in the CTM. In fact, as seen in Fig. 3b, initial and boundary conditions provide a relatively flat background and the patterns which appear clearly over the background are linked to surface emissions (Fig. 3a). To characterize the uncertainties in the concentration fields due to the initial and lateral boundary conditions, we performed a sensitivity test by using either climatological values from LMDZ-INCA (Fig. 4a) or a MACC reanalysis (Fig. 4b): maximum relative differences in concentrations of about 15 % over continental land are estimated (Fig. 4c). The errors assigned to initial and boundary conditions in Sect. 4.2.2 are based on this sensitivity test.

4.1.4 Comparison between CHIMERE and the observations

Large discrepancies (Fig. 5d) are found between the MOPITT CO observations (Fig. 5b) and the prior simulation by CHIMERE over Europe (Fig. 5a). For the first week of March 2015, CO concentrations are generally under-estimated by CHIMERE, particularly over central and eastern Europe (except in the south of Poland). On the contrary, CO concentrations seem to be over-estimated over Spain and Portugal. Large discrepancies are also found between the OMI NO2 super-observations and the prior simulation by PYVAR-CHIMERE (Fig. 6d), as already noticed by Huijnen et al. (2010), with an inter-comparison of NO2 OMI-DOMINO tropospheric columns with an ensemble of European regional air quality models including CHIMERE. Over Europe, the prior simulation strongly underestimates the tropospheric columns over industrial areas (e.g., over the Netherlands and Po Valley). These discrepancies might be due to different causes, which can all interact. A source of uncertainties is related to the observations. For example, satellite data inter-comparison studies reveal large differences between different retrievals of the same compound (Qu et al., 2020). This can be explained by uncertainties from the CTM (e.g., through the underestimation of the atmospheric production or the underestimation of the species lifetime). It could also be explained by an underestimation of the anthropogenic emissions in the BU inventory.

Figure 6NO2 collocated tropospheric columns (a) simulated by CHIMERE using the prior TNO-GHGco-v1 emissions and the climatological values from the LMDZ-INCA global model for initial and boundary conditions, (b) observed by OMI and (c) simulated by CHIMERE using the posterior emissions, in 1016 molec. cm−2, at the 0.5× 0.5 grid-cell resolution, for 19 February 2015. Relative differences between OMI and (d) the prior CHIMERE simulation or (e) the posterior CHIMERE simulation, in %. Statistics for the comparison between simulations and observations are given in Table 5 for the area in the purple box.

Table 3Description of the different sensitivity tests performed for the construction of the B matrix for the NOx inversion.

Download Print Version | Download XLSX

4.2 Inversions

4.2.1 Control vector x

For the CO inversion, the control vector x is as follows:

  • the CO anthropogenic emissions at a 7 d temporal resolution, at a 0.5×0.5 (longitude × latitude) horizontal resolution and over the first 8 vertical levels, i.e., for each of the corresponding 101×85×8 grid cells;

  • the CO lateral and top boundary conditions at a 7 d temporal resolution, at a 0.5×0.5 (longitude × latitude) resolution and over the 17 vertical levels of CHIMERE, i.e., (2×101+2×85)×17 grid cells;

  • the CO 3D initial conditions for 1 March 2015 at 00:00 UTC, at a 0.5×0.5 (longitude × latitude) resolution and over the 17 vertical levels of CHIMERE.

Considering its short lifetime, there are no boundary conditions for NO2. For the NOx inversion, the control vector x is as follows:

  • the NO and NO2 anthropogenic emissions at a 1 d temporal resolution, at a 0.5×0.5 (longitude × latitude) resolution, and over the first 8 vertical levels, i.e., for each of the corresponding 101×85×8 grid cells;

  • the NO and NO2 3D initial conditions for 19 February 2015 at 00:00 UTC, at a 0.5×0.5 (longitude × latitude) resolution, and over the 17 vertical levels of CHIMERE.

Figure 7(a) TNO-GHGco-v1 CO anthropogenic prior emissions, in kt CO per grid cell, and (b) increments provided by the inversion with constraints from MOPITTv8-NIR-TIR from 1 to 7 March 2015, in %.

4.2.2 Covariance matrices B and R

To our knowledge, there are few available studies dealing with the estimates of the uncertainties in gridded bottom-up emission inventories at the 0.5× 0.5 resolution or higher. The characterization of their statistics in the inversion configuration is consequently often based on crude assumptions from the inverse modelers. Defining the covariance matrices B and R is not an easy task, while incorrectly specifying these matrices has a very strong impact on the results of the inversion. In particular, the relative weights of B and R and the spatial and temporal correlations in B influence the degree of freedom and the structure for the adjustments attempted by the inversion in the optimization process. Consequently, as an example for the NOx inversion, different sensitivity tests described in Table 3 have been performed for the construction of the B matrix. For both the prior NO and NO2 emissions at 1 d and 0.5 resolution, the prior error standard deviations are first assigned to 50 % of the prior estimate of the emissions (test A), as in Souri et al. (2020). Sensitivity tests have also been performed with prior error standard deviations assigned to 80 % and 100 % of the prior estimate of the emissions (test C and test D, respectively; Fig. 8).

Figure 8(a) TNO-GHGco-v1 NOx anthropogenic prior emissions, in kt NO2 per grid cell and increments provided by the inversion (b) A, (c) B, (d) C, (e) D and (f) E with constraints from OMI 19 February 2015, in %. The description of the different inversions is given in Table 3.

With prior error standard deviations set at 15 % of the initial conditions, the changes in initial conditions are very small (not shown) and do not affect the posterior emissions (test B; Fig. 8). As indicated in Sect. 3.4 and in Table 1, it is possible to use correlations in B, as in Broquet et al. (2011, 2013) and in Kadygrov et al. (2015). We demonstrate the strong impact of spatial correlations, defined by an e-folding length of 50 km over land and over the sea, on our inversion results (test E; Fig. 8).

Table 4Statistics for the comparison between simulated and observed CO surface concentrations over central and eastern Europe (see the area in purple in Fig. 5). MB = mean bias, RMSE = root mean square error, SD = standard deviation are in ppbv. The spatial correlations r are presented with their p value.

Download Print Version | Download XLSX

Even though annual CO emissions in western Europe may be well known, with uncertainties of 6 % according to Super et al. (2020), larger uncertainties could affect eastern Europe. Moreover, large uncertainties still affect bottom-up emission inventories at the 0.5 resolution: spatial disaggregation of the national-scale estimates to provide gridded estimates causes a significant increase in the uncertainty for CO (Super et al., 2020). For the inversion of CO emissions, the error standard deviations assigned to the prior CO emissions at 7 d and 0.5 resolution are 100 %. This value of 100 % has already been chosen in Fortems-Cheiney et al. (2011) and in Fortems-Cheiney et al. (2012). For this CO illustration, the covariance matrix B of the prior errors is defined as diagonal (i.e., only variances in the individual control variables listed in Sect. 4.2.1 are taken into account). With such a setup, in theory, we could obtain negative posterior emissions since the inversion system does not impose a constraint of positivity in the results. Nevertheless, even an uncertainty of 100 % leads to a prior distribution mostly (> 80 %) on the positive side. The assimilation of data showing an increase above the background (at the edges of the domain; not shown) further drives the inversion towards positive emissions for both CO and NOx inversions. In practice, our inversion does not lead to negative posterior emissions (Fig. 7b). Spatial and temporal correlations in B would further limit the probability of getting negative emissions locally by smoothing the posterior emissions at a spatial scale at which the “aggregated” prior uncertainty is smaller than 100 %. However, a positivity constraint should be implemented in future versions of the system.

Based on the sensitivity test in Fig. 4, the errors assigned to the CO lateral boundary conditions and to their initial conditions are set at 15 %. As these relative errors are significantly lower than those for the emissions and as variations in the CO surface concentrations are mainly driven by emissions (Fig. 3), we assume a small relative influence of the correction of initial and boundary conditions on our results. The variance of the individual observation errors in R is defined as the quadratic sum of the measurement error reported in the MOPITT and the OMI data sets, and of the CTM errors (including chemistry and transport errors and representativity errors) set at 20 % of the retrieval values. The representativity errors could have been reduced with the choice of a finer CTM resolution (e.g., with a resolution closer to the size of the satellite pixel). Error correlations between the super-observations are neglected, so that the covariance matrix R of the observation errors is diagonal.

4.2.3 Inversion of CO emissions

A total of 10 iterations are needed to reduce the norm of the gradient of J by 90 % with the minimization algorithm M1QN3 and obtain the increments, i.e., the corrections provided by the inversion. The prior CO emissions over Europe for the first week of March 2015 and their increments are shown in Fig. 7. As expected from the large differences (Fig. 5d) between the prior surface concentrations (Fig. 5a) and the MOPITT observations (Fig. 5b), local increments can reach more than +50 % (Fig. 7b). CO emissions are increased over central and eastern Europe, except in the south of Poland. On the contrary, CO emissions are decreased over Spain and Portugal. The analyzed concentrations are the concentrations simulated by CHIMERE with the posterior fluxes: as expected, the optimization of the fluxes improves the fit of the simulated concentrations to the observations (Fig. 5c and e), particularly over central and eastern Europe. Over this area (see the purple box in Fig. 5), the mean bias between the simulation and the observations has been reduced by about 27 % when using the posterior emissions (mean bias of 11.6 ppbv against 15.9 ppbv with the prior emissions; Table 4). The RMSE and the standard deviation have been reduced by about 50 % and the correlation has been strongly improved (0.74 when using the posterior emissions against 0.02).

Figure 9(a) Bias ratio between CHIMERE simulations using the posterior emissions against prior TNO-GHGco-v1 emissions compared to the OMI-QA4ECV-v1.1 observations. All ratios lower than 1, in blue, demonstrate that posterior emissions improve the simulation compared to the prior ones. (b) OMI uncertainties, in %, for 19 February 2015.

Table 5Statistics for the comparison between simulated and observed NO2 tropospheric columns for the inversion E, mainly over the Netherlands (see the area in purple in Fig. 6). MB = mean bias, RMSE = root mean square error, SD = standard deviation are in molec. cm−2. The spatial correlations r are presented with their p value.

Download Print Version | Download XLSX

4.2.4 Inversion of NOx emissions

The prior NOx emissions and the corrections provided by the different sensitivity tests of Table 3 are shown in Fig. 8. Here, we analyzed the results from inversion E. As expected from the underestimation of the prior tropospheric columns in Fig. 6a, local increments may be large, for example over industrial areas (e.g., over the Po Valley) and over the Netherlands, with increments of more than +50 % (Fig. 8b). The analyzed NO2 tropospheric columns in Fig. 6c are the columns simulated by CHIMERE with the NO2 posterior fluxes: as expected, the optimization of the fluxes improves the fit of the simulated concentrations to the observations over the Netherlands (Fig. 6e). Over this area (see the purple box in Fig. 6), where the OMI uncertainties are lower than 50 % (Fig. 9b), the mean bias between the simulation and the observations has been reduced by about 24 % when using the posterior emissions (mean bias of 1.9 ×1015 molec. cm−2 against 2.6 ×1015 molec. cm−2 with the prior emissions; Table 5, Fig. 9a). The RMSE and the standard deviation have been reduced by about 7 %. The correlation has not been improved.

Even with high emission increments, the impact on the tropospheric columns is rather small. We have performed a test to explain this lack of sensitivity. We have simulated NO2 columns with anthropogenic emissions increased by a factor of 3 compared to the simulation in Fig. 6a. The ratio between these two simulations shows strong non-linearities, blurring the multiplicative effect of our increments and explaining the lack of sensitivity (not shown). By increasing NOx anthropogenic emissions, NO2 tropospheric columns can be strongly increased and can even exceed the observation values for particular pixels. NO2 tropospheric columns can also be decreased or only slightly increased. On average, it tends to increase the concentrations by a factor that is much smaller than the factor of increase in the anthropogenic emissions. However, the patterns where the posterior tropospheric columns exceed the observations or, on the contrary, are decreased or only slightly increased, explain why the inversion system does not attempt to increase further the average level of the concentration (to decrease further the general bias to the observations), even though it accounts for the impact of non-linearities in the chemistry through the use of the M1QN3 minimization algorithm. We can conclude that the strong non-linearities of the NOx chemistry mainly explain the lack of sensitivity between NOx emissions and satellite NO2 columns. Besides chemical effects, the lack of sensitivity could be also partly due to the contribution of emissions during the preceding days, and the assimilation window will be widened in the near future.

The posterior emissions and their uncertainties will have to be evaluated and may give hints as to the cause of the discrepancies between simulated and observed NO2 tropospheric columns. The biases between OMI and simulated NO2 tropospheric columns are a complex topic that is not only related to our CHIMERE simulations (Huijnen et al., 2010; Souri et al., 2020; Elguindi et al., 2020). Several studies have indeed already reported that strong non-linear relationships exist between NOx emissions and satellite NO2 columns (Lamsal et al., 2011; Vinken et al., 2014; Miyazaki et al., 2017; Li and Wang, 2019). This reveals that a fully comprehensive scientific study is required, by analyzing the NOx lifetime through processes such as the NO2+ OH reactions and/or the reactive uptake of NO2 and N2O5 by aerosols (e.g., Lin et al., 2012; Stavrakou et al., 2013).

5 Conclusion and discussion

This paper presents the Bayesian variational inverse system PYVAR-CHIMERE, which has been adapted to the inversion of reactive species such as CO and NOx, taking advantage of the previous developments for long-lived species such as CO2 (Broquet et al., 2011) and CH4 (Pison et al., 2018). We show the potential of PYVAR-CHIMERE, with inversions for CO and NOx illustrated over Europe. PYVAR-CHIMERE will now be used to infer CO and NOx emissions over long periods, e.g., first for a whole season or year and then for the recent decade 2005–2015 in the framework of the H2020 VERIFY project over Europe, and in the framework of the ANR PolEASIA project over China, to quantify their trend and their spatiotemporal variability. Nevertheless, as we have reported strong non-linear relationships between NOx emissions and satellite NO2 columns, a fully comprehensive scientific study is required, by analyzing the NOx lifetime through processes such as the NO2+ OH reactions and/or the reactive uptake of NO2 and N2O5 by aerosols (e.g., Lin et al., 2012; Stavrakou et al., 2013). Biogenic emissions will be also further studied to better understand the relationship between NOx emissions and NO2 spaceborne columns.

The PYVAR-CHIMERE system can handle any large number of both control parameters and observations. It will be able to cope with the dramatic increase in the number of data in the near future with, for example, the high-resolution imaging (pixel of 7×3.5 km2) of the new Sentinel-5P/TROPOMI program, launched in October 2017. These new space missions with high-resolution imaging have the ambition to monitor atmospheric chemical composition for the quantification of anthropogenic emissions. It will indeed entail using PYVAR-CHIMERE at higher spatiotemporal resolutions, but probably for smaller domains (i.e., over countries rather than over Europe) as a compromise between resolution and the computational cost. Moreover, a step forward in the joint assimilation of co-emitted pollutants will be possible with the PYVAR-CHIMERE system and the availability of TROPOMI co-localized images of CO and NO2. This should improve the consistency of the inversion results and can be used to inform inventory compilers and subsequently improve emission inventories. Moreover, this development will help in further understanding air quality problems and addressing air-quality-related emissions at the national to subnational scales.

Code and data availability

The OMI QA4ECV NO2 product can be found here: (last access: March 2019, Boersma et al., 2017).

The MOPITTv8-NIR-TIR CO product can be found here: (last access March 2019, Deeter et al., 2019).

The CHIMERE code is available here: (last access March 2019, CHIMERE, 2019).

The associated documentation of PYVAR-CHIMERE is available on the website (last access March 2019, PYVAR-CHIMERE, 2019). The documentation includes a whole description of PYVAR-CHIMERE and several tutorials on how to run a first PYVAR-CHIMERE simulation or how to run an inversion.

Author contributions

All authors have contributed to the writing the paper (main authors: AFC, GB, IP and GD) and to the development of the present version of the PYVAR-CHIMERE system (main developer: IP). IP and GD have parallelized the adjoint version from Menut et al. (2000, 2003) and Pison et al. (2007). IP has complemented the adjoint of new parameterizations since the CHIMERE release in 2011 and the tangent-linear model.

Competing interests

The authors declare that they have no conflict of interest.


We acknowledge L. Menut and C. Schmechtig for their contributions to the development work on the adjoint code of CHIMERE and its parallelization. We acknowledge the TNO team for providing their inventory of NOx and CO emissions over Europe. We also acknowledge the free use of tropospheric NO2 column data from the OMI sensor and the free use of CO surface concentrations from the MOPITT sensor. A large part of the development and analysis was conducted in the frame of the H2020 VERIFY project, funded by the European Commission Horizon 2020 research and innovation programme. We wish to thank all the persons involved in the preparation, coordination and management of this project. This work was also supported by the CNES (Centre National d’Etudes Spatiales), in the frame of the TOSCA ARGOS project. This work was granted access to the HPC resources of TGCC under the allocations A0050107232 and A0070102201 made by GENCI. Finally, we wish to thank François Marabelle (LSCE) and his team for computer support.

Financial support

For this study, Audrey Fortems-Cheiney was funded by the H2020 VERIFY project, funded by the European Commission Horizon 2020 research and innovation programme, under agreement number 776810. Lorenzo Costantino was funded by the PolEASIA ANR project under the allocation ANR-15-CE04-0005.

Review statement

This paper was edited by Ignacio Pisso and reviewed by three anonymous referees.


Aumont, B., Chervier, F., and Laval, S.: Contribution of HONO sources to the NOx/HOx/O3 chemistry in the polluted boundary layer, Atmos. Environ., 37, 487–498, 2003. 

Belikov, D. A., Maksyutov, S., Yaremchuk, A., Ganshin, A., Kaminski, T., Blessing, S., Sasakawa, M., Gomez-Pelaez, A. J., and Starchenko, A.: Adjoint of the global Eulerian–Lagrangian coupled atmospheric transport model (A-GELCA v1.0): development and validation, Geosci. Model Dev., 9, 749–764,, 2016. 

Boersma, K. F., Vinken, G. C. M., and Eskes, H. J.: Representativeness errors in comparing chemistry transport and chemistry climate models with satellite UV–Vis tropospheric column retrievals, Geosci. Model Dev., 9, 875–898,, 2016. 

Boersma, K. F., Eskes, H., Richter, A., De Smedt, I., Lorente, A., Beirle, S., Van Geffen, J., Peters, E., Van Roozendael, M., and Wagner, T.: QA4ECV NO2 tropospheric and stratospheric vertical column data from OMI (Version 1.1) [Data set], Royal Netherlands Meteorological Institute (KNMI),, 2017. 

Bousquet, P., Ciais, P., Peylin, P., Ramonet, M., and Monfray, P.: Inverse modeling of annual atmospheric CO2 sources and sinks: 1. Method and control inversion, J. Geophys. Res., 104, 26161–26178,, 1999. 

Broquet, G., Chevallier, F., Rayner, P., Aulagnier, C., Pison, I., Ramonet, M., Schmidt, M., Vermeulen, A. T., and Ciais, P.: A European summertime CO2 biogenic flux inversion at mesoscale from continuous in situ mixing ratio measurements, J. Geophys. Res., 116, D23303,, 2011. 

Broquet, G., Chevallier, F., Bréon, F.-M., Kadygrov, N., Alemanno, M., Apadula, F., Hammer, S., Haszpra, L., Meinhardt, F., Morguí, J. A., Necki, J., Piacentino, S., Ramonet, M., Schmidt, M., Thompson, R. L., Vermeulen, A. T., Yver, C., and Ciais, P.: Regional inversion of CO2 ecosystem fluxes from atmospheric measurements: reliability of the uncertainty estimates, Atmos. Chem. Phys., 13, 9039–9056,, 2013. 

Chevallier, F., Fisher, M., Peylin, P., Serrar, S., Bousquet, P., Bréon, F.-M., Chédin, A., and Ciais, P.: Inferring CO2 sources and sinks from satellite observations: method and application to TOVS data, J. Geophys. Res., 110, D24309,, 2005. 

Chevallier, F., Bréon, F.-M., and Rayner, P. J.: The contribution of the Orbiting Carbon Observatory to the estimation of CO2 sources and sinks: Theoretical study in a variational data assimilation framework, J. Geophys. Res., 112, D09307,, 2007. 

Chevallier, F., Ciais, P., Conway, T. J., Aalto, T., Anderson, B. E., Bousquet, P., Brunke, E. G., Ciattaglia, L., Esaki, Y., Fröhlich, M., Gomez, A., Gomez-Pelaez, A. J., Haszpra, L., Krummel, P. B., Langenfelds, R. L., Leuenberger, M., Machida, T., Maignan, F., Matsueda, H., Morguí, J. A., Mukai, H., Nakazawa, T., Peylin, P., Ramonet, M., Rivier, L., Sawa, Y., Schmidt, M., Steele, L. P., Vay, S. A., Vermeulen, A. T., Wofsy, S., and Worthy, D.: CO2 surface fluxes at grid point scale estimated from a global 21 year reanalysis of atmospheric measurements, J. Geophys. Res., 115, 1–17,, 2010. 

CHIMERE: CHIMERE code, available at:, last access: March 2019. 

Ciarelli, G., Theobald, M. R., Vivanco, M. G., Beekmann, M., Aas, W., Andersson, C., Bergstrom, R., Manders-Groot, A., Couvidat, F., Mircea, M., Tsyro, S., Fagerli, H., Mar, K., Raffort, V., Roustan, Y., Pay, M.-T., Schaap, M., Kranenburg, R., Adani, M., Briganti, G., Cappelletti, A., D'Isidoro, M., Cuvelier, C., Cholakian, A., Bessagnet, B., Wind, P., and Colette, A.: Trends of inorganic and organic aerosols and precursor gases in Europe: insights from the EURODELTA multi-model experiment over the 1990-2010 period, Geosci. Model Dev., 12, 4923-4954,, 2019. 

Deeter, M. N., Edwards, D. P., Francis, G. L., Gille, J. C., Mao, D., Martínez-Alonso, S., Worden, H. M., Ziskin, D., and Andreae, M. O.: Radiance-based retrieval bias mitigation for the MOPITT instrument: the version 8 product, Atmos. Meas. Tech., 12, 4561–4580,, 2019. 

de Foy, B., Lu, Z., and Streets, D. G.: Satellite NO2 retrievals suggest China has exceeded its NOx reduction goals from the twelfth Five-Year Plan, Nat. Sci. Rep., 6, 35912,, 2016. 

Ding, J., Miyazaki, K., van der A, R. J., Mijling, B., Kurokawa, J.-I., Cho, S., Janssens-Maenhout, G., Zhang, Q., Liu, F., and Levelt, P. F.: Intercomparison of NOx emission inventories over East Asia, Atmos. Chem. Phys., 17, 10125–10141,, 2017. 

EEA: Air quality in Europe – 2018 report, European Environment Agency, available at: (last access: March 2019), 2018. 

Elbern, H., Strunk, A., Schmidt, H., and Talagrand, O.: Emission rate and chemical state estimation by 4-dimensional variational inversion, Atmos. Chem. Phys., 7, 3749–3769,, 2007. 

Elguindi, N., Granier, C., Stavrakou, T., Darras, S., Bauwens, M., Cao, H., Chen, C., Denier van der Gon, H. A. C., Dubovik, O., Fu, T. M., Henze, D. K., Jiang, Z., Keita, S., Kuenen, J. J. P., Kurokawa, J., Liousse, C., Miyazaki, K., Müller, J.-F., Qu, Z., Solmon, F., and Zheng, B.: Intercomparison of magnitudes and trends in anthropogenic surface emissions from bottom-up inventories, top-down estimates, and emission scenarios, Earth's Future, 8, e2020EF001520,, 2020. 

EEA: EMEP/EEA air pollutant emission inventory guidebook, available at: (last access: 22 April 2021), 2016. 

Fortems‐Cheiney, A., Chevallier, F., Pison, I., Bousquet, P., Szopa, S., Deeter, M. N., and Clerbaux, C.: Ten years of CO emissions as seen from MOPITT, J. Geophys. Res., 116, D05304,, 2011. 

Fortems-Cheiney, A., Chevallier, F., Pison, I., Bousquet, P., Saunois, M., Szopa, S., Cressot, C., Kurosu, T. P., Chance, K., and Fried, A.: The formaldehyde budget as seen by a global-scale multi-constraint and multi-species inversion system, Atmos. Chem. Phys., 12, 6699–6721,, 2012. 

Gilbert, J. and Lemaréchal, C.: Some numerical experiments with variable storage quasi Newton algorithms, Math. Program., 45, 407–435, 1989. 

Guenther, A., Karl, T., Harley, P., Wiedinmyer, C., Palmer, P. I., and Geron, C.: Estimates of global terrestrial isoprene emissions using MEGAN (Model of Emissions of Gases and Aerosols from Nature), Atmos. Chem. Phys., 6, 3181–3210,, 2006. 

Hakami, A., Henze, D. K., Seinfeld, J. H., Singh, K., Sandu, A., Kim, S., Byun, D., and Li, Q.: The adjoint of CMAQ, Environ. Sci. Technol., 41, 7807–7817,, 2007. 

Hein, R., Crutzen, P. J., and Heimann, M.: An inverse modeling approach to investigate the global atmospheric methane cycle, Global. Biogeochem. Cy., 11, 43–76, 1997. 

Henze, D. K., Hakami, A., and Seinfeld, J. H.: Development of the adjoint of GEOS-Chem, Atmos. Chem. Phys., 7, 2413–2433,, 2007. 

Huijnen, V., Eskes, H. J., Poupkou, A., Elbern, H., Boersma, K. F., Foret, G., Sofiev, M., Valdebenito, A., Flemming, J., Stein, O., Gross, A., Robertson, L., D'Isidoro, M., Kioutsioukis, I., Friese, E., Amstrup, B., Bergstrom, R., Strunk, A., Vira, J., Zyryanov, D., Maurizi, A., Melas, D., Peuch, V.-H., and Zerefos, C.: Comparison of OMI NO2 tropospheric columns with an ensemble of global and European regional air quality models, Atmos. Chem. Phys., 10, 3273–3296,, 2010. 

Hooghiemstra, P. B., Krol, M. C., Bergamaschi, P., de Laat, A. T. J., van der Werf, G. R., Novelli, P. C., Deeter, M. N., Aben, I., and Rockmann, T.: Comparing optimized CO emission estimates usingMOPITT or NOAA surface network observations, J. Geophys. Res., 117, D06309,, 2012. 

Kadygrov, N., Broquet, G., Chevallier, F., Rivier, L., Gerbig, C., and Ciais, P.: On the potential of the ICOS atmospheric CO2 measurement network for estimating the biogenic CO2 budget of Europe, Atmos. Chem. Phys., 15, 12765–12787,, 2015. 

Konovalov, I. B., Beekmann, M., Richter, A., and Burrows, J. P.: Inverse modelling of the spatial distribution of NOx emissions on a continental scale using satellite data, Atmos. Chem. Phys., 6, 1747–1770,, 2006. 

Konovalov, I. B., Beekmann, M., Burrows, J. P., and Richter, A.: Satellite measurement based estimates of decadal changes in European nitrogen oxides emissions, Atmos. Chem. Phys., 8, 2623–2641,, 2008. 

Konovalov, I. B., Beekmann, M., Richter, A., Burrows, J. P., and Hilboll, A.: Multi-annual changes of NOx emissions in megacity regions: nonlinear trend analysis of satellite measurement based estimates, Atmos. Chem. Phys., 10, 8481–8498,, 2010. 

Koohkan, M. R., Bocquet, M., Roustan, Y., Kim, Y., and Seigneur, C.: Estimation of volatile organic compound emissions for Europe using data assimilation, Atmos. Chem. Phys., 13, 5887–5905,, 2013. 

Krol, M. C., Meirink, J. F., Bergamaschi, P., Mak, J. E., Lowe, D., Jöckel, P., Houweling, S., and Röckmann, T.: What can 14CO measurements tell us about OH?, Atmos. Chem. Phys., 8, 5033–5044,, 2008. 

Kuenen, J. J. P., Visschedijk, A. J. H., Jozwicka, M., and Denier van der Gon, H. A. C.: TNO-MACC_II emission inventory; a multi-year (2003–2009) consistent high-resolution European emission inventory for air quality modelling, Atmos. Chem. Phys., 14, 10963–10976,, 2014. 

Kurokawa, J., Ohara, T., Morikawa, T., Hanayama, S., Janssens-Maenhout, G., Fukui, T., Kawashima, K., and Akimoto, H.: Emissions of air pollutants and greenhouse gases over Asian regions during 2000–2008: Regional Emission inventory in ASia (REAS) version 2, Atmos. Chem. Phys., 13, 11019–11058,, 2013. 

Kurtenbach, R., Becker, K., Gomes, J., Kleffmann, J., LŽrzer, J., Spittler, M., Wiesen, P., Ackermann, R., Geyer, A., and Platt, U.: Investigations of emissions and heterogeneousformation of HONO in a road traffic tunnel, Atmos. Environ., 35, 3385–3394, 2001. 

Lamsal, L. N., Martin, R. V., Padmanabhan, A., van Donkelaar, A., Zhang, Q., Sioris, C. E., Chance, K., Kurosu, T. P., and Newchurch, M. J.: Application of satellite observations for timely updates to global anthropogenic NOx emission inventories, Geophys. Res. Lett., 38, L05810,, 2011. 

Lattuati, M.: Impact des émissions européennes sur le bilan de l'ozone troposphérique a l'interface de l'europe et de l'atlantique nord: apport de la modélisation lagrangienne et des mesures en altitude, Ph.D. thesis, Université Paris VI, 1997. 

Lelieveld, J., Klingmüller, K., Pozzer, A., Pöschl, U., Fnais, M., Daiber, A., and Münzel, T.: Cardiovascular disease burden from ambient air pollution in Europe reassessed using novel hazard ratio functions, Eur. Heart J., 40, 1590–1596,, 2019. 

Li, J. and Wang, Y.: Inferring the anthropogenic NOx emission trend over the United States during 2003–2017 from satellite observations: was there a flattening of the emission trend after the Great Recession?, Atmos. Chem. Phys., 19, 15339–15352,, 2019. 

Lin, J.-T., Liu, Z., Zhang, Q., Liu, H., Mao, J., and Zhuang, G.: Modeling uncertainties for tropospheric nitrogen dioxide columns affecting satellite-based inverse modeling of nitrogen oxides emissions, Atmos. Chem. Phys., 12, 12255–12275,, 2012. 

Liu, F., Beirle, S., Zhang, Q., Dörner, S., He, K., and Wagner, T.: NOx lifetimes and emissions of cities and power plants in polluted background estimated by satellite observations, Atmos. Chem. Phys., 16, 5283–5298,, 2016. 

Locatelli, R., Bousquet, P., Saunois, M., Chevallier, F., and Cressot, C.: Sensitivity of the recent methane budget to LMDz sub-grid-scale physical parameterizations, Atmos. Chem. Phys., 15, 9765–9780,, 2015. 

Mailler, S., Menut, L., Khvorostyanov, D., Valari, M., Couvidat, F., Siour, G., Turquety, S., Briant, R., Tuccella, P., Bessagnet, B., Colette, A., Létinois, L., Markakis, K., and Meleux, F.: CHIMERE-2017: from urban to hemispheric chemistry-transport modeling, Geosci. Model Dev., 10, 2397–2423,, 2017. 

Menut, L.: Adjoint modelling for atmospheric pollution processes sensitivity at regional scale during the ESQUIF IOP2, J. Geophys. Res.-Atmos., 108, 8562,, 2003. 

Menut, L., Vautard, R., Beekmann, M., and Honoré, C.: Sensitivity of photochemical pollution using the adjoint of a simplified chemistry-transport model, J. Geophys. Res., 105, 15379–15402, 2000. 

Menut, L., Goussebaile, A., Bessagnet, B., Khvorostiyanov, D., and Ung, A.: Impact of realistic hourly emissions profiles on air pollutants concentrations modelled with CHIMERE, Atmos. Envion., 49, 233–244,, 2012. 

Menut, L., Bessagnet, B., Khvorostyanov, D., Beekmann, M., Blond, N., Colette, A., Coll, I., Curci, G., Foret, G., Hodzic, A., Mailler, S., Meleux, F., Monge, J.-L., Pison, I., Siour, G., Turquety, S., Valari, M., Vautard, R., and Vivanco, M. G.: CHIMERE 2013: a model for regional atmospheric composition modelling, Geosci. Model Dev., 6, 981–1028,, 2013. 

Menut, L., Bessanet, B., Siour, G., Mailler, S., Pennel, R., and Cholakian, A.: Impact of lockdown measures to combat Covid-19 on air quality over western Europe, Sci. Total Environ., 140426, 741,, 2020. 

Mijling, B. and van der A, R. J.: Using daily satellite observations to estimate emissions of short-lived air pollutants on a mesoscopic scale, J. Geophys. Res., 117, D17302,, 2012. 

Mijling, B., van der A, R. J., and Zhang, Q.: Regional nitrogen oxides emission trends in East Asia observed from space, Atmos. Chem. Phys., 13, 12003–12012,, 2013. 

Miyazaki, K., Eskes, H., Sudo, K., Boersma, K. F., Bowman, K., and Kanaya, Y.: Decadal changes in global surface NOx emissions from multi-constituent satellite data assimilation, Atmos. Chem. Phys., 17, 807–837,, 2017. 

Muller, J.-P., Kharbouche, S., Gobron, N., Scanlon, T., Govaerts, Y., Danne, O., Schultz, J., Lattanzio, A., Peters, E., De Smedt, I., Beirle, S., Lorente, A., Coheur, P. F., George, M., Wagner, T., Hilboll, A., Richter, A., Van Roozendael, M., and Boersma, K. F.: Recommendations (scientific) on best practices for retrievals for Land and Atmosphere ECVs (QA4ECV Deliverable 4.2 version 1.0), 186 pp., available at: (last access: 12 April 2018), 2016. 

Owens, R. G. and Hewson, T.: ECMWF Forecast User Guide, ECMWF, Reading,, 2018. 

Pétron, G., Granier, C., Khattatov, B., Lamarque, J. F., Yudin, V., Muller, J. F., and Gille, J.: Inverse modeling of carbon monoxide surface emissions using Climate Monitoring and Diagnostics Laboratory network observations, J. Geophys. Res, 107, 4761,, 2002. 

Pison, I., Menut, L., and Bergametti, G.: Inverse modeling of surface NOx anthropogenic emission fluxes in the Paris area during the ESQUIF campaign, J. Geophys. Res.-Atmos., 112, D24302, doi:10,1029/2007JD008871, 2007. 

Pison, I., Bousquet, P., Chevallier, F., Szopa, S., and Hauglustaine, D.: Multi-species inversion of CH4, CO and H2 emissions from surface measurements, Atmos. Chem. Phys., 9, 5281–5297,, 2009. 

Pison, I., Berchet, A., Saunois, M., Bousquet, P., Broquet, G., Conil, S., Delmotte, M., Ganesan, A., Laurent, O., Martin, D., O'Doherty, S., Ramonet, M., Spain, T. G., Vermeulen, A., and Yver Kwok, C.: How a European network may help with estimating methane emissions on the French national scale, Atmos. Chem. Phys., 18, 3779–3798,, 2018. 

PYVAR-CHIMERE: PYVAR-CHIMERE documentation, available at:, last access: March 2019. 

Qu, Z., Henze, D. K., Cooper, O. R., and Neu, J. L.: Impacts of global NOx inversions on NO2 and ozone simulations, Atmos. Chem. Phys., 20, 13109–13130,, 2020. 

Rayner, P. J., Michalak, A. M., and Chevallier, F.: Fundamentals of data assimilation applied to biogeochemistry, Atmos. Chem. Phys., 19, 13911–13932,, 2019. 

Souri, A. H., Nowlan, C. R., González Abad, G., Zhu, L., Blake, D. R., Fried, A., Weinheimer, A. J., Wisthaler, A., Woo, J.-H., Zhang, Q., Chan Miller, C. E., Liu, X., and Chance, K.: An inversion of NOx and non-methane volatile organic compound (NMVOC) emissions using satellite observations during the KORUS-AQ campaign and implications for surface ozone over East Asia, Atmos. Chem. Phys., 20, 9837–9854,, 2020. 

Stavrakou, T. and Müller, J.-F.: Grid-based versus big region approach for inverting CO emissions using Measurement of Pollution in the Troposphere (MOPITT) data, J. Geophys. Res.-Atmos., 111, D15304,, 2006. 

Stavrakou, T., Müller, J.-F., Boersma, K. F., van der A, R. J., Kurokawa, J., Ohara, T., and Zhang, Q.: Key chemical NOx sink uncertainties and how they influence top-down emissions of nitrogen oxides, Atmos. Chem. Phys., 13, 9057–9082,, 2013. 

Super, I., Dellaert, S. N. C., Visschedijk, A. J. H., and Denier van der Gon, H. A. C.: Uncertainty analysis of a European high-resolution emission inventory of CO2 and CO to support inverse modelling and network design, Atmos. Chem. Phys., 20, 1795–1816,, 2020. 

Szopa, S., Foret, G., Menut, L., and Cozic, A.: Impact of large scale circulation on European summer surface ozone: consequences for modeling, Atmos. Environ., 43, 1189–1195,, 2008. 

Talagrand, O.: Assimilation of observations : an introduction, J. Met. Soc. Jpn., 75, 191–209, 1997. 

Tang, X., Zhu, J., Wang, Z. F., Wang, M., Gbaguidi, A., Li, J., Shao, M., Tang, G. Q., and Ji, D. S.: Inversion of CO emissions over Beijing and its surrounding areas with ensemble Kalman filter, Atmos. Environ., 81, 676–686, 2013. 

Valin, L. C., Russell, A. R., and Cohen, R. C.: Variations of OH rad-ical in an urban plume inferred from NO2 column measurements, Geophys. Res. Lett., 40, 1856–1860,, 2013. 

van der A, R. J., Eskes, H. J., Boersma, K. F., van Noije, T. P. C., van Roozendael, M., De Smedt, I., Peters, D. H. M. U., and Meijer, E. W.: Trends, seasonal variability and dominant NOx source derived from a ten year record of NO2 measured from space, J. Geophys. Res., 113, D04302,, 2008. 

Vestreng, V., Breivik, K., Adams, M., Wagner, A., Goodwin, J., Rozovskaya, O., and Oacyna, J.: Inventory Review 2005 – Emission Data reported to CLRTAP and under the NEC Directive – Initial review for HMs and POPs, EMEP Status report, Norwegian Meteorological Institute, Oslo, 2005. 

Vinken, G. C. M., Boersma, K. F., Maasakkers, J. D., Adon, M., and Martin, R. V.: Worldwide biogenic soil NOx emissions inferred from OMI NO2 observations, Atmos. Chem. Phys., 14, 10363–10381,, 2014. 

Yin, Y., Chevallier, F., Ciais, P., Broquet, G., Fortems-Cheiney, A., Pison, I., and Saunois, M.: Decadal trends in global CO emissions as seen by MOPITT, Atmos. Chem. Phys., 15, 13433–13451,, 2015.  

Yumimoto, K. and Uno, I.: Adjoint inverse modeling of CO emissions over Eastern Asia using four-dimensional variational data assimilation, Atmos. Environ., 40, 6836–6845,, 2006. 

Wang, Y., Broquet, G., Ciais, P., Chevallier, F., Vogel, F., Kadygrov, N., Wu, L., Yin, Y., Wang, R., and Tao, S.: Estimation of observation errors for large-scale atmospheric inversion of CO2 emissions from fossil fuel combustion, Tellus B, 69, 1325723,, 2017. 

Wang, Y., Broquet, G., Ciais, P., Chevallier, F., Vogel, F., Wu, L., Yin, Y., Wang, R., and Tao, S.: Potential of European 14CO2 observation network to estimate the fossil fuel CO2 emissions via atmospheric inversions, Atmos. Chem. Phys., 18, 4229–4250,, 2018. 

WHO: Ambient Air Pollution: a global assessment of exposure and burden of disease, World Health Organization, 2016. 

Zheng, T., French, N. H. F., and Baxter, M.: Development of the WRF-CO2 4D-Var assimilation system v1.0, Geosci. Model Dev., 11, 1725–1752,, 2018. 

Zheng, B., Chevallier, F., Yin, Y., Ciais, P., Fortems-Cheiney, A., Deeter, M. N., Parker, R. J., Wang, Y., Worden, H. M., and Zhao, Y.: Global atmospheric carbon monoxide budget 2000–2017 inferred from multi-species atmospheric inversions, Earth Syst. Sci. Data, 11, 1411–1436,, 2019. 

Short summary
Up-to-date and accurate emission inventories for air pollutants are essential for understanding their role in the formation of tropospheric ozone and particulate matter, for anticipating pollution peaks and for identifying the key drivers that could help mitigate their emissions. Complementarily with bottom-up inventories, the system described here aims at updating and improving the knowledge on the high spatiotemporal variability of emissions of air pollutants.