Estimation of OH in urban plumes using TROPOMI-inferred NO₂ ∕ CO

A new method is presented for estimating urban hydroxyl radical (OH) concentrations using the downwind decay of the ratio of nitrogen dioxide over carbon monoxide column-mixing ratios (${\mathrm{XNO}}_{\mathrm{2}}/\mathrm{XCO}$) retrieved from the Tropospheric Monitoring Instrument (TROPOMI). The method makes use of plumes simulated by the Weather Research and Forecast model (WRF-Chem) using passive-tracer transport, instead of the encoded chemistry, in combination with auxiliary input variables such as Copernicus Atmospheric Monitoring Service (CAMS) OH, Emission Database for Global Atmospheric Research v4.3.2 (EDGAR) NO_x and CO emissions, and National Center for Environmental Protection (NCEP)-based meteorological data. NO₂ and CO mixing ratios from the CAMS reanalysis are used as initial and lateral boundary conditions. WRF overestimates NO₂ plumes close to the center of the city by 15 % to 30 % in summer and 40 % to 50 % in winter compared to TROPOMI observations over Riyadh. WRF-simulated CO plumes differ by 10 % with TROPOMI in both seasons. The differences between WRF and TROPOMI are used to optimize the OH concentration, NO_x, CO emissions and their backgrounds using an iterative least-squares method. To estimate OH, WRF is optimized using (a) TROPOMI ${\mathrm{XNO}}_{\mathrm{2}}/\mathrm{XCO}$ and (b) TROPOMI-derived XNO₂ only.

For summer, both the ${\mathrm{NO}}_{\mathrm{2}}/\mathrm{CO}$ ratio optimization and the XNO₂ optimization increase the prior OH from CAMS by 32 ± 5.3 % and 28.3 ± 3.9 %, respectively. EDGAR NO_x and CO emissions over Riyadh are increased by 42.1 ± 8.4 % and 101 ± 21 %, respectively, in summer. In winter, the optimization method doubles the CO emissions while increasing OH by ∼ 52 ± 14 % and reducing NO_x emissions by 15.5 ± 4.1 %. TROPOMI-derived OH concentrations and the pre-existing exponentially modified Gaussian function fit (EMG) method differ by 10 % in summer and winter, confirming that urban OH concentrations can be reliably estimated using the TROPOMI-observed ${\mathrm{NO}}_{\mathrm{2}}/\mathrm{CO}$ ratio. Additionally, our method can be applied to a single TROPOMI overpass, allowing one to analyze day-to-day variability in OH, NO_x and CO emission.

1 Introduction

The rapidly growing urbanization has led to an increase in the number of big cities globally. More than 55 % of the global population resides in cities, and this fraction is projected to increase to 68 % in 2050 (United Nations, 2019). The associated rise in consumption of energy and materials leads to severe air pollution that is estimated to have caused premature deaths of 4 to 9 million people globally in 2015 (Sicard et al., 2021; Pascal et al., 2013; Burnett et al., 2018). Air pollution control measures and the application of cleaner technology reduced the NO₂ concentrations in developed cities such as Los Angeles and Paris by 1.5 % to 3.0 % yr⁻¹ between 1996 and 2017 (Georgoulias et al., 2019). The CO emission was reduced by 28.8 % to 60.7 % in these cities in the period 2000 to 2008 (Dekker et al., 2017). In developing cities such as Tehran and Baghdad, however, NO₂ concentrations increased by 8.6 % yr⁻¹ and 16.9 % yr⁻¹ between 1996 and 2017 (Georgoulias et al., 2019). The CO emission increased by 15 % in New Delhi in the period 2000 to 2008 (Dekker et al., 2017). As a consequence, air pollution monitoring and mitigation in developing cities are becoming increasingly important priorities.

Nowadays, urban air pollution can be studied using a combination of ground-based measurement networks and satellite observations (Sannigrahi et al., 2021; Ialongo et al., 2020). Satellite observations have helped to investigate urban air pollution, particularly in cities without a ground-based monitoring network (Beirle et al., 2019; Borsdorff et al., 2019). In past decades, improvements in the quality and spatial resolution of satellite measurements have allowed the detection of trends in air pollutants and the quantification of urban emissions (Lorente et al., 2019; Verstraeten et al., 2018; Borssdorff et al., 2019). Several studies have focused on NO_x, using NO₂ observations from the SCanning Imaging Absorption spectroMeter for Atmospheric CartograpHY (SCIAMACHY), the Ozone Monitoring Instrument (OMI) and TROPOMI (Ding et al., 2017; Lorente et al., 2019). At the resolution and sensitivity of TROPOMI, urban NO₂ enhancements can be detected readily, even in a single satellite overpass. OMI-derived NO₂ data have been used to quantify NO_x emissions, as well as the urban lifetime of NO₂, as demonstrated by Beirle et al. (2011) using the exponentially modified Gaussian function fit (EMG) method.

In the EMG method, the satellite-observed exponential decay of NO₂ downwind of the city center is used to quantify the first-order loss of NO₂, which is used to quantify the hydroxyl radical (OH), neglecting other NO_x removal pathways. Liu et al. (2016) modified the EMG method for application to complex emission patterns. The quantification of CO emissions from cities is more complicated compared with NO₂ because of its longer lifetime and the related importance of CO sources from the surroundings of cities. Nevertheless, a few studies have demonstrated the feasibility of quantifying relative changes in urban CO emission using Measurement of Pollution in the Troposphere (MOPPIT), Infrared Atmospheric Sounding Interferometer (IASI), Atmospheric Infrared Sounder (AIRS) and TROPOMI observations (Borsdorff et al., 2019; Dekker et al., 2017; Pommier et al., 2013).

In recent years, methods have been developed that combine satellite measurements of different trace gases (for example, the combined use of NO₂ and CO) to obtain specific information about pollutant sources (Lama et al., 2020; Hakkarainen et al., 2016; Miyazaki et al., 2017; Reuter et al., 2019; Silva and Arellano, 2017). The emission factors of CO and NO_x from fuel combustion are uncertain and vary strongly with the combustion efficiency (Flagan and Seinfeld, 1988). The satellite-observed ΔNO₂ $/$ ΔCO ratio is particularly sensitive to this fuel-burning efficiency, as demonstrated by Lama et al. (2020), and can be used to evaluate emission inventories. However, another important uncertainty arises from the removal of NO₂ by OH. OH is an important oxidant in the atmosphere, which determines the lifetime of trace gases such as CO, NO_x, sulfur dioxide (SO₂) and volatile organic compounds (VOCs) (Monks et al., 2009). OH plays an important role in atmospheric chemistry on scales ranging from urban air pollution to the global residence times of greenhouse gases. The direct measurement of OH is possible using spectroscopic methods, but the spatial representativeness of the data is limited due to its short lifetime (de Gouw et al., 2019). OH estimates from global chemical transport models (CTMs) have an uncertainty of > 50 % (Huijnen et al., 2019). Urban measurement campaigns point to large discrepancies between modeled and observed OH abundances – for example, in Lu et al. (2013), who found a factor 2.6 difference in a campaign in the suburbs of Beijing.

The aim of this study is therefore to estimate the average OH concentration in the urban plume of large cities (hereafter referred to as urban OH) from the downwind decay of the TROPOMI-observed ${\mathrm{NO}}_{\mathrm{2}}/\mathrm{CO}$ ratio. The proposed method makes use of the Weather Research and Forecast (WRF) model (Grell et al., 2005) to simulate the meteorological fields and atmospheric transport. The TROPOMI instrument (Veefkind et al., 2012), launched on 13 October 2017 on board the Sentinel-5 Precursor satellite, is particularly well suited for this task, as it measures both compounds with high sensitivity and spatial resolution. Our method uses CO because it has a longer lifetime than NO₂ (weeks to months compared to a few hours). Therefore, CO can be considered as an inert tracer at the timescale of urban plumes. The difference in the rate of decay between NO₂ and CO therefore provides information about the photochemical oxidation of NO₂ because atmospheric dispersion is expected to have a very similar impact on both tracers and therefore cancels out in their ratio. The use of the ${\mathrm{NO}}_{\mathrm{2}}/\mathrm{CO}$ ratio for estimating urban-scale OH is further compared to the EMG method using only satellite-retrieved NO₂ (Beirle et al., 2011).

The city of Riyadh (24.63^∘ N, 46.71^∘ E) is chosen as a test case. Riyadh is an isolated city and a strong source of CO and NO₂ pollution (Beirle et al., 2019; Lama et al., 2020). The frequent clear-sky conditions over Riyadh yield a large number of valid TROPOMI CO and NO₂ data. The signal to noise in TROPOMI is high enough to detect the enhancement of CO and NO₂ over Riyadh in a single overpass (Lama et al., 2020). Model results from the Copernicus Atmospheric Monitoring Service (CAMS) for Riyadh show a distinct seasonality in OH (see Fig. S1 in the Supplement), which we attempt to evaluate using TROPOMI data for summer and winter.

This paper is organized as follows: Sect. 2 describes the TROPOMI NO₂ and CO data, the WRF model setup that was used and the optimization method that is used for estimating OH. Optimization results and comparisons between TROPOMI and WRF are presented in Sect. 3, followed by a summary and conclusion of the main findings in Sect. 4. Additional figures and information about the optimization method are provided in the Supplement.

2 Data and method

2.1 TROPOMI NO₂ tropospheric column

We used the offline TROPOMI level-2 tropospheric column NO₂ [mol m⁻²] data from retrieval versions 1.2.x for 2018 and 1.3.x for 2019, available at https://s5phub.copernicus.eu; http://www.tropomi.eu (last access: 21 September, 2020). NO₂ data of versions 1.2.x and 1.3.x have minor processing differences such as removal of negative cloud fraction, better flagging and uncertainty estimation. However, they use the same retrieval algorithm applied to level-1b version 1.0.0 spectra (Babic et al., 2019) recorded by the TROPOMI UV-Vis module in the 405–465 nm spectral range. The TROPOMI NO₂ differential optical absorption spectroscopy (DOAS) software, developed at KNMI, is used for the processing of NO₂ slant column densities (van Geffen et al., 2019). The improved NO₂ DOMINO algorithm of Boersma et al. (2018) has been used to translate slant columns into tropospheric column densities. In this algorithm, stratospheric contributions are subtracted from the slant column densities, and the residual tropospheric slant column density is converted to tropospheric vertical column density using the air mass factor (AMF). The AMF depends on the surface albedo, terrain height, cloud height, cloud fraction and a priori NO₂ profiles from the TM5-MP model at $\mathrm{1}{}^{\circ }\times \mathrm{1}{}^{\circ }$ (Eskes et al., 2018; Lorente et al., 2017). The comparison of multi-axis differential optical absorption spectroscopy (MAX-DOAS) ground-based measurements in European cities shows that TROPOMI underestimates NO₂ columns by 7 % to 29.7 % (Lambert et al., 2019). To reduce the differences between satellite and model, we re-calculated the AMF by replacing the tropospheric AMF based on TM5-simulated vertical NO₂ columns with the WRF-Chem equivalent (Lamsal et al., 2010; Boersma et al., 2016; Visser et al., 2019; Huijnen et al., 2010) using the equation provided in Appendix A. After the AMF recalculation, the NO₂ vertical profiles are consistent between satellite and model. Furthermore, the use of WRF-Chem has the advantage that it resolves NO₂ gradients between urban and downwind regions better than the coarser-resolution TM5-MP model (Russell et al., 2011; McLinden et al., 2014; Kuhlmann et al., 2015). During summer, the AMF recalculation increases TROPOMI NO₂ by 5 % to 10 % and by 25 % to 30 % in winter in the urban plume over Riyadh, whereas background areas are less affected (see Fig. S2). The Sentinel-5P Product Algorithm Laboratory (S5P-PAL) reprocessed NO₂ data available at https://data-portal.s5p-pal.com/products/no2.html (last access: 1 September 2022) differ by 7.5 % to 10 % in summer (June to October 2018) and 13.5 % to 16 % in winter (November 2018 to March 2019) compared to the AMF-recalculated TROPOMI NO₂ data used in this study. These differences have been used to quantify the systematic uncertainty of the NO₂ data and its contribution to the uncertainty in the NO_x emission and lifetime derived using our method (see Tables S1, S2 and S3 in the Supplement).

2.2 TROPOMI CO

For CO, the offline level-2 CO data product version 1.2.2 has been used, available at https://s5phub.copernicus.eu (last access: 20 September 2020). The Shortwave Infrared Carbon Monoxide Retrieval (SICOR) algorithm is applied to TROPOMI 2.3 µm spectra to retrieve CO total column density [molec. cm⁻²] (Landgraf et al., 2016). The retrieval method is based on a profile-scaling approach, in which TROPOMI-observed spectra are fitted by scaling a reference vertical profile of CO using the Tikhonov regularization technique (Borsdorff et al., 2014). The reference CO profile is obtained from the TM5 transport model (Krol et al., 2005). The averaging kernel (A) quantifies the sensitivity of the retrieved total CO column to variations in the true vertical profile (ρ_true) as follows Borsdorff et al., 2018a):

$$\begin{array}{}\text{(1)}& C_{\mathrm{retrieval}}=A\cdot {\mathit{\rho }}_{\mathrm{true}}+{\in }_{\mathrm{CO}},\end{array}$$

where C_retrieval is the retrieved column-average CO mixing ratio, and ∈_CO is the retrieval error, statistically represented by the retrieval uncertainty that is provided for each CO retrieval.

The comparison of TROPOMI-derived XCO to the 28 different Total Carbon Column Observing Network (TCCON) ground-based stations suggests that the difference between TCCON and TROPOMI is in the range of 9.1 ± 3.3 % (Shah et al., 2020). Such difference is used to estimate the uncertainty in the emission and lifetime (see Tables S1, S2, S3 and Text S6).

2.3 Satellite data selection and filtering criteria

As NO₂ and CO are retrieved from different channels of TROPOMI using different retrieval algorithms, the filtering criteria and spatial resolutions of CO and NO₂ are different. The data filtering makes use of the quality assurance value (qa) and is provided with the CO and NO₂ retrievals, ranging from 0 (no data) to 1 (high-quality data). We selected NO₂ retrievals with qa ≥ 0.75 (clear-sky condition) and CO retrievals with qa ≥ 0.7 (clear sky or low-level cloud), as in Lama et al. (2020). The SICOR algorithm was originally developed for SCIAMACHY to account for the presence of low-elevation clouds, increasing the number of valid measurements (Borsdorff et al., 2018a). In addition, the CO stripe-filtering technique is applied as described by Borsdorff et al. (2018a). Using dry-air column density derived from the surface pressure data in CO and NO₂ TROPOMI files, the total CO column and tropospheric NO₂ column densities are converted to dry column-mixing ratios XCO (ppb) and XNO₂ (ppb). The spatial resolution of the NO₂ data is finer compared to the CO data (3.5×7 km² versus 5.5×7 km²). After the CO and NO₂ retrievals pass the filtering criteria, their co-location is approximated by assigning the center coordinates of an NO₂ retrieval to the CO footprint in which it is located (Lama et al., 2020).

2.4 Weather Research Forecast model (WRF)

We have used WRF chemistry model (http://www.wrf-model.org/, last access: 22 August 2019), version 3.9.1.1 to simulate NO₂ and CO mixing ratios over Riyadh. WRF is a non-hydrostatic model designed by the National Center for Environmental Protection (NCEP) for both atmospheric research and operational forecasting applications. For this study, we have set up three nested domains in the model at resolutions of 27, 9 and 3 km, centered at 24.63^∘ N, 46.71^∘ E. The first and second domain cover Saudi Arabia and provide the boundary conditions for the nested third domain (see Fig. S3). The analysis in this paper uses the 500×500 km² sub region around Riyadh in the third domain, containing 161 by 161 grid cells. All domains are extended vertically from the Earth's surface to 50 hPa, using 31 vertical layers, with 17 layers in the lowermost 1500 m. WRF simulations are performed using a time step of 90 s for the period June 2018 to March 2019, using a spin-up time of 10 d.

We have used the Unified Noah land surface model for surface physics (Ek et al., 2003; Tewari et al., 2004), an updated version of the Yonsei University (YSU) boundary layer scheme (Hu et al., 2013) for the boundary layer processes, and the rapid radiative transfer method (RRTM) for shortwave and longwave radiation (Mlawer et al., 1997). Cloud physics is solved with the new Tiedtke cumulus parameterization scheme (Zhang and Wang, 2017). The WRF Single-Moment 6-class scheme is used for microphysics (Hong et al., 2006). The WRF coupling with chemistry (WRF-Chem) allows the simulation of tracer transport and the chemical transformation of trace gases and aerosols. Here, we used the passive tracer transport function instead of the encoded chemistry in WRF to speed up the model simulation and to reduce the computational cost. In addition, the passive-tracer option helps in separating the influences of wind, OH and the rate constant of the NO₂+OH reaction (K${}_{{\mathrm{NO}}_{\mathrm{2}}\cdot \mathrm{OH}}$) on the ${\mathrm{NO}}_{\mathrm{2}}/\mathrm{CO}$ ratio in the downwind city plume. Compared to previously used methods (Beirle et al., 2011; Valin et al., 2013) which did not use a transport model at all, we consider this an important improvement. The function of different tracers, their acronyms and an explanation of different WRF simulations are provided in Table 1.

The meteorological initial and boundary conditions are based on NCEP data at $\mathrm{1}{}^{\circ }\times \mathrm{1}{}^{\circ }$ spatial and 6 h temporal resolutions, available at https://rda.ucar.edu/datasets/ds083.2/ (last access: 10 September 2019). Nitrogen oxides (NO_x = NO₂+NO) and CO anthropogenic emissions have been taken from the Emission Database for Global Atmospheric Research v4.3.2 (EDGAR) 2012 at $\mathrm{0.1}{}^{\circ }\times \mathrm{0.1}{}^{\circ }$ spatial resolution (Crippa et al., 2016). The EDGAR 2012 data have been re-gridded to the resolution of the WRF domains, and hourly, weekly and monthly emission variations are taken into account using the temporal emission factors provided by van der Gon et al. (2011). The chemical boundary conditions for CO and NO_x are based on the CAMS chemical reanalysis product at $\mathrm{0.75}{}^{\circ }\times \mathrm{0.75}{}^{\circ }$ spatial and 3 h temporal resolutions (Inness et al., 2019), retrieved from https://ads.atmosphere.copernicus.eu/cdsapp#!/dataset/cams-global-reanalysis-eac4?tab=_form, last access: 1 November 2020). XCO and XNO₂ boundary conditions based on CAMS are assumed to be representative as background values within the domain. Since we do not explicitly compute the sources and sinks of background NO₂ inside the domain, we decide to transport the boundary conditions as background passive tracers.

Table 1 Summary of WRF simulations and the definition of tracers and acronyms used.

The atmospheric transport in WRF causes the influence of NO_x and CO emissions from Riyadh on their column-average mixing ratios to be linear. Instead of a simplified photochemistry solver, we make use of a WRF-Chem module for passive tracer transport for transporting NO_x. This WRF module has been modified to account for the first order loss of NO_x in the reaction of NO₂ with OH, using ${\mathrm{NO}}_x/{\mathrm{NO}}_{\mathrm{2}}$ ratios from CAMS to translate NO_x into NO₂ and CAMS OH fields to compute the chemical transformation of NO₂ to HNO₃ (see Text S1 for detail).

This is a simplified treatment of the lifetime of NO_x as other photochemical pathways play a role, such as

• the oxidation of NO₂ in reaction with organic radicals (RO₂) to form the alkyl and multifunctional nitrates (RONO₂) (Romer Present et al., 2019);

• NO_x loss due to the formation of dinitrogen pentoxide (N₂O₅) followed by heterogeneous transformation to HNO₃ (Shah et al., 2020);

• peroxyacetyl nitrate (PAN) formation in equilibrium between NO₂ and the peroxyacetyl radical (Moxim, 1996);

• the dry deposition of NO₂ on the surface and plant stomata (Delaria et al., 2020).

The loss of NO₂ by OH to HNO₃ accounts for 60 % of the global NO_x emission (Stavrakou et al., 2013). Macintyre and Evans (2010) showed that the N₂O₅ pathway reduces NO_x concentrations by 10 % in the tropics (30^∘ N to 30^∘ S) and 40 % at northern latitudes. The NO_x loss through N₂O₅ hydrolysis is largest at northern latitudes during winter (50 % to 150 %), unlike in the tropics where, its seasonality is small. Moreover, the removal of N₂O₅ is primarily important during nighttime because of its photolysis during daytime, whereas our analysis focuses on the midday overpass time (13:30 LT) of TROPOMI when OH abundances are highest. For these reasons, we consider it safe to neglect the loss of NO_x through N₂O₅ in our analysis for Riyadh. The dry deposition flux is also expected to be low, as it is controlled largely by stomatal uptake, which is assumed to be insignificant for the low vegetation cover of Riyadh. The same is expected to be true for PAN formation because of its thermal decomposition at increasing temperatures. We acknowledge that our OH estimates should be regarded as upper limits due to the neglect of other NO_x transformation pathways. A quantification of the combined effect would require full chemistry simulations, which we consider to be outside of the scope of this paper.

Note that, in this study, OH is only applied to the urban NO_x emission tracer (XNO_x emis). The CAMS NO_x background tracer (XNO_x Bg) is transported in WRF without OH decay, since it already represents the balance between regional sources and sinks. CAMS hydroxyl radical (OH) data at $\mathrm{0.75}{}^{\circ }\times \mathrm{0.75}{}^{\circ }$ spatial and 3 h temporal resolutions (Inness et al., 2019), retrieved at https://ads.atmosphere.copernicus.eu/cdsapp#!/dataset/cams-global-reanalysis-eac4?tab=_form, last access: 1 July 2020), are spatially, temporally and vertically interpolated to the WRF grid. The NO_x lifetime is derived as follows:

$$\begin{array}{}\text{(2)}& {\displaystyle \frac{{\mathrm{dNO}}_{\mathrm{2}}}{\mathrm{d}t}}={\mathrm{K}}_{{\mathrm{NO}}_{\mathrm{2}}\mathrm{OH}}\cdot \left[\mathrm{OH}\right]\cdot \left[{\mathrm{NO}}_{\mathrm{2}}\right],\text{(3)}& {\displaystyle }\mathrm{fact}={\displaystyle \frac{{\mathrm{NO}}_x}{{\mathrm{NO}}_{\mathrm{2}}}},\text{(4)}& {\displaystyle }{\mathit{\tau }}_{{\mathrm{NO}}_x}={\displaystyle \frac{\mathrm{1}}{\frac{{\mathrm{K}}_{{\mathrm{NO}}_{\mathrm{2}}\mathrm{OH}}}{\mathrm{fact}}\cdot \left[\mathrm{OH}\right]}},\end{array}$$

where, K${}_{{\mathrm{NO}}_{\mathrm{2}}\mathrm{OH}}$ is the International Union of Pure and Applied Chemistry's (IUPAC's) second order rate constant for the reaction of NO₂ with OH; “fact” represents the fractional contribution of NO₂ to NO_x (${\mathrm{NO}}_x/{\mathrm{NO}}_{\mathrm{2}}$). This NO_x-to-NO₂ conversion factor is derived from the CAMS reanalysis and is re-gridded to WRF to account for its spatial and temporal variation. ${\mathit{\tau }}_{{\mathrm{NO}}_x}$ is the lifetime of NO_x.

Figure 1 TROPOMI-derived XCO (a) and average wind speed and wind direction from the surface to the top of the boundary layer (b) derived from the CAMS global reanalysis eac4 data at the TROPOMI overpass time over Riyadh for 4 August 2018. The white star represents the center of Riyadh. The black box (B1), with a dimension of 300×100 km², is rotated in the average wind direction at a 50 km radius from the center of Riyadh at the TROPOMI overpass time, resulting in the red box. For the calculation of cross-directional averaged NO₂ and CO, the red box is divided into 29 smaller cells with the width (Δx) ∼ 11 km. For this, TROPOMI-derived XCO is gridded at $\mathrm{0.1}{}^{\circ }\times \mathrm{0.1}{}^{\circ }$.

The components of NO_x (NO and NO₂) have short lifetimes during daytime because of the photo-stationary equilibrium exchanging NO and NO₂ into each other. For this reason, we estimate the lifetime of their sum (NO_x), which is determined largely by the reaction with OH. In earlier work with satellite NO₂ data, the Jet Propulsion Laboratory (JPL) high-pressure limit was used as the rate constant to represent the first order loss of NO₂ (Beirle et al., 2011; Lama et al., 2020; Lorente et al., 2019). However, we found this approximation to be too crude and therefore apply the full IUPAC-recommended pressure-dependent formula for the second order rate constant. Supplement Fig. S4 shows the difference between the three rate constants, i.e., JPL high-pressure limit, JPL second order and IUPAC second order, confirming the importance of accounting for the pressure dependence.

WRF output for the third domain is interpolated spatially and temporally to the footprints of TROPOMI. The interpolated WRF NO_x tracers are converted to NO₂ using the conversion factor derived from the CAMS reanalysis, accounting for its spatial and temporal variation (for the names and functions of tracers, see Table 1). The averaging kernel available for each TROPOMI CO and NO₂ observation is applied to the WRF output after interpolation to the vertical layers of the TROPOMI retrieval. To compare WRF output to TROPOMI, WRF-derived XNO₂(XNO_2 WRF) is calculated by combining the NO₂ tracer that accounts for the OH effect (XNO_2 (emis,OH)) and the CAMS NO₂ background (XNO_{2 Bg}) (see Figs. S5 and S6). Similarly, the CO emission tracer (XCO_emis) is added to the CAMS CO background (XCO_Bg) to calculate WRF-simulated XCO (XCO_WRF) (see Figs. S7 and S8).

2.5 ${\mathrm{NO}}_{\mathrm{2}}/\mathrm{CO}$ ratio calculation using box rotation

The variation of the ${\mathrm{NO}}_{\mathrm{2}}/\mathrm{CO}$ ratio in the downwind city plume is calculated as a function of distance x from the city center in a downwind direction. We select days with an average wind speed (U) in the range of 3.0 m s⁻¹ (Beirle et al., 2011) $<U<\mathrm{8}$.5 m s⁻¹ (Valin et al., 2013) within a 50 km radius from the center of Riyadh (24.63^∘ N, 46.71^∘ E). The horizontal distribution of EDGAR emissions over Riyadh is used within this 50 km radius (Fig. S9). A total of 95 d in summer and 70 d in winter meet the wind speed criteria over Riyadh for the ratio calculation. The boundary layer average wind speed and direction are calculated using the CAMS global reanalysis eac4 (retrieved at https://ads.atmosphere.copernicus.eu/cdsapp#!/dataset/cams-global-reanalysis-eac4?tab=_form, last access: 1 August 2020) at a $\mathrm{0.75}{}^{\circ }\times \mathrm{0.75}{}^{\circ }$ spatial and 3 h temporal resolution. For this, the CAMS wind vector is spatially and temporally interpolated to the central coordinate of TROPOMI pixels.

To compute the ${\mathrm{NO}}_{\mathrm{2}}/\mathrm{CO}$ ratio as function of the downwind distance x, TROPOMI and WRF data have been re-gridded at $\mathrm{0.1}{}^{\circ }\times \mathrm{0.1}{}^{\circ }$. A box (B1) is selected with a width of 100 km, from 100 km in upwind to 200 km in downwind direction of the city center (see Fig. 1a). The dimension of the box is motivated by multiple TROPOMI overpasses over Riyadh, showing NO₂ and CO enhancements advected downwind over a ∼ 200 km distance without other large sources of NO₂ and CO within a 100 km radius of the city center (see Fig. 1a). Figure 1b shows the boundary layer averaged wind speed and wind direction over Riyadh, indicating flow towards the northeast on 4 August 2018. The box is rotated for every TROPOMI overpass, depending upon the daily average wind direction within a 50 km radius from the center of Riyadh, as shown in Figs. 1a and S10. The rotated box B1 is divided into N rectangular boxes, orthogonal to the wind direction with length (Δx) ∼ 11 km (see Figs. 1 and S10). The XNO₂ and XCO grid cells that fall within the N rectangular boxes are selected to derive zonally averaged XNO₂ and XCO for summer and winter.

Unlike the enhancements over the city, ΔXNO₂ and ΔXCO become smaller than retrieval uncertainties at large distance from the city, where the ratio ΔXNO${}_{\mathrm{2}}/\mathrm{\Delta }$XCO becomes ill defined. Therefore, we decided to use the ratio of mean XNO₂ and XCO instead of enhancements over the background. To analyze the influence of atmospheric transport and the OH sink on the WRF-derived ${\mathrm{XNO}}_{\mathrm{2}}/\mathrm{XCO}$ ratio, two different ratios are derived: (1) $\frac{{\mathrm{XNO}}_{\text{2\hspace{0.17em}emis}}}{{\mathrm{XCO}}_{\mathrm{emis}}}$, named “ratio_without OH”, and (2) $\frac{{\mathrm{XNO}}_{\text{2 (emis,OH)}}}{{\mathrm{XCO}}_{\mathrm{emis}}}$, named “ratio_{with OH}” (see Table 1). The CAMS background accounts for the balance between regional source and sink in CTMs, so it is excluded to analyze the influence of atmospheric transport on the ratio. For the comparison between TROPOMI and WRF, the CAMS backgrounds are included in “WRF ratio” ($\frac{{\mathrm{XNO}}_{\text{2 WRF}}}{{\mathrm{XCO}}_{\mathrm{WRF}}}$) (see Table 1). The comparison of WRF ratio to TROPOMI ratio and the contribution of its components are presented in Sect. 3.2.

2.6 OH estimation: satellite data only

In the EMG method, following Beirle et al. (2011), 2D NO₂ column density maps are assigned to eight equal wind sectors, spanning 360^∘ for summer and winter; 1D column densities per wind sector are computed by averaging in a cross-wind direction. This way, average NO₂ column density functions of the downwind distance to the city center have been constructed for summer and winter (see Fig. S11). Using the EMG method, as in Beirle et al. (2011), the e-folding distance (x₀) and NO₂ emissions have been estimated. The NO₂ lifetime is derived by dividing x₀ by the average wind speed (5.46 and 5.24 m s⁻¹ for winter and summer, respectively) and is provided in Table 2. The OH concentration is derived from the inferred NO₂ lifetime using the IUPAC second order rate constant (for details, see Text S2 and S3 in the Supplement). Rate constants at the time of TROPOMI overpasses are obtained from WRF by averaging the IUPAC second order rate constant from the surface to the top of the planetary boundary layer (PBL). The PBL height at the time TROPOMI overpass has been taken from WRF. EMG-derived NO₂ emissions are also converted to NO_x emissions using the CAMS-derived conversion factor. Summer- and winter-averaged CAMS-derived conversion factors for the box of 300 km × 100 km are 1.28 and 1.31, respectively.

2.7 OH estimation: WRF optimization

To jointly estimate the NO_x and CO emissions as well as the OH concentration from the TROPOMI data, a least-squares optimization method is used. This method fits the model to the data by minimizing a cost function (J) (see Text S4 for details). The reaction of NO₂ with OH introduces a non-linearity in the OH optimization. To account for this non-linearity, we linearize the problem around the a priori starting point using the small perturbations (10 %) Δbackground, Δemission and ΔOH. The non-linear model is fitted to the observations by optimizing scaling factors f_Bg, f_emis, f_OH to the perturbation functions Δbackground, Δemission and ΔOH, respectively. This process is repeated iteratively, updating the linearization point and re-computing the perturbation functions. The scaling factors f_emis, f_oh and f_bg represent the modification of the prior in percentage change.

We estimate OH by optimizing WRF with TROPOMI in two ways: (1) optimizing the simulated ${\mathrm{NO}}_{\mathrm{2}}/\mathrm{CO}$ ratio using TROPOMI-derived ratios, named as “ratio optimization”, and (2) optimizing NO₂ and CO separately using TROPOMI-derived XCO and XNO₂, named as “component-wise optimization”. First, the ratio optimization is described, followed by the component-wise optimization. Optimized ratios are derived as follows:

$$\begin{array}{}\text{(5)}& {\displaystyle }\begin{array}{rl}{F}_{\mathrm{TROPOMI}}& =F+{\displaystyle \frac{\mathrm{\Delta }F}{\mathrm{\Delta }\mathrm{emis}}}\cdot {\displaystyle \frac{f_{\mathrm{emis}}}{\mathrm{10}}}+{\displaystyle \frac{\mathrm{\Delta }F}{\mathrm{\Delta }\mathrm{OH}}}\cdot {\displaystyle \frac{f_{\mathrm{OH}}}{\mathrm{10}}}\\ & +{\displaystyle \frac{\mathrm{\Delta }F}{\mathrm{\Delta }\mathrm{Bg}}}\cdot {\displaystyle \frac{f_{\mathrm{Bg}}}{\mathrm{10}}},\end{array}\text{(6)}& {\displaystyle }F={\displaystyle \frac{{\mathrm{XNO}}_{\text{2 WRF}}}{{\mathrm{XCO}}_{\mathrm{WRF}}}},\text{(7)}& {\displaystyle }{\mathrm{XNO}}_{\text{2 WRF}}={\mathrm{XNO}}_{\text{2 (emis, OH)}}+{\mathrm{XNO}}_{\text{2 Bg}},\text{(8)}& {\displaystyle }{\mathrm{XCO}}_{\mathrm{WRF}}={\mathrm{XCO}}_{\mathrm{emis}}+{\mathrm{XCO}}_{\mathrm{Bg}},\text{(9)}& {\displaystyle }{\displaystyle \frac{\mathrm{\Delta }F}{\mathrm{\Delta }\mathrm{emis}}}={\displaystyle \frac{{\mathrm{XNO}}_{\text{2 (emis, OH)}}\cdot \mathrm{1.05}+{\mathrm{XNO}}_{\text{2 Bg}}}{{\mathrm{XCO}}_{\mathrm{emis}}\cdot \mathrm{0.95}+{\mathrm{XCO}}_{\mathrm{Bg}}}}-F,\text{(10)}& {\displaystyle \frac{\mathrm{\Delta }F}{\mathrm{\Delta }\mathrm{OH}}}={\displaystyle \frac{{\mathrm{XNO}}_{\mathrm{2}(\mathrm{emis},\mathrm{OH}\times \mathrm{1.1})}+{\mathrm{XNO}}_{\text{2 Bg}}}{{\mathrm{XCO}}_{\mathrm{emis}}+{\mathrm{XCO}}_{\mathrm{Bg}}}}-F,\text{(11)}& {\displaystyle \frac{\mathrm{\Delta }F}{\mathrm{\Delta }\mathrm{Bg}}}={\displaystyle \frac{{\mathrm{XNO}}_{\text{2 (emis, OH)}}+{\mathrm{XNO}}_{\mathrm{2}\mathrm{Bg}}\cdot \mathrm{1.05}}{{\mathrm{XCO}}_{\mathrm{emis}}+{\mathrm{XCO}}_{\mathrm{Bg}}\cdot \mathrm{0.95}}}-F,\end{array}$$

Here, F_TROPOMI is the TROPOMI-derived ${\mathrm{NO}}_{\mathrm{2}}/\mathrm{CO}$ ratio, F is the WRF ratio, $\frac{\mathrm{\Delta }F}{\mathrm{\Delta }\mathrm{emis}}$ is the change in F due to an increase in the NO₂ emission by 5 % and a decrease in the CO emission by 5 % ($\mathrm{1.05}/\mathrm{0.95}=\sim$ 10 %), $\frac{\mathrm{\Delta }F}{\mathrm{\Delta }\mathrm{OH}}$ is the change in F due to an increase in OH by 10 %, and $\frac{\mathrm{\Delta }F}{\mathrm{\Delta }\mathrm{Bg}}$ is the change in F due to an increase in the XNO₂ background by 5 % and a decrease in the CO background by 5 %. XNO_{2 (emis, OH)} is the contribution of city NO_x emissions to XNO₂, accounting for the OH sink. XNO_{2 Bg} is the NO₂ background. XCO_emis is the contribution of the EDGAR city CO emissions to XCO, and XCO_Bg is the CO background derived from CAMS. XNO_{2 WRF} and XCO_WRF are the WRF-derived XNO₂ and XCO, respectively. ${\mathrm{XNO}}_{\mathrm{2}(\mathrm{emis},\mathrm{OH}\times \mathrm{1.1})}$ is the contribution of city NO_x emissions to XNO₂ after increasing CAMS OH by 10 %. The scaling factors f_emis, f_OH and f_Bg obtained from the ratio optimization have been divided by 10 because $\frac{\mathrm{\Delta }F}{\mathrm{\Delta }\mathrm{emis}},\frac{\mathrm{\Delta }F}{\mathrm{\Delta }\mathrm{OH}}$ and $\frac{\mathrm{\Delta }F}{\mathrm{\Delta }\mathrm{Bg}}$ are defined as the change in F due to modification of emission, OH and background by 10 %.

Although the ratio optimization is sensitive to the emission ratio and the OH sink of NO₂, it is not sensitive to the absolute emissions of CO and NO₂. Therefore, we performed component-wise optimizations for XCO and XNO₂ to optimize absolute emissions. We also compare the OH factor obtained from the ratio optimization and component-wise optimization to test the robustness of the method. The optimized XNO₂ is derived using Eq. (12). XCO is optimized using the same equation but without considering the OH sink (see Appendix B).

$$\begin{array}{}\text{(12)}& {\displaystyle }\begin{array}{rl}& {\text{XNO}}_{\text{2 TROPOMI}}={\text{XNO}}_{\text{2 WRF}}+\mathrm{\Delta }X{\text{NO}}_{\text{2 emis}}\cdot {\displaystyle \frac{f_{\text{emis}}}{\mathrm{10}}}\\ & +\mathrm{\Delta }{\text{XNO}}_{\text{2 OH}}\cdot {\displaystyle \frac{f_{\text{OH}}}{\mathrm{10}}}+\mathrm{\Delta }{\text{XNO}}_{\text{2 Bg}}\cdot {\displaystyle \frac{f_{\text{Bg}}}{\mathrm{10}}}\end{array}\text{(13)}& {\displaystyle }\begin{array}{rl}\mathrm{\Delta }{\text{XNO}}_{\text{2 emis}}& ={\text{XNO}}_{\mathrm{2}\left(\text{emis, OH}\right)}\cdot \mathrm{1.10}\\ & -{\text{XNO}}_{\mathrm{2}\left(\text{emis, OH}\right)}\end{array}\text{(14)}& {\displaystyle }\begin{array}{rl}\mathrm{\Delta }{\text{XNO}}_{\text{2 OH}}& ={\text{XNO}}_{\mathrm{2}\left(\text{emis, OH}\times \mathrm{1.1}\right)}\\ & -{\text{XNO}}_{\mathrm{2}\left(\text{emis, OH}\right)}\end{array}\text{(15)}& {\displaystyle }\mathrm{\Delta }{\text{XNO}}_{\text{2 Bg}}={\text{XNO}}_{\text{2 Bg}}\cdot \mathrm{1.10}-{\text{XNO}}_{\text{2 Bg}}\end{array}$$

Here, XNO_{2 TROPOMI} is the TROPOMI-derived XNO₂, XNO_{2 WRF} is the WRF XNO₂. ΔXNO_{2 emis} is the change in XNO₂ due to an increase in emissions by 10 %, ΔXNO_{2 OH} is the change in XNO₂ due to an increase in CAMS OH by 10 %, and ΔXNO_{2 Bg} is a change in the background XNO₂ by 10 %. The scaling factors f_emis, f_OH and f_Bg are divided by a factor 10 because ΔXNO_{2 emis}, ΔXNO_{2 OH} and ΔXNO_{2 Bg} are defined as 10 % changes in NO_x emission, OH and background level.

Figure 2 Comparison between XNO₂ (a, c) and XCO (b, d) from TROPOMI and WRF over Riyadh, averaged over June to October 2018. Panels (a) and (b) show TROPOMI data, and panels (c) and (d) show the corresponding co-located WRF results. XNO_{2 WRF} is derived by adding XNO_{2(emis, OH)} and XNO_{2 Bg}. XCO_WRF is derived by adding XCO_emis and XCO_Bg. The white star represents the center of the city. TROPOMI and WRF results are gridded at $\mathrm{0.1}{}^{\circ }\times \mathrm{0.1}{}^{\circ }$.

3 Results and discussion

3.1 XNO₂ and XCO over Riyadh

In this subsection, we compare WRF-derived XCO_WRF and XNO_{2 WRF} with TROPOMI for summer (see Fig. 2) and winter (see Fig. S12) over Riyadh. TROPOMI and WRF-derived XCO and XNO₂ are averaged from June to October 2018 for summer and November 2018 to March 2019 for winter in a domain of 500×500 km² centered around Riyadh. The comparison for summer in Fig. 2 shows TROPOMI NO₂ after replacing the TM5-based tropospheric AMF with WRF profiles, as described in Visser et al. (2019). The enhancement of XNO₂ and XCO over Riyadh due to urban emissions is clearly separated from the background for TROPOMI and WRF, showing that the city of Riyadh is well suited to investigating the use of the ${\mathrm{NO}}_{\mathrm{2}}/\mathrm{CO}$ ratio to quantify OH in urban plumes. Due to the longer lifetime of CO, the TROPOMI-observed XCO plume extends further in the southeast direction compared to XNO₂. Figure 2 shows that our WRF simulations are able to reproduce the TROPOMI-retrieved XNO₂ (r²=0.96) and XCO (r²=0.78) plumes, confirming that WRF-derived $\frac{{\mathrm{XNO}}_{\text{2 WRF}}}{{\mathrm{XCO}}_{\mathrm{WRF}}}$ is suitable for the optimization of CTM-derived OH concentrations using TROPOMI data. XNO_{2 WRF} is higher by 25 % compared to TROPOMI in the city center. In the background, XCO_WRF shows a similar spatial distribution as TROPOMI XCO, but the values are higher by 5 % to 10 % (see Fig. 2). Close to the city center, XCO_WRF is ∼ 5.7 % higher than TROPOMI XCO. In EDGAR 2011, emission sources are located in the center of Riyadh (see Fig. S9). However, as noted by Beirle et al. (2019), they extend to a larger part of the city in reality. This difference in spatial distribution leads to higher XNO_2WRF and XCO_WRF close to the center of Riyadh compared to TROPOMI.

In winter, the wind direction is predominantly from the southeasterly sector in WRF and TROPOMI (see Fig. S12). The spatial distribution of XCO_WRF (r²=0.73) and XNO_{2 WRF} (r²=0.88) matches quite well with TROPOMI. Therefore, the difference between summer and winter should offer the opportunity to quantify the seasonality in emissions and OH concentrations over Riyadh. In winter, XCO_WRF is ∼ 5 % to 10 % higher than TROPOMI, while XNO_2WRF is higher by 40 % to 50 %. The difference could either point to uncertainties in the ${\mathrm{NO}}_{\mathrm{2}}/\mathrm{CO}$ emission ratio, uncertainties in the NO₂ lifetime or inaccuracies in the background. By quantifying OH, we can evaluate these explanations (see Sect. 3.3). XNO_{2 WRF} is higher by 20 % in winter than in summer. Contrarily, TROPOMI NO₂ is lower by ∼ 30 % in winter (Fig. S12) compared to summer (Fig. 2). Again, to disentangle the role of changing sources and sinks, we need an independent estimate of OH.

3.2 The ${\mathrm{XNO}}_{\mathrm{2}}/\mathrm{XCO}$ ratio and OH

Before comparing TROPOMI- and WRF-derived ${\mathrm{XNO}}_{\mathrm{2}}/\mathrm{XCO}$ ratios, we first analyze the influence of atmospheric transport and the OH sink on the WRF-derived ${\mathrm{XNO}}_{\mathrm{2}}/\mathrm{XCO}$ ratio. To do this, three ratios are used: (1) ratio_{without OH}, (2) ratio_{with OH}, and (3) WRF ratio (see Table 1). As seen in Figs. 3, S13 and S14, WRF is able to reproduce the TROPOMI-observed downwind evolution of XNO₂ and XCO in summer and winter. The peak of the XNO₂ and XCO plumes is shifted away from the city center due to the balance between the accumulation of urban emissions in the atmospheric column and atmospheric transport (Lorente et al., 2019).

As expected, ratio_without OH shows an approximately straight line when the background is removed because transport influences NO₂ and CO in the same way and therefore cancels out in the ratio (see Fig. 3b). The ratio_{with OH} however, shows an approximately Gaussian relation with distance due to the influence of the sink on NO₂. This comparison demonstrates the sensitivity of the relation between the ${\mathrm{XNO}}_{\mathrm{2}}/\mathrm{XCO}$ ratio and downwind distance to the NO₂ lifetime, which we want to exploit to quantify OH. When including the background, the shapes of the functions in Fig. 3c change (not shown) because the relative weights of the background and city contributions to the ratio vary with distance from the city center. In summer, the WRF ratio is higher by ∼ 15 % close to the center-of-city TROPOMI due to the overestimation of XNO_{2 WRF} in WRF (see Fig. 3d). However, in the downwind plume, at a distance of 100 km, WRF ratio is higher by 20 % to 50 % compared to TROPOMI.

Figure 3 Comparison of WRF and TROPOMI averaged across the wind for each small box (a) XNO₂, (b) XCO, (c) WRF ratio (${\mathrm{XNO}}_{\mathrm{2}}/\mathrm{XCO}$) without CAMS background and (d) WRF ratio (${\mathrm{XNO}}_{\mathrm{2}}/\mathrm{XCO}$), with background and TROPOMI as a function of distance to the center of Riyadh for summer (June to October 2018).

In winter, ratio_{without OH} and ratio_{with OH} show relations with downwind distance that are similar to summer, confirming that an OH sink leads to a Gaussian structure of the ratio (see Fig. S14). The winter WRF ratio is 40 % to 60 % higher than TROPOMI due to the overestimation of XNO₂ by 40 % to 50 %. The WRF ratio close to the center of the city is also 20 % higher in winter than in summer due to higher winter XNO_{2 WRF} (see Figs. S12 and S15). In contrast, TROPOMI shows a higher ratio in summer compared to winter (see Fig. S15). These differences between TROPOMI- and WRF-derived ratios offer an opportunity to address uncertainties in CTM-computed urban OH and emission inventories, which will be explored next.

3.3 WRF optimization using synthetic data

To translate the discrepancies between TROPOMI- and WRF-derived ratios of Sect. 3.2 into implied differences in emissions, OH and background, the least-squares optimization method has been used as described in Sect. 2.7. Before optimizing WRF using TROPOMI, pseudo-data experiments in WRF have been carried out to test if the optimization method is capable of recovering true emissions and OH levels. To this end, changes in OH concentrations, emissions and background by known scaling factors have been applied to the WRF prior simulation to create a synthetic dataset. This process is repeated multiple times to create thousands of synthetic datasets. Subsequently, the scaling factors are obtained in the inversion procedure. These tests reveal that the estimation errors for f_emis, f_OH and f_Bg are less than 2.5 % (see Fig. S16). This confirms that the least-squares optimization method works, with two iterations leading to a sufficient accuracy, and can be used to estimate emissions and OH from TROPOMI data. Using TROPOMI data, estimation errors for f_emis, f_OH and f_Bg are expected to be higher due to atmospheric transport errors, simplified chemistry, and XCO and XNO₂ retrieval uncertainties. These errors did not play a role in the pseudo-data experiments, in which perfect transport and sampling were assumed.

To obtain a more realistic estimate of the uncertainty in least-squares-optimized OH, TROPOMI data have been replaced by NO₂, CO and the ${\mathrm{NO}}_{\mathrm{2}}/\mathrm{CO}$ ratio derived from WRF-Chem using the Carbon Bond Mechanism Z (CBM-Z) gas-phase chemical mechanism (Zaveri and Peters, 1999). EDGAR-based VOCs, NO_x and CO emissions have been used in combination with boundary conditions for NO, NO₂, CO and ozone (O₃) from CAMS to run WRF-Chem for 17 August and 18 November 2018, representing a summer and winter day, respectively. For 17 August 2018, the ratio and XNO₂ optimization increase the CAMS-based prior OH of 1.19×10⁷ molec. cm⁻³ by 15.7 % and 13.4 %, respectively (see Fig. S17). In the fully coupled online chemistry with WRF simulation, the boundary layer averaged OH for the box of 300 km × 100 km amounts to 1.33×10⁷ molec. cm⁻³, which is < 5 % lower than the optimized OH value that is derived using our method. The optimized NO_x and CO emissions differ by < 11 % from the emission input used in the full-chemistry-version WRF. In winter, the optimization increases CAMS-based OH of 1.03×10⁷ molec. cm⁻³ by 19.4 %. The OH derived from WRF with full online chemistry is 1.07×10⁷ molec. cm⁻³ and lower by 15.2 % than the optimized OH value. The component-wise optimization increases the EDGAR NO_x and CO emissions by 23.1 % and 10.5 %, respectively (see Fig. S18). Overall, the uncertainty in optimized NO_x, CO emission and OH derived from this test is < 11 % in summer and 10 % to 23 % in winter. Since the lifetime of NO_x is determined by other reactions in addition to the oxidation to HNO₃ considered in our method, it is expected to overestimate the real OH value. The test using WRF full chemistry confirms that this is indeed the case. The uncertainties for OH, NO_x emission and CO emission are in good agreement with the CLASS computations explained in detail in Text S6.

3.4 WRF optimization using seasonally averaged TROPOMI data

The results for summer are summarized in Fig. 4, showing the optimized fit to the TROPOMI data as well as the corresponding scaling factors f_emis, f_OH and f_Bg that are estimated. The optimized emission, OH and Bg obtained from the second iteration are divided by the prior to derive the f_emis, f_OH and f_Bg (see Text S5 for details). The convergence of the iterative procedure is shown in Figs. S19 and S20. The estimated uncertainties for the scaling factors f_emis, f_OH and f_Bg are derived by summing the contribution of wind speed, length and width of the box, NO₂ bias, CO bias, and the different pathways of NO_x loss in quadrature (see Text S6, Tables S1 and S2). For summer and winter, the uncertainties of the optimized OH concentrations are < 17 % and < 29 %, respectively. For NO_x and CO emissions, the uncertainty is < 29 % in summer and winter. Figure 4a shows WRF ratios for summer in comparison to TROPOMI, before and after optimizing the OH concentration. The optimized WRF ratios fit the TROPOMI ratios well, with X²=0.1 (for the derivation of X², see Text S7). The prior and optimized emission ratios, OH concentrations and background ratios obtained from component and ratio optimizations for summer and winter are provided in Table S4. According to the ratio optimization, the emission ratio and CAMS OH are underestimated by 155 ± 26 % and 32 ± 5.3 %, respectively (see Table S4). The optimized CAMS background ratio is lower by 70 ± 6.5 % compared to prior. It should be realized here that the ratio optimization does not estimate the absolute emission of NO₂ and CO, only their ratio.

Figure 4 Comparison between TROPOMI and WRF, before and after optimization for summer (averaged over June to October 2018) (a) ${\mathrm{XNO}}_{\mathrm{2}}/\mathrm{XCO}$ ratio, (b) XNO₂ and (c) XCO in comparison to TROPOMI; f_OH, f_emis and f_Bg are optimized scaling factors obtained iteratively for OH, emissions and background by the least-squares optimization method; f_emis, f_OH and f_Bg are derived by accounting the total change in emission, OH and background using the corresponding scaling factors obtained from the first and second iterative steps. The unit of the scaling factor is in percent (%).

To derive the absolute emission, we performed component-wise optimizations of WRF-derived XCO_WRF and XNO_{2 WRF}. Optimized XCO_WRF and XNO_{2 WRF} fit well to the TROPOMI data (see Fig. 4b and c). In the XNO₂ optimization, the EDGAR NO_x emission is increased by 42.1 ± 8.4 %, and the CAMs background is reduced by 75.9 ± 10.0 %. CAMS OH is increased by 28.3 ± 3.9 %, which is close to the results obtained from the ratio optimization (see Table S4). In the XCO optimization, EDGAR CO emissions are roughly doubled, and the background is reduced by 4.5 ± 0.7 % compared to CAMS (see Table S4).

Figure 5 As Fig. 4, for winter (averaged over November 2018 to March 2019).

The summer optimized ${\mathrm{NO}}_x/\mathrm{CO}$ emission ratio derived from the component-wise optimization is 0.55 ± 0.09. The optimized emission ratio from ratio optimization is larger by factor 3.6 compared to component-wise optimization (see Table S4). The difference between two estimates can be explained by different constraints on the solution in the two methods. In particular, the ratio inversion allows emission adjustment in a fixed relation between NO₂ and CO emissions, whereas the component-wise has the full flexibility to adjust CO and NO₂ emissions. The ${\mathrm{NO}}_{\mathrm{2}}/\mathrm{CO}$ ratio over a city is the sum of the contributions of the background and the city emissions. The relative weight of the two is determined by the absolute background levels and absolute emissions of CO and NO₂. Therefore, the emission ratio estimated by ratio optimization is sensitive to the XNO_{2 Bg}. However, the difference between the two estimates is larger than expected but does not affect the OH estimation. Lama et al. (2020) inferred an ${\mathrm{NO}}_{\mathrm{2}}/\mathrm{CO}$ emission ratio over Riyadh of 0.47 ± 0.1 for 2018 from TROPOMI, favoring the Monitoring Atmospheric Chemistry and Climate and CityZen (MACCity) emission ratio over that of EDGAR. The optimized emission ratio obtained from component-wise optimization is consistent with Lama et al. (2020) and with MACCity summer emissions. This shows that, for the accurate estimation of the emission and emission ratio, the component-wise optimization method is preferable.

Figure 5 presents optimization results for winter, where optimized WRF is in similar good agreement with TROPOMI as for summer, with X²=0.11. For winter, the ratio optimization increases the emission ratio by 58.8±33 % and OH by 52.0 ± 14 %. The ratio and component-wise optimizations again show similar OH adjustments, demonstrating the robustness of our method. The background ratio is reduced by 66.8 ± 11 %. The XNO₂ optimization reduces the EDGAR NO_x emission by 15.45 ± 4.1 % and the CAMS background by 70.2 ± 6.1 %. For XCO, the WRF XCO_Bg is reduced by 1.74 ± 0.1 % in combination with a doubling of the EDGAR CO emission. The optimized emission ratio (${\mathrm{NO}}_x/\mathrm{CO}$) derived from component-wise optimization is 0.36, which is 4 times lower than the optimized emission ratio obtained from ratio optimization (see Table S4).

Table 2 Overview of WRF-optimized OH and NO_x emissions for Riyadh and comparison to the EMG method. The estimated uncertainty for EMG- and WRF-derived NO_x emissions and OH concentrations is the sum of the contribution of wind speed, length and width of the box, NO₂ bias correction, CO bias, and the different pathways of NO_x loss provided in Tables S1, S2 and S3.

3.5 WRF optimization using a single TROPOMI overpass

To demonstrate the application of our WRF optimization method to a single TROPOMI overpass, results are presented in this subsection for 18 August 2018. This date was selected for clear-sky conditions, with most of the TROPOMI NO₂ and CO pixels passing the data quality filter. During this day, the urban plume was transported in a southwestern direction over Riyadh. The spatial distribution of XNO_{2 WRF} (r²=0.76) and XCO_WRF (r²=0.65) matches quite well with TROPOMI (see Fig. S21). The optimized ratio, XNO₂ and XCO for a single day fit well with TROPOMI (X²=0.1, 0.3 and 0.7), comparable to the summer-averaged plumes, indicating that the optimization method can be applied to a single TROPOMI overpass. The ratio optimization increases the emission ratio and CAMS OH by 111 ± 18.4 % and 37.9 ± 6.2 %, respectively, whereas the background is reduced by 51.5 ± 5.2 % (see Fig. S22a). The XNO₂ optimization increases the EDGAR NO_x emission by 25.5 ± 5.1 % and CAMS OH by 32.3 ± 4.4 %, whereas the NO_x background is reduced by 54.4 ± 7.0 % (see Fig. S22b). The CO optimization doubles the EDGAR CO emission and reduces the background by 6.1 ± 0.97 % (see Fig. S20c). The optimized NO_x and CO emissions for 18 August are 8.9 ± 1.7 and 18.9 ± 4.0 kg s⁻¹, respectively, and differ by < 25 % from the summer optimized emission (see Tables 2 and S5). The optimized OH derived from a single TROPOMI overpass is 1.73×10⁷ ± 0.3 molec. cm⁻³ and differs by < 5 % from the summer-averaged OH, i.e., 1.7×10⁷ ± 0.3 molec. cm⁻³, confirming that the method yields realistic results for a single overpass.

3.6 WRF optimization vs. the EMG method

To investigate the consistency between our method and the EMG method, the derived NO_x lifetimes, emissions and OH concentrations using both methods are listed in Table 2 for winter and summer. Our optimization and the EMG method agree well on the seasonal change in NO_x emission and OH concentration. Both methods result in higher NO_x emissions and shorter lifetimes in summer and lower NO_x emissions and longer lifetimes in winter. Riyadh has dry and warm summer days, and the increase in power consumption due to the use of air conditioning contributes to the higher emissions in summer than in winter (Lange et al., 2022). During the summer, EMG and the WRF optimization method both increase the NO_x emission and OH concentration compared with the prior. The sizes of the NO_x emission and OH concentration increase, obtained using the WRF optimization method, which is higher than the EMG method by 10 % to 29 %. However, the differences between the EMG method and the component optimization method are smaller compared to the uncertainties of the emission and OH concentration derived for the optimization method. For winter, the difference between the EMG and WRF-optimized results are smaller than the difference between the EMG results and the prior. The NO_x emission after optimization differs from the EMG method by 33 %. Optimized OH concentration and NO_x lifetime differ by < 10 % compared to the EMG method. In general, the difference between the EMG and optimization results is within the uncertainty range of 20 % to 30 %, confirming their consistency and strengthening the confidence in the estimates that are obtained from TROPOMI data. In contrast to the EMG method, the optimization method can be used for a single TROPOMI overpass (see Sect. 3.5) and does not require yearly averaged NO₂ data, allowing analysis of day-by-day OH, NO_x and CO emissions (see Sect. 3.3). Segregation and averaging of NO₂ urban plumes by wind sector is not required in the optimization method. The effect of transport cancels out in taking the ${\mathrm{NO}}_{\mathrm{2}}/\mathrm{CO}$ ratio, and loss of NO₂ is mostly governed by OH during the midday. In this study, NO_x emission and OH concentration are estimated iteratively, whereas the EMG method arrives at the solution in a single step. However, since our optimization method requires a WRF model simulation, it is computationally more expensive. Uncertainties in transport may create mismatches with the satellite observations, leading to errors in the optimized fit. This influences the quality of derived emission estimates (Dekker et al., 2017). Therefore, finding a simplified approach using satellite data to derive the emission ratio and to estimate OH concentration in urban plumes will be our focus in the future. In the future, the accuracy of our method can be further improved by accounting for other NO_x removal pathways.

3.7 WRF optimized emissions and emission trends

It should be realized that the a priori EDGAR emissions and TROPOMI optimized estimates represent different years (2012 and 2018, respectively). To check whether the emission differences that are found may be explained by trends in emissions, we compare EDGARv5.0 2012 NO_x and CO emissions with 2018, accounting for seasonal and diurnal emission variations using temporal emission factors from van der Gon et al. (2011). EDGAR 2018 NO_x and CO emissions are derived by linear extrapolation using emissions from 2000 to 2015 (see Fig. S23). For summer, midday NO_x emissions, the EDGAR emissions increased by 16.7 % from 2012 to 2018, which is lower than our optimization results. For winter, midday NO_x emissions increased in EDGAR by 15.2 % from 2012 to 2018, whereas the WRF optimization yielded reductions of 15.6 %. In EDGAR, summer and winter CO emissions increased from 2012 to 2018 by 38.5 %. However, the WRF optimization suggests that the EDGAR CO emissions for summer and winter need to be doubled (see Table S5). Borsdorff et al. (2018b) mentioned that EDGAR CO emissions have to be increased significantly to match with TROPOMI CO observations over middle eastern cities such as Tehran, Yerevan, Tabriz and Urmia. Overall, this points to a significant uncertainty in the EDGAR emission inventory at the city scale.

To test the accuracy of the linear extrapolation of EDGAR data, we compare the relative change in NO_x and CO emissions in 2012 to 2018 using CAMS Global (CAMS-GLOB) anthropogenic v4.2 emission datasets (https://ads.atmosphere.copernicus.eu/cdsapp#!/dataset/cams-global-emission-inventories?tab=overview, last access: 15 June 2021). CAMS-GLOB shows that, for summer and winter, NO_x emissions increased by 26 % from 2012 to 2018, which is higher than EDGAR by a factor of 1.7. CAMS-GLOB-based summer and winter CO emissions increase by 20 % from 2012 to 2018, which differs by ∼ 40 % compared to EDGAR. In general, the relative increase in CO and NO_x emissions from EDGAR and CAMS-GLOB is much smaller compared to the difference with our optimization method.

4 Discussion

The TROPOMI-retrieved ${\mathrm{XNO}}_{\mathrm{2}}/\mathrm{XCO}$ ratio is useful for estimating midday OH over isolated localized sources, such as the city of Riyadh, showing a clear contrast between the urban plume and the background. Such TROPOMI-derived OH estimates offer a new opportunity to evaluate urban photochemistry in chemistry transport models. OH depends non-linearly on NO_x and VOC emissions, meteorological conditions, etc. (Sillman et al., 1990), which vary substantially between cities that are monitored by TROPOMI. Therefore, the application of our method to the global and multi-year dataset that is available could contribute substantially to the understanding of urban photochemistry and the development of effective pollution mitigation strategies. In addition, the method requires local sources with NO₂ and CO emissions that are large enough to be detected by TROPOMI. Especially in European cities with lower CO emission, where TROPOMI cannot detect the CO enhancement along with NO₂, this method cannot be applied.

We realize that our method only considers the first order loss of NO₂ by OH-forming HNO₃. In reality, the NO₂ lifetime is influenced by more spatially and temporally varying factors such as temperature, ozone and radiation (Lang et al., 2015; Romer et al., 2018). In cities, the loss of NO₂ via the formation of alkyl and multifunctional nitrates (RONO₂) is also an important reaction influencing the lifetime of NO₂ (Browne et al., 2013; Sobanski et al., 2017). For CO, secondary production from short-lived volatile organic compounds can also play an important role in urban pollution plumes. The application of full chemistry that includes all the sources and losses of NO₂ and CO could therefore further improve the accuracy of OH estimates.

For cities at higher latitudes, especially in winter, it becomes more critical to account for the contribution of other pathways of NO_x loss than OH oxidation. Isolated tropical and subtropical cities are therefore best suited for application of our current method.

A sensitivity test has been performed in which XNO_x, Bg is lost by OH. In this case, the optimized NO_x emission and OH for summer and winter differ by < 7 % from the default method, where the background is treated as an inert tracer (see Table S6). Furthermore, a sensitivity test has been performed in which the prior emission has been changed. The optimized emission varied by < 5 %, demonstrating robustness of the method to the choice of prior (see Fig. S24). This also indicates that the optimization method can be used to study emission changes. Figure S25 shows that road transport, power plants and manufacturing industries are the largest pollutant emitter over Riyadh (Beirle et al., 2019). In this study, NO_x and CO anthropogenic emissions are introduced at the surface, whereas the emission height of different sources is expected to vary in reality. The different emission heights for NO_x and CO emission sources can also influence the result. In the future, realistic emission heights should also be incorporated in WRF for accurate estimation of OH. Moreover, the temporal emission factors that have been used by van der Gon et al. (2011) are based on European countries. The comparison of van der Gon et al. (2011) with the Copernicus Atmosphere Monitoring Service TEMPOral profiles (CAMS-TEMPO) (Guevara et al., 2021) suggests that temporal emission factors for weekend road transport and monthly residential combustion are different in Riyadh compared to European countries. CAMS-TEMPO is expected to provide a more accurate representation of emission variation due to the information on temporal and spatial variations that is included. Road transport CO emissions are the largest contributor by ∼ 75 % to the total emissions over Riyadh, whereas NO_x emissions from the road contribute by 24 % to the total NO_x emission. Residential combustion has the smallest contribution of ∼ 0.3 % to 0.4 % to total NO_x and CO emissions (see Fig. S25). In the future, the application of accurate diurnal emission factors for road transport (see Fig. S26) can further improve the accuracy of urban OH concentrations estimated using TROPOMI-derived ${\mathrm{XNO}}_{\mathrm{2}}/\mathrm{CO}$ ratios. In addition, the seasonality for NO_x and CO emissions is different in Riyadh than in Europe, which should also be accounted for in future studies.

5 Conclusions

In this study, a new method is presented for estimating OH concentrations in urban plumes using TROPOMI-observed ${\mathrm{XNO}}_{\mathrm{2}}/\mathrm{XCO}$ ratios in combination with WRF simulations of the downwind pollution plume of large cities. Our new method has been tested for the city of Riyadh using synthetic as well as real TROPOMI data. Seasonal emissions and OH concentrations have been estimated for summer (June to October 2018) and winter (November 2018 to March 2019). WRF is well able to reproduce the spatial distribution of TROPOMI-retrieved XNO₂ and XCO plumes over Riyadh during the summer and winter seasons. However, before the optimization, WRF overestimates XNO₂ by 15 % to 30 % in summer and 40 % to 50 % in winter compared to TROPOMI. In both seasons, TROPOMI XCO agrees within 10 % with WRF. The WRF-derived ${\mathrm{XNO}}_{\mathrm{2}}/\mathrm{XCO}$ ratio is higher by 15 % to 30 % in summer and 40 % to 60 % in winter compared to TROPOMI, explained mostly by differences in XNO₂.

The differences between WRF and TROPOMI observations have been used to optimize emissions and the NO₂ lifetime. To this end, scaling factors for the city emissions, OH and the background level have been optimized iteratively using a least-squares method. Ratio and component-wise optimizations have been compared to test the overall consistency of the method. In summer, the ratio and XNO₂ optimization for XNO₂ suggest that the OH prior from CAMS is underestimated by 32 ± 5.3 % and 28.3 ± 3.9 %, respectively . The OH estimates obtained from the ratio and NO₂-only optimization differ by < 10 %, demonstrating the robustness of the method. Summertime emissions of NO_x and CO from EDGAR are increased by 42.1 ± 8.4 % and 101 ± 21 %. For winter, the ratio and component-wise optimizations increase OH by ∼ 52 ± 14  % to fit TROPOMI-inferred ratios. In the optimization of winter data, NO_x emissions are reduced by 15.5 ± 4.1 %, and CO emissions are doubled. In the future, the remaining differences between TROPOMI observations and WRF simulations could be reduced further by the use of precise temporal and monthly emission factors, emission heights and full chemistry to account for secondary sources and sinks of CO and NO₂.

TROPOMI-inferred OH concentrations obtained from the least-squares optimization method have been compared to the EMG method. For the summer and winter, the optimized OH concentrations differ by 10 % between the two methods. These results confirm that urban emissions and OH concentrations can be estimated robustly from TROPOMI data. With our method, single TROPOMI overpasses can be used to estimate OH, whereas the EMG method requires averaging of urban NO₂ plumes by wind sector. The iterative approach allows one to test the factors, i.e., f_emis, f_oh and f_bg, obtained from the optimization method, whereas the EMG method does not allow such flexibility.

An important remaining uncertainty is the bias correction of the TROPOMI XNO₂ retrieval. Following the recommended procedure, the air mass factor AMF is recalculated by replacing the tropospheric AMF based on TM5, which is provided with the data, with WRF-Chem. The TROPOMI XNO₂ bias correction increases the mixing ratio in the urban plume of Riyadh by 5 % to 10 % in summer and 25 % to 30 % in winter. The background is less affected by the bias correction. Without TROPOMI XNO₂ bias correction, the uncertainty in the scaling factor for OH can vary up to 20 % and up to 60 % for NO_x emissions over Riyadh.

Appendix A: AMF recalculation

The air mass factor (AMF) used in the retrieval of TROPOMI XNO₂ has been re-calculated by replacing the tropospheric AMF, calculated from the NO₂ column simulated by TM5, with its WRF-Chem equivalent, as described by Lamsal et al. (2010) and Boersma et al. (2016), using the following Eq. (A1):

$$\begin{array}{}\text{(A1)}& M_{\mathrm{trop},\mathrm{WRF}}=M_{\mathrm{trop},\mathrm{TM}\mathrm{5}}\times {\displaystyle \frac{\sum _{l=\mathrm{1}}^L{A}_{\mathrm{trop},l}{x}_{l,\mathrm{WRF}}}{\sum _{l=\mathrm{1}}^L{x}_{l,\mathrm{WRF}}}},\end{array}$$

where M_trop,WRF and M_trop,TM5 are the tropospheric air mass factors derived from WRF and TM5, respectively. A_trop,l is the tropospheric averaging kernel, ranging from the surface to the uppermost layer of the troposphere in the TM5 model (l). x_l,WRF is the equivalent NO₂ column density in model layer l, based on WRF. A_trop in Eq. (A1) is derived using $A_{\mathrm{trop}}=\mathrm{A}\times \frac{M}{M_{\mathrm{trop}}}$, where M and M_trop are the total and tropospheric AMFs, respectively. Finally, the bias-corrected NO₂ vertical column density is computed using

$$\begin{array}{}\text{(A2)}& {\mathrm{NO}}_{\text{2, bias corrected}}={\displaystyle \frac{M_{\mathrm{trop},\mathrm{TM}\mathrm{5}}}{M_{\mathrm{trop},\mathrm{WRF}}}}\times {\mathrm{NO}}_{\mathrm{2}},\end{array}$$

where NO₂ is the TROPOMI tropospheric NO₂ vertical column density, and NO_{2, bias-corrected} is the bias-corrected TROPOMI tropospheric NO₂ vertical column density.

Appendix B: XCO component-wise optimization

The component-wise optimization of XCO_WRF to estimate the emission and background of CO uses the following equations:

$$\begin{array}{}\text{(B1)}& {\displaystyle }\begin{array}{rl}{\mathrm{XCO}}_{\mathrm{TROPOMI}}& ={\mathrm{XCO}}_{\mathrm{WRF}}+\mathrm{\Delta }{\mathrm{XCO}}_{\mathrm{emis}}\cdot {\displaystyle \frac{f_{\mathrm{emis}}}{\mathrm{10}}}\\ & +\mathrm{\Delta }{\mathrm{XCO}}_{\mathrm{Bg}}\cdot {\displaystyle \frac{f_{\mathrm{Bg}}}{\mathrm{10}}},\end{array}\text{(B2)}& {\displaystyle }{\mathrm{XCO}}_{\mathrm{WRF}}={\mathrm{XCO}}_{\mathrm{emis}}+{\mathrm{XCO}}_{\mathrm{Bg}},\text{(B3)}& {\displaystyle }\mathrm{\Delta }{\mathrm{XCO}}_{\mathrm{emis}}=\mathrm{0.10}\cdot {\mathrm{XCO}}_{\mathrm{emis}},\text{(B4)}& {\displaystyle }\mathrm{\Delta }{\mathrm{XCO}}_{\mathrm{Bg}}=\mathrm{0.10}\cdot {\mathrm{XCO}}_{\mathrm{Bg}}.\end{array}$$

Here, XCO_TROPOMI is TROPOMI XCO, XCO_WRF is the WRF-simulated XCO accounting for emissions and background CO, XCO_emis is the XCO contribution from the urban CO emission, and XCO_Bg is the CAMS-derived XCO background. ΔXCO_emis is the change in XCO due to emission, and ΔXCO_Bg is the change in the XCO background level.