Monitoring CO emissions of the metropolis Mexico City using TROPOMI CO observations

The Tropospheric Monitoring Instrument (TROPOMI) on the ESA Copernicus Sentinel-5 satellite (S5-P) measures carbon monoxide (CO) total column concentrations as one of its primary targets. In this study, we analyze TROPOMI observations over Mexico City in the period 14 November 2017 to 25 August 2019 by means of collocated CO simulations using the regional Weather Research and Forecasting coupled with Chemistry (WRF-Chem) model. We draw conclusions on the emissions from different urban districts in the region. Our WRF-Chem simulation distinguishes CO emissions from the districts Tula, Pachuca, Tulancingo, Toluca, Cuernavaca, Cuautla, Tlaxcala, Puebla, Mexico City, and Mexico City Arena by 10 separate tracers. For the data interpretation, we apply a source inversion approach determining per district the mean emissions and the temporal variability, the latter regularized to reduce the propagation of the instrument noise and forward-model errors in the inversion. In this way, the TROPOMI observations are used to evaluate the Inventario Nacional de Emisiones de Contaminantes Criterio (INEM) inventory that was adapted to the period 2017–2019 using in situ ground-based observations. For the Tula and Pachuca urban areas in the north of Mexico City, we obtain 0.10±0.004 and 0.09±0.005 Tg yr⁻¹ CO emissions, which exceeds significantly the INEM emissions of <0.008 Tg yr⁻¹ for both areas. On the other hand for Mexico City, TROPOMI estimates emissions of 0.14±0.006 Tg yr⁻¹ CO, which is about half of the INEM emissions of 0.25 Tg yr⁻¹, and for the adjacent district Mexico City Arena the emissions are 0.28±0.01 Tg yr⁻¹ according to TROPOMI observations versus 0.14 Tg yr⁻¹ as stated by the INEM inventory. Interestingly, the total emissions of both districts are similar (0.42±0.016 Tg yr⁻¹ TROPOMI versus 0.39 Tg yr⁻¹ adapted INEM emissions). Moreover, for both areas we found that the TROPOMI emission estimates follow a clear weekly cycle with a minimum during the weekend. This agrees well with ground-based in situ measurements from the Secretaría del Medio Ambiente (SEDEMA) and Fourier transform spectrometer column measurements in Mexico City that are operated by the Network for the Detection of Atmospheric Composition Change Infrared Working Group (NDACC-IRWG). Overall, our study demonstrates an approach to deploying the large number of TROPOMI CO data to draw conclusions on urban emissions on sub-city scales for metropolises like Mexico City. Moreover, for the exploitation of TROPOMI CO observations our analysis indicates the clear need for further improvements of regional models like WRF-Chem, in particular with respect to the prediction of the local wind fields.

1 Introduction

Carbon monoxide (CO) is an atmospheric trace gas emitted by incomplete combustion to the atmosphere (e.g., biomass burning, industrial activity, and traffic). Its background concentration is relatively low with an atmospheric residence time varying from days to month (Holloway et al., 2000) depending on the atmospheric concentration of the hydroxyl radical (Spivakovsky et al., 2000). These characteristics established CO as a tracer for air pollution and transport processes in the atmosphere (e.g., Gloudemans et al., 2009; Pommier et al., 2013; Schneising et al., 2020).

The Tropospheric Monitoring Instrument (TROPOMI) launched in 2017 as the single payload of ESA's Copernicus Sentinel-5 Precursor mission aims to monitor CO as one of its primary targets. The operational CO column product is inferred from TROPOMI's shortwave infrared measurements with daily global coverage and a high spatial resolution of 7×7 km² (Veefkind et al., 2012). Early in the mission, the TROPOMI CO dataset was validated with ground-based measurements of the Total Carbon Column Observing Network (TCCON) (Borsdorff et al., 2018a) and inter-compared with simulated CO fields of the European Centre for Medium-Range Weather Forecasts (ECMWF) Integrated Forecasting System (Borsdorff et al., 2018b). On 11 July 2018, it was concluded that the TROPOMI CO data quality is fully compliant with the mission requirements of 15 % precision and 10 % accuracy and so it was released for public usage (https://scihub.copernicus.eu, last access: 16 December 2020).

Borsdorff et al. (2018a, 2019) illustrated the capability of TROPOMI to detect CO emissions from pollution hot spots of medium-sized to large cites (e.g., Yerevan, Tabriz, Urmia, and Tehran), industrial areas (e.g., Po Valley in Italy), and even from pollution along arterial roads in Armenia. To monitor the emissions of metropolises, data interpretation of multi-annual datasets is required. The different inversion techniques discussed by Varon et al. (2018) for plume inversions, i.e., the source pixel method, the mass balance method, and the inversion of a Gaussian plume model, are appropriate for interpreting emissions of point sources but are less suitable for flux inversion of spatially extended sources. Therefore, in this study we estimate CO emissions by inverting simulations of the regional atmospheric modeling Weather Research and Forecasting coupled with Chemistry (WRF-Chem) as an atmospheric tracer transport model, which allows us to simulate the CO column at the same spatial resolution as TROPOMI. Possible error sources of this type of flux inversion are the limited validity of the simulated wind fields, prior assumptions on the spatial distribution of emissions, and the simulated atmospheric dispersion (Borsdorff et al., 2019).

Mexico City is a prime example of a CO pollution hot spot that is clearly detectable by TROPOMI. It is a fast-growing megacity located at an altitude of 2240 m on the Central Plateau which is surrounded by mountains. The urban area is divided into 10 different urban districts (Tula, Pachuca, Tulancingo, Toluca, Cuernavaca, Cuautla, Tlaxcala, Puebla, Mexico City (CdMx), and Mexico City Arena (ACdMx)), and the metropolis has a long history of atmospheric pollution measurements. More than 29 in situ CO measurement stations are distributed over the city, operated by the Secretaría del Medio Ambiente (SEDEMA, Mexican Ministry of the Environment). About every 2 years, the ministry reports on the CO emissions of Mexico City. Based on a bottom-up approach using the in situ measurements, it is concluded that a major part of Mexico City's CO emissions is caused by light-duty motor vehicles with a significant decline in recent years. For the zona metropolitana del valle de México (ZMVM, Greater Mexico City), SEDEMA reported a reduction in CO emissions from 2.04 to 0.7 and 0.28 Tg yr⁻¹ in the years 2000, 2014, and 2016 (SMA-GDF, 2018).

These in situ measurements are complemented by ground-based Fourier transform infrared (FTIR) observations of the NDACC (Network for the Detection of Atmospheric Composition Change) IRWG (Infrared Working Group) network, which among other products regularly provides CO total column concentrations. Using NDACC and IASI satellite observations of CO, Stremme et al. (2013) estimated the overall annual CO emissions of Mexico City to be about 2.15 Tg yr⁻¹ for the year 2008. Building on this, TROPOMI CO observations add new possibilities for air quality monitoring due to its regional coverage and daily overpass combined with the high precision of its data.

In this study, we analyze about 2 years of TROPOMI CO measurements using collocated WRF-Chem CO simulations for Mexico to get more insight into the emissions of Mexico City. To this end, in Sect. 1 we introduce the TROPOMI CO dataset and the simulation of the WRF-Chem model, and Sect. 2 describes our methodology to fit the WRF-Chem model to the TROPOMI data for emission estimates. Section 3 discusses our findings, and Sect. 4 gives the summary and conclusion.

2 TROPOMI CO dataset

To investigate CO emissions of the Mexico City metropolis, we select the TROPOMI dataset of CO total column observations between 14 November 2017 and 25 August 2019 over Mexico. The data are processed with the shortwave infrared CO retrieval (SICOR) algorithm that was developed for Copernicus operational data processing (Landgraf et al., 2016a). Algorithm settings like the spectral windows, a priori profiles, and other auxiliary data are reported by Landgraf et al. (2016b). The SICOR algorithm accounts for atmospheric scattering by retrieving effective cloud parameters (altitude, optical thickness) together with the total column concentrations of CO and of the interfering gases H₂O, HDO, and CH₄. The radiative transfer simulation uses the HITRAN2016 database for all species as described by Borsdorff et al. (2019), and the inversion deploys the profile-scaling approach that scales a reference profile to fit the spectral measurement (Borsdorff et al., 2014). Here, the a priori profile is taken from spatio-temporally resolved atmospheric transport simulations of the TM5 model (Krol et al., 2005). The TROPOMI CO data product includes the total column averaging kernel a_col that relates the true vertical CO profile ρ_true to the retrieved total column concentration c_ret following the equation

$$\begin{array}{}\text{(1)}& c_{\mathrm{ret}}={\mathit{a}}_{\mathrm{col}}{\mathit{\rho }}_{\mathrm{true}}+\mathit{ϵ},\end{array}$$

with the noise contribution ϵ. This study limits the analysis to scenes under clear-sky and low-cloud atmospheric conditions. This corresponds to a quality assurance value q>0.5 which is also provided by the S5P data product. Finally, the individual TROPOMI CO orbits show an artificial striping in flight direction, probably due to a deficient instrument calibration. To reduce this feature, we apply an a posteriori data correction as discussed by Borsdorff et al. (2019) based on frequency filtering in the Fourier space. Finally, on 5 August 2019, the spatial sampling of the data product with satellite nadir geometry was improved from 7×7 to 7×5.6 km² due to a shorter readout time of the detectors. This event is covered by our dataset.

Figure 1 Urban districts surrounding Mexico City. For each of the color-coded domains a separate WRF-Chem tracer run was performed based on the emissions within the polygons. The elevation map in the background is under copyright © Esri, Airbus DS, USGS, NGA, NASA, CGIAR, N Robinson, NCEAS, NLS, OS, NMA, Geodatastyrelsen, Rijkswaterstaat, GSA, Geoland, FEMA, Intermap, and the GIS user community.

3 Methodology

3.1 The WRF-Chem model

We simulate the TROPOMI CO column concentrations by deploying the WRF-Chem model version 3.9.1.1. The simulation covers the time period of TROPOMI measurements on the regional domain shown in Fig. 1. We ignore photo-chemical oxidation and secondary production of CO in the atmosphere (chem_opt option 106 (RADM2-KPP), as a tracer with gaschem off), which is justified by the long lifetime of CO compared with the size of the model domain as discussed by Dekker et al. (2017). Especially for the region of Mexico City, the contribution of atmospheric chemistry to the total CO concentration is less than 3 % as presented by Mejia (2020). Hence, WRF-Chem simulates the transport of CO surface emissions as traces as done by e.g., Borsdorff et al. (2019) and Dekker et al. (2017, 2019). The spatial resolution of the simulation is chosen to be comparable with the TROPOMI CO product sampling. Each grid cell of the considered simulation domain (270×270 km²) is 3×3 km². The WRF-Chem simulation employs the emission inventory Inventario Nacional de Emisiones de Contaminantes Criterio (INEM) for the year 2013 but scaled by a factor of 0.48 to make it applicable to the years 2017 to 2019. Here the scaling factor is based on recent SEDEMA surface measurements (García-Reynoso et al., 2018). The inventory includes diurnal, week-to-week, and monthly variations in the CO emissions, where weekly and daily temporal profiles are derived from traffic counts in Mexico. The inventory is described in more detail by García-Reynoso et al. (2018). Finally, the model run is constrained by National Centers for Environmental Prediction (NCEP) North American Mesoscale (NAM) 12 km analysis wind fields (NCEP, 2015) and yields vertical CO concentration profiles for every latitude and longitude grid cell and every model time step and tracer run.

The WRF-Chem simulation uses 10 independent tracers to estimate the CO emissions of the areas Tula, Pachuca, Tulancingo, Toluca, Cuernavaca, Cuautla, Tlaxcala, Puebla, the metropolitan area of Mexico City CdMx, and the adjoint urban area ACdMx. The total simulated CO field is given by the sum of the simulated CO fields of the tracer. Since no atmospheric chemistry is accounted for, each CO tracer field is linear in terms of the corresponding emissions per district:

$$\begin{array}{}\text{(2)}& {\mathrm{F}}_{\mathrm{WRF}}({\mathit{\alpha }}_{\mathrm{1}},\mathrm{\cdots },{\mathit{\alpha }}_{\mathrm{10}})=\sum _{i=\mathrm{1}}^{\mathrm{10}}{\mathbf{k}}_i{\mathit{\alpha }}_i,\end{array}$$

where α_i is the corresponding scaling factor and k_i represents the CO tracer field for the reference emissions (adapted INEM data) for α_i=1.

Before using our model to simulate TROPOMI data, we interpolate the model fields to the geolocation and time of the TROPOMI observations. Subsequently, we integrate the model CO profiles to total column densities by applying the total column averaging kernel of the TROPOMI CO retrieval following Eq. (1). We summarize this numerical step in the observation operator 𝒪, which transforms the forward model into

$$\begin{array}{}\text{(3)}& {\mathrm{F}}_{\mathrm{sat}}({\mathit{\alpha }}_{\mathrm{1}},\mathrm{\cdots },{\mathit{\alpha }}_{\mathrm{10}})=\sum _{i=\mathrm{1}}^{\mathrm{10}}\mathcal{O}\left({\mathbf{k}}_i\right){\mathit{\alpha }}_i.\end{array}$$

Hence, the operator 𝒪 accounts for the TROPOMI-specific vertical sensitivity, which can change from measurement to measurement and so ensures that the comparison between TROPOMI and WRF-Chem is free of the null space or smoothing error (Rodgers, 2000; Borsdorff et al., 2014). Here, the scaling factors α_i are not affected by the operation. In a next step, we transform Eq. (3) to

$$\begin{array}{}\text{(4)}& {\mathrm{F}}_{\mathrm{sat}}(E_{\mathrm{1}},\mathrm{\cdots },E_{\mathrm{10}})=\sum _{i=\mathrm{1}}^{\mathrm{10}}\mathcal{O}\left({\stackrel{\mathrm{̃}}{\mathbf{k}}}_i\right)E_i,\end{array}$$

where ${\stackrel{\mathrm{̃}}{\mathbf{k}}}_i=\frac{{\mathbf{k}}_i}{E_{i,\mathrm{INEM}}}$ and $E_i={\mathit{\alpha }}_i{E}_{i,\mathrm{INEM}}$ with the corresponding emissions E_i,INEM of the INEM inventory interpolated to the TROPOMI overpass time.

To improve the capability of our forward model to fit TROPOMI observations, we introduce a spatially constant CO background field k_bg and an altitude dependence term ${\mathbf{k}}_{\mathrm{elv}}=z-z_{\mathrm{ref}}$ with corresponding scaling factors α_bg and α_elv. Here, z is the respective elevation of the TROPOMI CO ground pixels and z_ref=2240 m is an arbitrary reference altitude set to the elevation of Mexico City:

$$\begin{array}{}\text{(5)}& \begin{array}{rl}& {\mathrm{F}}_{\mathrm{sat}}(E_{\mathrm{1}},\mathrm{\cdots },E_{\mathrm{10}},{\mathit{\alpha }}_{\mathrm{bg}},{\mathit{\alpha }}_{\mathrm{elv}})\\ & =\sum _{i=\mathrm{1}}^{\mathrm{10}}\mathcal{O}\left({\stackrel{\mathrm{̃}}{\mathbf{k}}}_i\right)E_i+{\mathbf{k}}_{\mathrm{bg}}{\mathit{\alpha }}_{\mathrm{bg}}+{\mathbf{k}}_{\mathrm{elv}}{\mathit{\alpha }}_{\mathrm{elv}}.\end{array}\end{array}$$

These two effective model components account for the CO contribution over the Mexico City area originating from outside the model domain such as from fires, power plants, biogenic production, and other cities as well as long-range transport (Borsdorff et al., 2019) and an altitude-dependent linear vertical gradient of the CO columns. Both do not interfere with any localized emission sources. They mitigate shortcomings of the WRF-Chem simulations ignoring CO boundary conditions in the model domain.

Finally, for the interpretation of our CO forward simulations, we make an important assumption. Although the WRF-Chem simulations account for the temporal accumulation of the localized CO emissions over days and weeks, we allocate an emission estimate of the corresponding overpass time to each TROPOMI overpass. Here, we assume that a TROPOMI CO image is dominated by the emissions of the urban districts for the particular observation day, where the temporal accumulation of CO from previous days is partly described by the WRF-Chem simulation due to the corresponding scaling of the inventory and partly mitigated by fitting the nuisance parameters α_bg and α_elv.

3.2 Inversion methodology

Interpreting a series of n TROPOMI CO images,

$$\begin{array}{}\text{(6)}& \mathit{y}=({\mathit{y}}_{\mathrm{1}},\mathrm{\cdots },{\mathit{y}}_n),\end{array}$$

at overpass times $t_{\mathrm{0}},\mathrm{\cdots },t_n$ means estimating the corresponding emissions given by the state vector,

$$\begin{array}{}\text{(7)}& \mathit{x}=({\mathit{x}}_{\mathrm{1}},\mathrm{\cdots },{\mathit{x}}_n),\end{array}$$

where each element comprises

$$\begin{array}{}\text{(8)}& {\mathit{x}}_i=(E_{\mathrm{1},i},\mathrm{\cdots },E_{\mathrm{10},i},{\mathit{\alpha }}_{\mathrm{bg},i},{\mathit{\alpha }}_{\mathrm{elv},i})\end{array}$$

at the corresponding time t_i. Our linear forward model in Eq. (5) describes the measurement vector y by

$$\begin{array}{}\text{(9)}& \left(\begin{array}{l}{\mathit{y}}_{\mathrm{1}}\\ {\mathit{y}}_{\mathrm{2}}\\ \mathrm{⋮}\\ {\mathit{y}}_n\end{array}\right)=\left(\begin{array}{cccc}{\mathbf{K}}_{\mathrm{1}}& \mathrm{0}& \mathrm{\cdots }& \mathrm{0}\\ \mathrm{0}& {\mathbf{K}}_{\mathrm{2}}& \mathrm{\cdots }& \mathrm{0}\\ \mathrm{⋮}& \mathrm{⋮}& \mathrm{\ddots }& \mathrm{⋮}\\ \mathrm{0}& \mathrm{0}& \mathrm{\cdots }& {\mathbf{K}}_n\end{array}\right)\left(\begin{array}{l}{\mathit{x}}_{\mathrm{1}}\\ {\mathit{x}}_{\mathrm{2}}\\ \mathrm{⋮}\\ {\mathit{x}}_n\end{array}\right),\end{array}$$

with the forward-model Jacobian matrix ${\mathbf{K}}_i=\left(\mathcal{O}\right({\stackrel{\mathrm{̃}}{\mathbf{k}}}_{\mathrm{1},i}),\mathrm{\cdots }\mathcal{O}({\stackrel{\mathrm{̃}}{\mathbf{k}}}_{\mathrm{10},i}),{\mathbf{k}}_{\mathrm{bg},i},{\mathbf{k}}_{\mathrm{elv},i})$, in short y=Kx. Equation (9) can be inverted by

$$\begin{array}{}\text{(10)}& {\mathit{x}}_{\mathrm{est}}=\underset{\mathit{x}}{min}\left\{\left|\right|\mathit{y}-\mathbf{K}\mathit{x}|{|}_{{\mathbf{S}}_{\mathrm{e}}}^{\mathrm{2}}\right\},\end{array}$$

which is equivalent to the solution $({\mathit{x}}_{\mathrm{est},\mathrm{1}},\mathrm{\cdots },{\mathit{x}}_{\mathrm{est},n})$ of the individual problems y_i=K_ix_i due to the block diagonal form of Eq. (9). Here, the norm of an arbitrary vector p is defined by $\left|\right|\mathit{p}|{|}_{{\mathbf{S}}_{\mathrm{e}}}^{\mathrm{2}}={\mathit{p}}^T{\mathbf{S}}_{\mathrm{e}}^{-\mathrm{1}}\mathit{p}$ and S_e is the measurement error covariance matrix with the variance of the TROPOMI retrieval error on the diagonal.

Due to measurement noise and forward-model errors, the least-squares inversion of Eq. (10) results in unfavorable error propagation and so requires regularization. Because our problem is linear in the state vector x, regularization can be performed as part of the fitting approach or a posteriori to the least-squares solution, without loss of generality. To regularize the noise propagation, we first derive the temporal mean,

$$\begin{array}{}\text{(11)}& {\overline{\mathit{x}}}_{\mathrm{est}}={\displaystyle \frac{\mathrm{1}}{n}}\sum _{i=\mathrm{1}}^n{\mathit{x}}_{\mathrm{est},i},\end{array}$$

from the unregularized solution in Eq. (10). This modifies our cost function to

$$\begin{array}{}\text{(12)}& {\mathit{x}}_{\mathrm{est}}=\underset{\mathit{x}}{min}\left\{\left|\right|\mathit{y}-\mathbf{K}\mathit{x}|{|}_{{\mathbf{S}}_{\mathrm{e}}}^{\mathrm{2}}\right\}\text{with}\overline{\mathit{x}}={\displaystyle \frac{\mathrm{1}}{n}}\sum _{i=\mathrm{1}}^n{\mathit{x}}_i.\end{array}$$

In this way, we divided the solution of the original inversion problem (10) into two steps: first we determine the mean emissions from the individual least-squares solutions x_est,i, which yields the constrained least-squares problem in Eq. (12) to describe the temporal variability. The side constraint guarantees that measurement information is not used twice. Finally, we add an additional Tikhonov side constraint to Eq. (12) to regularize the error propagation,

$$\begin{array}{}\text{(13)}& {\mathit{x}}_{\mathrm{est}}=\underset{\mathit{x}}{min}\left\{\left|\right|\mathit{y}-\mathbf{K}\mathit{x}|{|}_{{\mathbf{S}}_{\mathrm{e}}}^{\mathrm{2}}+||\mathit{x}-\overline{\mathit{x}}|{|}_{\mathbf{\Gamma }}^{\mathrm{2}}\right\},\end{array}$$

with

$$\begin{array}{}\text{(14)}& \overline{\mathit{x}}={\displaystyle \frac{\mathrm{1}}{n}}\sum _{i=\mathrm{1}}^n{\mathit{x}}_i.\end{array}$$

Here, Γ is an appropriate regularization matrix. For a block diagonal form of Γ analogous to the Jacobian matrix in Eq. (9), namely

$$\begin{array}{}\text{(15)}& \mathbf{\Gamma }=\left(\begin{array}{cccc}{\mathbf{\Gamma }}_{\mathrm{1}}& \mathrm{0}& \mathrm{\cdots }& \mathrm{0}\\ \mathrm{0}& {\mathbf{\Gamma }}_{\mathrm{2}}& \mathrm{\cdots }& \mathrm{0}\\ \mathrm{⋮}& \mathrm{⋮}& \mathrm{\ddots }& \mathrm{⋮}\\ \mathrm{0}& \mathrm{0}& \mathrm{\cdots }& {\mathbf{\Gamma }}_n\end{array}\right),\end{array}$$

the minimization problem (13) decomposes into n problems,

$$\begin{array}{}\text{(16)}& {\mathit{x}}_{\mathrm{est},\mathrm{i}}=\underset{{\mathit{x}}_{\mathit{i}}}{min}\left\{\left|\right|{\mathit{y}}_{\mathit{i}}-{\mathbf{K}}_{\mathbf{i}}{\mathit{x}}_{\mathit{i}}|{|}_{{\mathbf{S}}_{\mathrm{e}}}^{\mathrm{2}}+||{\mathit{x}}_{\mathit{i}}-\overline{\mathit{x}}|{|}_{\mathbf{\Gamma }}^{\mathrm{2}}\right\},\end{array}$$

which are only coupled by the external side constraint (14). A closer look at our inversion problem shows that the two constraints have similar effects. The Tikhonov constraint $\left|\right|\mathit{x}-\overline{\mathit{x}}|{|}_{{\mathbf{S}}_{\mathrm{e}}}$ minimizes the variation in the state vector around its mean depending on the regularization parameter λ, whereas the external constraint requires strict conservation of the mean.

Therefore, in practice, we solve the inversion (16) and evaluate the external constraint on the mean afterwards to confirm proper use of the measurement information. Its solution is given by

$$\begin{array}{}\text{(17)}& {\mathit{x}}_{\mathrm{est},\mathrm{i}}={\mathbf{G}}_{\mathrm{i}}({\mathit{y}}_i-{\mathbf{K}}_i{\overline{\mathit{x}}}_{\mathit{i}})+{\overline{\mathit{x}}}_i\end{array}$$

with the gain matrix

$$\begin{array}{}\text{(18)}& {\mathbf{G}}_i=({\mathbf{K}}_i^T{\mathbf{S}}_{\mathrm{e},\mathrm{i}}^{-\mathrm{1}}{\mathbf{K}}_i+{\mathbf{\Gamma }}_i{)}^{-\mathrm{1}}{\mathbf{K}}_i^T{\mathbf{S}}_{\mathrm{e},\mathrm{i}}^{-\mathrm{1}}.\end{array}$$

The inversion's averaging kernel relates the “true” state vector x_true,i to x_est,i; namely

$$\begin{array}{}\text{(19)}& {\mathit{x}}_{\mathrm{est},i}={\mathbf{A}}_i\left({\mathit{x}}_{\mathrm{true},i}-{\mathit{x}}_i\right)+\overline{\mathit{x}}\end{array}$$

with

$$\begin{array}{}\text{(20)}& {\mathbf{A}}_i={\mathbf{G}}_i{\mathbf{K}}_i.\end{array}$$

A_i represents the derivative ${\mathbf{A}}_{\mathrm{i},\mathrm{kl}}=\frac{\partial {\mathit{x}}_{\mathrm{est},j}}{\partial {\mathit{x}}_{\mathrm{true},l}}$, where its diagonal elements describe the retrieval sensitivity of a state vector element to its true value. The degree of freedom for signal

$$\begin{array}{}\text{(21)}& {\mathrm{DFS}}_i=\mathrm{trace}\left({\mathbf{A}}_i\right),\end{array}$$

indicates the total number of independent pieces of information.

To evaluate the fit quality for each overpass, we consider the fit residuals ${\mathit{\delta }}_i={\mathit{y}}_i-{\mathbf{K}}_i{\mathit{x}}_{\mathrm{est},i}$. Additionally, we evaluate the goodness of the fit described by the reduced chi-squared value,

$$\begin{array}{}\text{(22)}& {\mathit{\chi }}_i^{\mathrm{2}}={\displaystyle \frac{\mathrm{1}}{{\mathit{\nu }}_i}}\sum _{l=\mathrm{1}}^L({\mathit{\delta }}_i,l/{\mathit{y}}_{\mathrm{err},ik}{)}^{\mathrm{2}}.\end{array}$$

Here L is the number of observations of a single overpass, y_err,i is the retrieval error, and ${\mathit{\nu }}_i=L-{\mathrm{DFS}}_i$.

3.3 Estimate of the mean emissions

The first step of the inversion described in the previous section means determining the a priori emissions from a set of TROPOMI data with highest information content using a unregularized least-squares fit, Γ=0. Here, the individual emission estimates may be noisy due to enhanced error propagation in the inversion; however, averaging all inversions reduces noise contribution and so gives a reliable estimate of mean emissions for the different districts. The validity of this approach highly depends on the selected dataset of TROPOMI overpasses. On the one hand, the ensemble should be large enough to estimate mean emissions for the considered time period, but on the other hand it should be strictly filtered for cases where the forward model is in good agreement with the measurement such that a stable inversion of all the emissions is possible. The information content of a single overpass varies and depends on several aspects: (1) the number of useful measurements and their cloud coverage changes between different TROPOMI overpasses. Here, clouds shield the lower troposphere, where atmospheric measurements are particularly sensitive to the surface emissions E_i. (2) The pixel size at the swath edge is about 32 km and so about 5 times larger than at the sub-satellite point. This reduces not only the number of pixels covering a certain area but also the sensitivity of the individual TROPOMI observations. (3) The quality of the forward model depends on the meteorological situation, where we consider model simulations for low wind speeds more reliable. These considerations lead to the criteria of the data filtering to determine the mean emissions for each district. We only select overpasses which meet both filter criteria:

• TROPOMI observations cover 70 % of the data domain.

• For all observations the across-track pixel size is <15 km.

The filter criteria reduce the original set of 551 overpasses to 199, which we consider to be sufficient to estimate the overall average emission rate per district, yielding $\overline{\mathit{x}}$. For this we use the median instead of the mean because of its robustness against outliers. With the same reasoning we define the percentile difference,

$$\begin{array}{}\text{(23)}& \mathit{\delta }{P}_j=\mathrm{|}{\displaystyle \frac{P_j\left(\mathrm{84.1}\right)-P_j\left(\mathrm{15.9}\right)}{\mathrm{2}}}\mathrm{|},\end{array}$$

to describe scattering in the data, which corresponds to the standard deviation of normally distributed parameters. Finally, we calculate the error in the mean using the percentile difference.

3.4 Final data product

Subsequently, the final data reduction step is performed solving the inversion problem (13). For all overpasses, we choose Γ_i to be a diagonal matrix with

$$\begin{array}{}\text{(24)}& \mathrm{diag}{\mathbf{\Gamma }}_i=[{\mathit{\gamma }}_{\mathrm{1}},{\mathit{\gamma }}_{\mathrm{2}},\mathrm{\cdots },{\mathit{\gamma }}_{\mathrm{10}},\mathrm{0},\mathrm{0}],\end{array}$$

where the zeros ensure that the elements of the state vector α_bg and α_elv are not regularized. Obviously, the regularization parameter γ_k must be well-chosen to optimize the balance between minimum error propagation on the fit parameter and maximum information content inferred from the measurement. If γ_k is chosen to be too small, the propagation of the TROPOMI measurement noise as well as retrieval biases and forward-model errors dominates the inversion. If γ_k is chosen to be too large, the estimated state vector reproduces the a priori estimate without appropriate use of the information content of the measurement. For our application, we fix the regularization parameter γ_k for $k=\mathrm{1},\mathrm{\cdots },\mathrm{10}$ to constant values such that the scatter of the retrieved emissions stays within predefined boundaries.

Considering the temporal variation in the INEM emissions to be about 40 %, we adjusted the regularization parameter ${\mathit{\gamma }}_{\mathrm{1}},\mathrm{\cdots },{\mathit{\gamma }}_{\mathrm{10}}$ such that the retrieved emissions vary within 60 % of their average. The value 60 % is empirically chosen to balance information content against noise propagation. It puts a moderate constraint on the inversion, ensuring on the one hand a stable inversion and on the other hand a realistic variation in the retrieved emissions around the a priori.

One great advantage of the final retrieval product is that it includes the averaging kernel A_i. This can be used to filter the data with respect to the information provided by the TROPOMI measurements. For the emissions of each tracer, we filter on the individual emissions E_l, considering averaging-kernel values ${\mathbf{A}}_{l,l}>\mathrm{0.3}$. This form of data mining optimizes the data use, keeping in mind that TROPOMI overpasses may be appropriate for determining one specific source but not all sources simultaneously. The concept of information-content-based filtering turned out to be very useful and enhances the data exploitation compared to the unregularized least-squares fitting used to determine the mean emission values.

Figure 2 Background CO concentration for the domain shown in Fig. 1 estimated by fitting the WRF-Chem simulation to the TROPOMI data. (a) Background CO for individual collocations from 9 November 2017 to 25 August 2019. (b) Monthly mean background CO based on the individual collocations. The error bars are the standard error in the mean, and the light-blue line is the time-invariant a priori used in the fit.

4 Results

Figure 2 shows the CO background that was fitted as an auxiliary parameter during the inversion described in Sect. 3.2. The concentration and its annual cycle are shown. Here, the biomass burning season between February and June causes the corresponding CO enhancement, whereas lower CO concentrations are observed during the rainy season between June and November. The extremely high CO column values on 15 May 2019 are due to the transport of CO-enriched air from wildfires in the southwest of Mexico into the model domain. Figure 3 shows the CO concentration in the state of Mexico under normal conditions and after the fires, which caused a serious health hazard in Mexico City. These types of fires outside the model domain create an inhomogeneous background CO field over Mexico City, which cannot be described by our forward model. Only fitting a constant background is not sufficient in these extreme cases, and so during the fire season many data cannot be used (we excluded the months May and June 2019).

Figure 3 TROPOMI CO data over Mexico City averaged on a 0.1 by 0.1^∘ lat–long grid. (a) Averaged from 12 to 18 April 2019 showing undisturbed background CO levels. (b) Averaged from 12 to 18 May 2019 showing high CO concentrations in Mexico City caused by fires in the southeast. The street map in the background is under copyright © 2009 ESRI, AND, TANA, ESRI Japan, UNEP-WCMC.

Figure 4 Example cases for fitting the WRF-Chem simulation to the TROPOMI data deploying the “final-fit” approach for (a) 20 September, (b) 7 November, (c) 19 November 2018, and (d) 17 August 2019. TROPOMI CO retrievals (left column), WRF-Chem simulation fitted to the TROPOMI data (middle column), and the residual (right column, TROPOMI – WRF-Chem).

Figure 4 shows three examples of TROPOMI overpasses, which include a pixel resolution of 7×7 km² (panels a, b, and c) and the enhanced spatial resolution of 7×5.6 km² (panel d), where the latter represents the TROPOMI instrument baseline since 6 August 2019. Focusing on the dry season, the TROPOMI instrument can detect distinct CO enhancements over the different emission areas in central Mexico with the retrievals from single-orbit overpasses (see left column of Fig. 4). After fitting our forward model to the TROPOMI measurements, simulated data and observations are brought into good agreement as illustrated in Fig. 4. Particularly for low-wind-speed conditions in Fig. 4a, TROPOMI and WRF-Chem show distinct CO enhancements over the different emission areas of Mexico. Furthermore, the transport of CO-enhanced air from Mexico City towards the south following the mountain orography and the accumulation of CO in the south is seen by TROPOMI in agreement with the WRF-Chem simulation (Fig. 4c). This clearly shows that regional models like WRF-Chem have a great potential to reproduce the large-scale patterns seen by the TROPOMI instrument. However, we also found clear localized residuals in the difference δ_j between observations and the forward model (right column of Fig. 4). For atmospheric conditions under high wind speeds the WRF-Chem simulations can deviate more from the TROPOMI measurements as shown in Fig. 4c. Here, the plume of CO-enriched air extending from Mexico City towards the north is simulated as very narrow compared to the more dispersed plume seen by TROPOMI. This points to a possible underestimation of the atmospheric dispersion in the WRF-Chem simulation. A very prominent residual between TROPOMI and WRF-Chem not only is shown in Fig. 4d but also is present in Fig. 4a and b. Here TROPOMI measures a strong CO enhancement in the north of Mexico City that is not reproduced by the WRF-Chem model. This points at a deficient spatial distribution of INEM emissions.

Figure 5 Statistics of CO emissions averaged from 9 November 2017 to 25 August 2019 for the tracer domains shown in Fig. 1. (a) Median of the a priori emissions (adapted INEM inventory) used for the WRF-Chem simulation (grey) and retrieved from the TROPOMI data (unregularized fit in blue, regularized fit in green). The error bars indicate the standard error in the mean calculated from the delta percentiles (b) used as a robust estimation of the standard deviation; (c) the median of the goodness of the fit (χ²), and (d) the number of collocations. The number of collocations and the χ² values of the a priori simulation and unregularized fit are the same for all tracer domains (blue and grey line), but in the final regularized fit they change due to the information content filtering. Here, a collocation corresponds to a specific day because TROPOMI overpasses the region only once.

For each tracer domain Fig. 5a shows the mean emissions of the adapted INEM inventory, the unregularized least-squares fit, and the final data product. The mean emissions agree very well between the last two approaches, indicating that the final inversion in Eq. (17) satisfies the constraint of the predefined mean value. This supports our assumption that the external constraint does not have to be accounted for in the chosen Tikhonov constraint of the inversion. The scatter of the least-squares product is high and in most cases exceeds 100 % (see Fig. 5b), which is expected for the non-constrained inversion. Moreover, we find significant differences between the emissions of the a priori and the final data product. The retrieved emissions for the urban districts Tula (0.10±0.004 Tg yr⁻¹) and Pachuca (0.09±0.005 Tg yr⁻¹) in the north of Mexico City seem to be underestimated by the emission inventory (both were less then 0.008 Tg yr⁻¹). It is not yet clear what sources are missing in the inventory; this needs to be addressed in future studies. However, we identified an oil refinery and a power plant near to Tula and cement and lime kilns near to Pachuca that could have contributed to the CO emissions. Furthermore, we found that the emissions of the central part of Mexico City (CdMx) are assumed to be too high in the adapted INEM inventory (0.25 Tg yr⁻¹), where TROPOMI measurements indicate lower values for CdMx (0.14±0.006 Tg yr⁻¹). This is accompanied by higher values for the adjacent district ACdMx (0.28±0.01 Tg yr⁻¹). The sum of both emissions (0.42±0.016 Tg yr⁻¹) is similar to the a priori emissions (0.39 Tg yr⁻¹). This may mean that the total emissions of the domain including CdMx and ACdMx is well represented in the emission inventory but only the spatial distribution of the source intensity needs refinement.

Figure 6 (a) The goodness of the fit (χ²) of the TROPOMI CO measurements and the simulation of the WRF-CHEM model for single-orbit overpasses. We distinguish the cases, only fitting background parameters (a priori) and additionally fitting the 10 tracer fields (unregularized and regularized). (b) Improvement of the goodness of the fit (χ²) when fitting the 10 tracer fields (unregularized and regularized) relative to the (a priori) case. For the regularized fit we only show the χ² of the urban district ACdMx because it provides good data coverage when filtering for the degree of freedom.

The χ² values in Fig. 5 clearly show that the agreement between TROPOMI and WRF-Chem can be improved by fitting the emissions of the different city districts (blue line) instead of using the INEM inventory (grey line). The regularization approach increases the χ² values (green bars) because the inversion cannot compensate so well for differences between TROPOMI and WRF-Chem by choosing unrealistic emissions. However, the χ² values of the final fit are still lower than the ones for the a priori INEM emissions (grey line). Overall, the χ² values exceeds 1 which indicates that the difference between TROPOMI and WRF-Chem is dominated by systematic errors in the WRF-Chem simulation. Figure 6 shows the χ² values for the difference between the TROPOMI CO observations of single overpasses and the WRF-Chem simulations over the considered time range of the study. The χ² values follow a seasonal pattern with enhanced values during the biomass burning season between February and June and low values during the rainy season between June and November. As mentioned before, in the vicinity of other pollution sources (e.g., wildfires) the background variability in CO becomes more complex and can interfere with the retrieved local emissions of Mexico City. Hence, a better model of the CO background concentration and its variability is needed to cope with this effect. However, Fig. 6 also shows that fitting the emissions of the different city districts significantly improves the χ² between TROPOMI and WRF-Chem over the whole time range and the improvement can be even more than 50 % for single overpasses compared to cases of simulations only using the emission values provided by the INEM inventory.

Figure 7 Averaging-kernel matrices showing the sensitivity and cross-sensitivities for the scaling of the different tracer fields. The same cases as in Fig. 4 are shown for the dates (a) 20 September, (b) 7 November, (c) 19 November 2018, and (d) 17 August 2019 but deploying the regularized retrieval.

For a correct interpretation of the retrieved emissions the averaging kernel, as shown in Fig. 7 for four example cases, offers several advantages. The figure shows that generally the averaging kernels have high values on the diagonal, indicating high sensitivity to the quantity to be retrieved. It indicates that TROPOMI measurements can be used to distinguish emissions of the different urban districts of Mexico, with the exception of the emissions of the Tulancingo district. Due to the small mean emissions, the averaging kernel indicates a low sensitivity of the data product. Furthermore, the averaging kernel shows cross-correlations between the different elements of the state vector due to the regularization. Although these interdependencies exist, e.g., between the emissions of CdMx and ACdMx as shown in panel (d) of Fig. 7, these are still small compared the diagonal. The averaging-kernel information is very useful to filter the emission product with respect to the information provided by the TROPOMI measurements. Using the sensitivity of individual sources, this results in a different number of coincidences for the different districts (panel c of Fig. 5). This form of data mining optimizes the data use, keeping in mind that TROPOMI overpasses may be appropriate for determining one specific source but not all sources simultaneously. In this manner, error propagation in the inversion can be minimized.

Figure 8 Retrieved CO emissions from the TROPOMI data for the tracers CdMx (left column) and ACdMx (right column). (a) Time series of individually retrieved CO emissions. The error bars indicate the error in the fit, and the black line is the time-invariant a priori used in the fit. (b) Degree of freedom of the scaling factor for the tracer field. Only data with degrees of freedom (dofs) >0.3 are accounted for. (c) Weekly cycle of the CO emissions. Median values are shown, and the error bars are the standard error in the mean deploying the delta percentile as a robust estimation of the standard deviation.

Due to the little scatter and the higher number of data of the final data product for the suburbs CdMx and ACdMx, we can draw conclusions about the time-dependent variability in emissions. Figure 8a shows the time series of the emissions for CdMx and ACdMx, which vary around the a priori value. This temporal variation is determined from the measurements as all a priori information is assumed to be time invariant. The scatter of the data is still high and even includes negative values. Even though negative emissions are not physical, we need to keep them in our analyses because filtering negative noise can induce a positive bias in the mean. Panel (b) of the figure shows relatively high values of the diagonal elements of the averaging kernel for the emissions of the two urban districts. Finally, panel (c) of the figure indicates a clear weekly CO cycle in the data with low values during weekends. During the week, the CO emissions of the two districts do not differ significantly due to the error estimates, and more TROPOMI data are required to further constrain the weekly cycle. We found that the CO values on Saturday and Sunday are equally low. An explanation for this could be that the main source of CO in Mexico City during the week is traffic which is responsible for the weekly cycle and the remaining sources of cooking, water heating, etc. should not change much during the weekend.

Figure 9 Weekly cycle of the CO concentration. (a) Based on 29 in situ measurement station operated by SEDEMA. (b) Ground-based FTIR vertical column measurements of an instrument located in Mexico. Median values are shown, and the error bars are the standard error in the mean deploying the delta percentile as a robust estimation of the standard deviation.

A similar weekly cycle is observed by Mexico City in situ measurements provided by 29 SEDEMA ground stations. For each of the sites, we use data from 2017 to 2018 for the overpass time of TROPOMI (12:00–15:00 local time), calculate a weekly cycle, and group the data into the stations located in the CdMx urban area and those located in the wider area of the metropolis. Figure 9a depicts the median of all weekly cycles and the standard error in the mean with a clear minimum during weekends. The error bars indicate that the overall shape of the weekly cycles for the remaining days varies a lot from station to station.

The lower CO concentrations during the weekend are also detectable with column retrievals from ground-based FTIR measurements in Mexico City, 2280 m a.s.l.; 19.32^∘ N, −99.18^∘ E, at the campus of the National Autonomous University of Mexico by the Center of Atmospheric Sciences (CCA). The spectra used are recorded in the mid-infrared with a resolution of 0.075 cm⁻¹ (Bezanilla et al., 2014; Plaza-Medina et al., 2017), and the CO column and profile is retrieved using the standard NDACC retrieval strategy (García-Franco et al., 2018; Borsdorff et al., 2018a). Figure 9b shows the averaged weekly cycle with the standard error derived from the column measurements. Due to the low data density at weekends we used the full time range from 5 December 2010 to 10 September 2019 without filtering for the overpass time of TROPOMI. These independent ground-based measurements confirm the weekly CO cycle found in the TROPOMI data. In general, the TROPOMI CO data product agrees very well with retrievals from ground-based FTIR measurements performed by the TCCON network worldwide with an averaged bias of less than 6±3.8 ppb, and the bias with retrievals from NDACC measurements is even lower (Borsdorff et al., 2018a).

5 Conclusions

In this study, we analyzed TROPOMI CO retrieval from 551 overpasses of the instrument over central Mexico, which corresponds to about 2 years of measurements starting from the 14 November 2017 until the 25 August 2019. We found that urban pollution can be monitored by the TROPOMI CO data. The high signal-to-noise ratio of the measurements allowed us to distinguish distinct CO enhancements over the various urban districts of central Mexico using single-orbit overpasses of TROPOMI with a high spatial resolution of 7×7 km² that is enhanced to 7×5.6 km² from 6 August 2019 onwards.

With a dedicated WRF-Chem tracer simulation for the full-time range of the current TROPOMI data record, we could distinguish the contribution of 10 urban districts: Tula, Pachuca, Tulancingo, Toluca, Cuernavaca, Cuautla, Tlaxcala, Puebla, CdMx, and ACdMx. The model data were collocated with the TROPOMI measurements and convolved with the total column averaging kernel to account for the vertical sensitivity of the instrument. Here, the WRF-Chem tracer simulation does not account for atmospheric chemistry, and so the simulated CO tracer fields are linear in the emission rates of the different districts. The model is extended by two effective parameters describing a spatially constant CO background and a dependency of the simulated column on terrain height.

The CO emissions are determined in two steps. First, we apply an unregularized least-squares fit of the model to the TROPOMI observations to determine the averaged emissions per district. A strict data screening based on the measurements and WRF-Chem model simulation reduced the TROPOMI dataset from 551 to 199 overpasses. Second, we solve a regularized least-squares problem, which minimizes the variation in the emissions around their mean to reduce the noise propagation in the inversion. By means of appropriate regularization parameters, we reduce the scatter of the retrieved emissions to about 60 % of the median for all urban districts. For data interpretation and screening, the use of the averaging kernel is of great advantage. The final retrieval product includes an averaging kernel as a retrieval diagnostic, which allows us to analyze retrieval sensitivities and cross-correlations between the inferred emission rates.

The derived averaged emissions for the various urban districts of Mexico deviate significantly from emission estimates of the Inventario Nacional de Emissions de Contaminantes Criterio (INEM) inventory adapted to the period 2017–2019. The TROPOMI emissions from the urban districts Tula (0.10±0.004 Tg yr⁻¹) and Pachuca (0.09±0.005 Tg yr⁻¹) in the north of Mexico City deviate significantly from the INEM inventory with 0.008 Tg yr⁻¹ for both areas. For the emissions of the central part of Mexico City (CdMx), TROPOMI indicates 0.14±0.006 Tg yr⁻¹ versus 0.25 Tg yr⁻¹ INEM emissions, and for the ACdMx district, TROPOMI indicates 0.28±0.01 Tg yr⁻¹ versus 0.14 Tg yr⁻¹ INEM emissions. Together, both districts have similar emissions with 0.42 Tg yr⁻¹ seen by TROPOMI versus 0.39 Tg yr⁻¹ from the inventory, pointing to a different relative distribution of the CO emissions seen by TROPOMI. Moreover, using a posteriori data screening to optimize data selection per emission source allows us to distill a weekly cycle of CO emissions at the districts CdMx and ACdMx from the dataset with a clear minimum during weekends. This finding is in agreement with in situ observations and ground-based FTIR measurement in the metropolis.

Our study shows the need for and the potential of regional atmospheric transport modeling for the interpretation of TROPOMI CO data over metropolitan areas like Mexico City. Here, the CO pollution is a composite of emissions from different districts, and its transport leads to complex CO enhancement patterns in the atmosphere. The WRF-Chem tracer model could simulate the TROPOMI measurements to a great extent; however model errors are still significant and further improvement is required to fully explore the TROPOMI CO observations of urban sources. Another potential error source of our method is the accuracy of the week-to-week and monthly variations in the emissions in the INEM inventory considering the fixed overpass time of TROPOMI. Furthermore, basin cities can be problematic with low wind speed for days, which could lead to the accumulation of signals from more than 1 d in the basins, and this is not yet covered by our approach. To account for this effect in our inversion needs major adjustments, which will be investigated in follow-up studies.