Estimation of OH in urban plumes using TROPOMI inferred NO 2 / 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 (XNO 2 / XCO) retrieved from the Tropospheric Monitoring Instrument (TROPOMI). The method makes use of plumes simulated by the Weather Research and Forecast 15 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) NOx and CO emissions, and National Center for Environmental Protection (NCEP) based meteorological data. NO 2 and CO mixing ratios from the CAMS reanalysis are used as initial and lateral boundary conditions. WRF overestimates NO 2 plumes close to the center of the city by 15 % to 30 % in summer and 40 % to 50 % in 20 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, NOx, CO emissions and their backgrounds using a iterative least square method. To estimate OH, WRF is optimized using a) TROPOMI

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 NO2/CO ratio. The proposed method makes use of the 75 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 NO2 (weeks-months compared to a few hours). Therefore, CO can be considered as an inert tracer at the time-scale of urban plumes. The difference in the rate of decay between NO2 and CO provides therefore information about the 80 photochemical oxidation of NO2, 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 NO2/CO ratio for estimating urban scale OH is further compared to the Exponentially Modified Gaussian function fit (EMG) method, using only satellite retrieved NO2 (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 NO2 pollution (Beirle et al., 2019;Lama et al., 2020). The frequent clear sky conditions over Riyadh yield a large number of valid 85 TROPOMI CO and NO2 data. The signal to noise in TROPOMI is high enough to detect the enhancement of CO and NO2 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), which we attempt to evaluate using TROPOMI data for summer and winter. This paper is organized as follows: Section 2 describes the TROPOMI NO2 and CO data, the WRF model setup that was used, 90 and the optimization method that is used for estimating OH. Optimization results and comparisons between TROPOMI and WRF are presented in section 3, followed by a summary and conclusion of the main finding in section 4. Additional figures and information about the optimization method are provided in the Supplement. We used the offline TROPOMI level 2 tropospheric column NO2 [mole m -2 ] 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). NO2 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-465nm spectral range. The TROPOMI NO2 DOAS software, 100 developed at KNMI, is used for the processing of NO2 slant column densities (van Geffen et al., 2019). The improved NO2 DOMINO algorithm of  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 NO2 profiles from TM5-MP at 1 • × 1 • (Eskes et 105 al., 2018;Lorente et al., 2017). The comparison of MAX-DOAS ground based measurements in European cities shows that TROPOMI underestimates of NO2 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 NO2 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 the Appendix A. After the AMF recalculation, the NO2 vertical profiles are consistent between satellite 110 and model. Furthermore, the use of WRF-Chem has the advantage that it resolves NO2 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 NO2 by 5 % to 10 % and in winter by 25 % to 30 % in the urban plume over Riyadh, whereas background areas are less affected (see Fig S2 ). The S5P-PAL reprocessed NO2 data based on version 2.3.1 available at https://data-portal.s5p-pal.com/products/no2.html differs by 7.5 % to 10 % in summer ( June to 115 October, 2018) and 13.5 % to 16 % in winter (November, 2018to March, 2019 compared to the AMF recalculated TROPOMI NO2 data used in this study. These differences have been used to quantify the systematic uncertainty of the NO2 data and its contribution to the uncertainty in the NOx emission and lifetime derived using our method (see Table S1, S2 and S3).

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 120 access: 20 September, 2020). The SICOR algorithm is applied to TROPOMI 2.3 μm spectra to retrieve CO total column density [molec cm −2 ] (Landgraf et al., 2016). The retrieval method is based on a profile scaling approach, in which TROPOMIobserved 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) Formatted: Subscript quantifies the sensitivity of the retrieved total CO column to variations in the true vertical profile (ρ true ), as follows (Borsdorff 125 et al., 2018a): where, Cretrieval is the retrieved column average CO mixing ratio, ∈ 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 TCCON ground based station suggest that difference between 130 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 NOx emission and life time (see Table S1, S2, S3 and Text S6).

Satellite Data Selection and Filtering Criteria
As NO2 and CO are retrieved from different channels of TROPOMI using different retrieval algorithms, the filtering criteria 135 and spatial resolutions of CO and NO2 are different. The data filtering makes use of the quality assurance value (qa) and is provided with the CO and NO2 retrievals, ranging from 0 (no data) to 1 (high quality data). We selected NO2 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 140 by Borsdorff et al. (2018). Using dry air column density derived from the surface pressure data in CO and NO2 TROPOMI files, the total CO column and tropospheric NO2 column densities are converted to dry column mixing ratios XCO (ppb) and XNO2 (ppb). The spatial resolution of the NO2 data is finer compared to the CO data (3.5x7 km 2 versus 5.5x7 km 2 ). After the CO and NO2 retrievals pass the filtering criteria, their co-location is approximated by assigning the centre coordinates of an NO2 retrieval to the CO footprint in which it is located (Lama et al., 2020). 145

Weather Research Forecast model (WRF)
We have used WRF-chemistry model (http://www.wrf-model.org/ ), version 3.9.1.1 to simulate NO2 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 setup three nested domains in the model at resolutions of 27 km, 9 km and 3 km, centred at 24.63°N, 46.71°E. The first and second domain cover Saudi Arabia and 150 provide the boundary conditions for the nested third domain (see Fig. S3). The analysis in this paper uses the 500 x 500 km 2 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 seconds for the period June 2018 to March 2019, using a spin-up time of 10 days. We have used the Unified Noah land surface model for surface physics (Ek et al., 2003;Tewari et al., 2004), an updated 155 version of the Yonsei University (YSU) boundary layer scheme  for the boundary layer processes, and the Rapid Radiative Transfer Method (RRTM) for short-wave and long-wave 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 and Lim, 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 160 transport function instead of the encoded chemistry in WRF to speed up the model simulation and reduce the computational cost. In addition, the passive tracer option helps in separating the influences of wind, OH and the rate constant of the NO2+OH reaction (KNO2.OH) on the NO2/CO ratio in the downwind city plume. Compared to previously used methods (Beirle et al., 2011b;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 acronym and explanation of different WRF simulations is provided in Table 1. 165 The meteorological initial and boundary conditions are based on NCEP data at 1°x1° spatial and 6-hr temporal resolution available at https://rda.ucar.edu/datasets/ds083.2/. Nitrogen Oxides (NOx = NO2 +NO) and CO anthropogenic emissions have been taken from the Emission Database for Global Atmospheric Research v4.3.2 (EDGAR) 2012 at 0.1°x0.1° 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 170 der Gon et al. (2011). The chemical boundary conditions for CO and NOx are based on the CAMS chemical reanalysis product at 0.75°x0.75° spatial, and 3-hourly temporal resolution (Inness et al., 2019), retrieved from https://ads.atmosphere.copernicus.eu/cdsapp#!/dataset/cams-global-reanalysis-eac4?tab=form, last access: 1 st November, 2020). XCO and XNO2 boundary condition based on CAMS is assumed to be representative as background value within the domain. Since we do not explicitly compute the sources and sinks of background NO2 inside the domain, we decide to 175 transport the boundary conditions as background passive tracers. 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 180 transport for transporting NOx. This WRF module has been modified to account for the first order loss of NOx in reaction of NO2 with OH, using NOx/NO2 ratios from CAMS to translate NOx into NO2 and CAMS OH fields to compute the chemical transformation of NO2 to HNO3 (see Text S1 for detail ).
In addition, we account for the chemical transformation of NOx to HNO3 in the reaction of NO2 with OH. This is a simplified treatment of the lifetime of NOx as other photochemical pathways play a role, such as: 185 -The oxidation of NO2 in reaction with organic radicals (RO2) to form the alkyl and multifunctional nitrates (RONO2) (Romer Present et al., 2019) -NOx loss due to the formation dinitrogen pentoxide (N2O5) followed by heterogeneous transformation to HNO3 (Shah et al., 2020).
The loss of NO2 by OH to HNO3 accounts for 60 % of the global NOx emission (Stavrakou et al., 2013). Macintyre and Evans.,(2010) showed that the N2O5 pathway reduces NOx concentrations by 10 % in the tropics (30 o N to 30 o S) and 40 % at northern latitudes. The NOx loss through N2O5 hydrolysis is largest at northern latitudes during winter (50 % to 150 %), unlike the tropics where its seasonality is small. Moreover, the removal of N2O5 is primarily important during night time 195 because of its photolysis during daytime, whereas our analysis focuses on the midday overpass time (13:30) of TROPOMI when OH abundances are highest. For these reasons, we consider it safe to neglect the loss of NOx through N2O5 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 200 as upper limits due to the neglect of other NOx transformation pathways. A quantification of the combined effect would require full chemistry simulations, which we consider outside of the scope of this paper.
Note that in this study, OH is only applied to the urban NOx emission tracer (XNO x emis ). The CAMS NOx 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 a resolution of 0.75° x 0.75° spatial and 3 hourly temporal resolution (Inness 205 et al., 2019) retrieved at https://ads.atmosphere.copernicus.eu/cdsapp#!/dataset/cams-global-reanalysis-eac4?tab=form, last access: 1 st July, 2020) is spatially, temporally and vertically interpolated to the WRF grid. The NOx lifetime is derived as follows: where, KNO2 OH is the International Union of Pure and Applied Chemistry (IUPAC) 2 nd order rate constant for the reaction of NO 2 with OH. "fact" represents the fractional contribution of NO2 to NOX (NOx/NO2). This NOx to NO2 conversion factor is 215 derived from the CAMS reanalysis and re-gridded to WRF, to account for its spatial and temporal variation. τ NOx is the lifetime of NOx.

220
The components of NOx (NO and NO2) have short lifetimes during daytime because of the photo stationary equilibrium exchanging NO and NO2 into each other. For this reason, we estimate the lifetime of their sum (NOX) which is determined largely by the reaction with OH. In earlier work with satellite NO2 data, the Jet Propulsion Laboratory (JPL) high pressure limit was used as rate constant to represent the first order loss of NO2 (Beirle et al., 2011;Lama et al., 2020;Lorente et al., 225 2019). However, we found this approximation to be too crude, and therefore apply the full IUPAC recommended pressure dependent formula for the 2 nd order rate constant. Supplement Figure S4 shows the difference between the three rate constants, i.e. JPL high pressure limit, JPL 2 nd order and IUPAC 2 nd 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 230 WRF-NOx tracers are converted to NO2 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 NO2 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 XNO2 (XNO 2 WRF ) is calculated by combining the NO2 tracer that accounts for the OH effect (XNO 2 (emis,OH) ) and the CAMS NO2 background ( XNO 2 Bg ) (see Fig. S5 and S6) . 235 Similarly, the CO emission tracer (XCO emis ) is added to the CAMS CO background (XCO Bg ) to calculate WRF simulated XCO (XCO WRF ) (see Fig. S7 and S8).

NO2/CO ratio calculation using box rotation
The variation of the NO2/CO ratio in the downwind city plume is calculate as a function of distance x from the city centre in downwind direction. We select days with an average wind speed (U) in the range of 3.0 ms -1 (Beirle et al., 2011) < U < 8.5 240 ms -1 (Valin et al., 2013) within a 50 km radius from the centre 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). Ninety five days in summer and 70 days in winter meet the wind speed criteria over Riyadh for the ratio calculation. The boundary layer average wind speed and direction is calculated using the CAMS global reanalysis eac4 (retrieved at https://ads.atmosphere.copernicus.eu/cdsapp#!/dataset/camsglobal-reanalysis-eac4?tab=form , last access : 1 st August, 2020) at a resolution of 0.75°x0.75° spatial and 3 hourly temporal 245 resolution. For this, the CAMS wind vector is spatially and temporally interpolated to the central coordinate of TROPOMI pixels.
To compute the NO2/CO ratio as function of the downwind distance x, TROPOMI and WRF data have been re-gridded at 0.1°x0.1°. 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 centre (see Fig 1a). The dimension of the box is motivated by multiple TROPOMI overpasses over Riyadh showing NO2 and 250 CO enhancements advected downwind over a ~200 km distance, without other large sources of NO2 and CO within a 100 km radius of the city centre (see Fig. 1a). Figure 1(b) shows the boundary layer averaged wind speed and wind direction over Riyadh indicating flow towards the northeast on 4 th of August, 2018. The box is rotated for every TROPOMI overpass depending upon the daily average wind direction within a 50 km radius from centre of Riyadh as shown in Figure 1(a) and Figure S10. The rotated box B1 is divided into N rectangular boxes, orthogonal to the wind direction with length (∆x) ~11 km 255 (see Fig. 1 and Fig. S10). The XNO2 and XCO grid cells that fall within the N rectangular boxes are selected to derive zonally averaged XNO2 and XCO for summer and winter.
Unlike the enhancements over the city, ∆XNO2 and ∆XCO become smaller than retrieval uncertainties at large distance from the city, where the ratio ∆XNO2/∆XCO becomes ill-defined. Therefore, we decided to use the ratio of mean XNO2 and XCO instead of enhancements over the background. To analyse the influence of atmospheric transport and the OH sink on the WRF 260 derived XNO2/XCO ratio two different ratios are derived: 1.
XNO 2 (emis,OH) XCO 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 analyse the influence of atmospheric transport on the ratio. For the comparison between TROPOMI and WRF, the CAMS backgrounds are included in "WRF RATIO" ( Table 1). The comparison of WRF RATIO to TROPOMI ratio, and the contribution of its components, is presented in section 3.2. 265

OH estimation: satellite data only
In the EMG method, following Beirle et al. (2011), 2D NO2 column densities 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 cross wind direction.
This way, average NO2 column density functions of the downwind distance to the city centre have been constructed for summer and winter (see Fig. S11). Using the EMG method as in Beirle et al., (2011), the e-folding distance x0 and NO2 emissions have 270 been estimated. The NO2 lifetime is derived by dividing x0 by the average wind speed (5.46 ms -1 and 5.24 ms -1 for winter and summer, respectively) and is provided in Table 2. The OH concentration is derived from the inferred NO2 lifetime using the IUPAC second order rate constant (for details see section Text S2 and S3). Rate constants at the time of TROPOMI overpasses are obtained from WRF by averaging the IUPAC second order rate constant from the surface to top of the planetary boundary layer. The PBL height at the time TROPOMI overpass has been taken from WRF. EMG derived NO2 emissions are also 275 converted to NOx emissions using the CAMS-derived conversion factor. Summer and winter averaged CAMS derived conversion factors for the box of 300 km x 100 km are 1.28 and 1.31, respectively.

OH estimation: WRF optimization
To jointly estimate the NOx 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). 280 The reaction of NO2 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 small perturbations (10 %) ∆background, ∆emission, ∆OH. The non-linear model is fitted to the observations, by optimizing scaling factors fBg, femis , fOH 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 factor femis, foh and fbg represent the modification of the prior in percentage change. 285 We estimate OH by optimizing WRF with TROPOMI in two ways 1) optimizing the simulated NO2/CO ratio using TROPOMI-derived ratios, named as "Ratio optimization" and 2) optimizing NO2 and CO separately using TROPOMI derived XCO and XNO2 named as "Component wise optimization". First the ratio optimization is described followed by the component wise optimization. Optimized ratios are derived as follows: is the change in F due to an increase in the XNO2 background by 5 % and a decrease in the CO background by 5 %. XNO 2 (emis,OH) is the contribution of city NOx emissions to XNO2 accounting for the OH sink, XNO 2 Bg 300 is the NO2 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 is the WRF derived XNO2 and XCO respectively. XNO 2 (emis,OH * 1.1) is the contribution of city NOx emissions to XNO2 after increasing CAMS OH by 10 %. The scaling factors femis ,fOH and fBg obtained from the ratio optimization have been multiplied by 10 % (i.e. divided by 10) as they represent changes in emission ratio, OH and Bg by 10 %. The scaling factors femis ,fOH and fBg obtained from the ratio optimization have been divided by 10 because 305 ∆F ∆emis , ∆F ∆OH and ∆F ∆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 NO2, it is not sensitive to the absolute emissions of CO and NO2. Therefore, we performed component-wise optimizations for XCO and XNO2 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 XNO2 is derived using Eq. (11). XCO is optimized using the same equation but 310 without considering the OH sink (see Appendix B).

13
Here, XNO 2 TROPOMI is the TROPOMI derived XNO2, XNO 2 WRF is the WRF XNO2. ∆XNO 2 emis is the change in XNO2 due to an increase in emission by 10 %, ∆XNO 2 OH is change in XNO2 due to an increase in CAMS OH by 10 % and ∆XNO 2 Bg is a change in the background XNO2 by 10 %. The scaling factors femis ,fOH and fBg are divided by a factor 10, because ∆XNO 2 emis , ∆XNO 2 OH and ∆XNO 2 Bg are defined as 10 % changes in NOx emission, OH and background level. 320

XNO2 and XCO over Riyadh
In this subsection, we compare WRF-derived XCO WRF and XNO 2 WRF with TROPOMI for summer (see  14 The comparison for summer in Figure 2 shows TROPOMI NO2 after replacing the TM5-based tropospheric AMF with WRF profiles as described in Visser et al. (2019). The enhancement of XNO2 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 investigate the use of the NO2/CO ratio to quantify OH in urban plumes. Due to the longer life-time of CO, the TROPOMI-observed XCO plume 330 extends further in the southeast direction compared to XNO2. Figure 2 shows that our WRF simulations are able to reproduce the TROPOMI retrieved XNO2 (r 2 = 0.96) and XCO (r 2 =0.78) plumes, confirming that WRF-derived 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 centre. 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 centre, XCO WRF is ~5.7 % higher than TROPOMI XCO. In 335 EDGAR 2011, emission sources are located in the centre 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 2 WRF and XCO WRF close to centre of Riyadh compared to TROPOMI.
In winter, the wind direction is predominantly from the south easterly sector in WRF and TROPOMI (see Fig S12). The spatial distribution of XCO WRF (r 2 = 0.73) and XNO 2 WRF (r 2 = 0.88 ) matches quite well with TROPOMI. Therefore, the difference 340 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 2 WRF is higher by 40 % to 50 %. The difference could either point to uncertainties in the NO2/CO emission ratio, uncertainties in the NO2 lifetime, or inaccuracies in the background. By quantifying OH, we can evaluate these explanations (see section 3.3). XNO 2 WRF is higher by 20 % in winter than in summer. Contrary, TROPOMI NO2 is lower by ~30 % in winter (Fig S12.) compared to summer (Fig. 2). Again, to 345 disentangle the role of changing sources and sinks, we need an independent estimate of OH.

The XNO2/XCO ratio and OH
Before comparing TROPOMI and WRF-derived XNO2/XCO ratios, we first analyse the influence of atmospheric transport and the OH sink on the WRF derived XNO2/XCO ratio. To do this three ratios are used 1. Ratio without OH 2. Ratio with OH 3. 350 WRF RATIO (see Table 1). As seen in Figure 3, S13 and S14, WRF is able to reproduce the TROPOMI-observed downwind evolution of XNO2 and XCO in summer and winter. The peak of the XNO2 and XCO plumes is shifted away from the city centre 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 355 influences NO2 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 NO2. This comparison demonstrates the 16 sensitivity of the relation between XNO2/XCO ratio and downwind distance to the NO2 lifetime, which we want to exploit to quantify OH. When including the background, the shapes of the functions in Figure 3c change (not shown), because the relative weights of the background and city contributions to the ratio vary with distance of the city centre. In summer, the WRF 360 RATIO is higher by ~15 % close to centre 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.
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 XNO2 by 40 % to 50 %. The WRF RATIO close to the centre of city is also 20 % 365 higher in winter than in summer, due to higher winter XNO 2 WRF than in summer (see Fig 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.

WRF optimization using synthetic data
To translate the discrepancies between TROPOMI and WRF derived ratios of section 3.2 into implied differences in emissions, OH and background, the least squares optimization method has been used as described in section 2.6. 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 375 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 femis , fOH and fBg are less than 2.5 % (see Fig. S16). This confirms that the least square 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 femis , fOH and fBg are expected to be higher due to 380 atmospheric transport errors, simplified chemistry, and XCO and XNO2 retrieval uncertainties . These errors did not play a role in the pseudo-data experiments, in which perfect transport and sampling was assumed.
To obtain a more realistic estimate of the uncertainty in least squares optimizated OH, TROPOMI data have been replaced by NO2 , CO and NO2/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, NOx and CO emission have been used in combination with 385 boundary condition for NO, NO2, CO, ozone (O3) from CAMS to run WRF-chem for August 17 th , 2018 and November 18 th , 2018 representing a summer and winter day, respectively.
For August 17 th , 2018, the ratio and XNO2 optimization increase the CAMS based prior OH of 1.19x10 7 moleculescm -3 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 x 100 km amounts to 1.33x10 7 moleculescm -3 , which is <5 % lower than the optimized 390 OH value that is derived using our method. The optimized NOx 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 x 10 7 moleculescm -3 by 19.4 %. The OH derived from WRF with full online chemistry is 1.07x10 7 moleculescm -3 and lower by 15.2 % than the optimized OH value. The component wise optimization increases the EDGAR NOx and CO emissions by 23.1 % and 10.5 %, respectively (see Fig S18). Overall, the uncertainty in optimized NOx, CO emission and OH derived from this test is <11 % 395 in summer and 10 % to 23 % in winter. Since the lifetime of NOx is determined by other reactions in addition to the oxidation to HNO3 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 uncertainty for OH, NOx emission and CO emission are in good agreement with the CLASS computations explained in detail in Text S6.

WRF optimization using seasonally averaged TROPOMI data 400
The results for summer are summarized in Figure 4, showing the optimized fit to the TROPOMI data as well as the corresponding scaling factors femis , fOH and fBg that are estimated. The optimized emission, OH and Bg obtained from the 2 nd iteration is divided by prior to derive the femis , fOH and fBg (see Text S5 for details). The convergence of the iterative procedure is shown in Fig S19 and S20 . The estimated uncertainties for the scaling factors femis, fOH and fBg are derived by summing the contribution of wind speed, length and width of the box, NO2 bias, CO bias and the different pathways of NOx loss in 405 quadrature (see Text S6, Tables S1 and S2). For summer and winter, the uncertainties of the optimized OH concentrations is <17 % and < 29 % respectively. For NOx 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 Χ 2 = 0.1 (for the derivation of Χ 2 see section Text S7) .The prior and optimized  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 NO2 and CO, but only their ratio.

415
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 4c). In the XNO2 optimization, the EDGAR NOx emission is increased by 42.1 ± 8.4 % and the CAMs background is reduced by 75.9 ± 10.0 %. CAMS OH is increased  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). 420 The summer optimized NOx/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 NO2 and CO emissions whereas the component wise has the full flexibility to adjust CO and NO2 emission. The NO2/CO ratio over a city is the sum of the contributions of the 425 background and the city emission. The relative weight of the two is determined by the absolute background levels and absolute emissions of CO and NO2.Therfore, the emission ratio estimated by ratio optimization is sensitive to the XNO2Bg. However, the difference between the two estimates is larger than expected but does not affect the OH estimation. Lama et al., (2020) inferred an NO2/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 430 from component wise optimization is consistent to Lama et al., (2020) and 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 Χ 2 = 0.11 . For winter, the ratio optimization increases 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. 435 The background ratio is reduced by 66.8 ± 11 %. The XNO2 optimization reduces the EDGAR NOx 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 (NOx/CO) derived from component wise optimization is 0.36 which is lower by 4.0 times than optimized emission ratio obtained from ratio optimization (see Table S4) .  uncertainty for EMG and WRF derived NOx emission and OH concentration is the sum of the contribution of wind speed, length and width of box and NO2 bias correction provided in Table S1, S2 and S3.

WRF optimization using a single TROPOMI overpass
To demonstrate the application of our WRF optimization method to single TROPOMI overpasses, results are presented in this 445 subsection for August 18 th , 2018. This date was selected for clear sky conditions with most of the TROPOMI NO2 and CO pixels passing the data quality filter. During this day, the urban plume is transported in southwestern direction over Riyadh.
The optimized ratio, XNO2 and XCO for a single day fit well with TROPOMI ( Χ 2 = 0.1, 0.3 and 0.7) comparable to the summer averaged plumes indicating that the optimization method can be applied to single TROPOMI overpass. The ratio 450 optimization increases the emission ratio and CAMS OH respectively by 111 ± 18.4 % and 37.9 ± 6.2 % respectively, whereas the background is reduced by 51.5 ± 5.2% (see Fig S22 a). The XNO2 optimization increases the EDGAR NOx emission by 25.5 ± 5.1 % and CAMS OH by 32.3 ± 4.4 %, whereas the NOx background is reduced by 54.4 ± 7.0 % (see Fig S20 b). The CO optimization doubles the EDGAR CO emission and reduces the background by 6.1 ± 0.97 % (see Fig S20 c). The optimized NOx and CO emission for August 18 th is 8.9 ± 1.7 kg/s and 18.9 ± 4.0 kg/s respectively and differs by <25 % with the summer 455 optimized emission (see Table 2 and S5). The optimized OH derived from a single TROPOMI overpass is 1.73x10 7 ± 0.3 molecules cm -3 differs by < 5 % from the summer averaged OH i.e. 1.7 x 10 7 ± 0.3 moleculescm -3 confirming that the method yields realistic results for a single overpass.

WRF optimization Vs the EMG method
To investigate the consistency between our method and the EMG method, the derived NOx lifetimes, emissions and OH 460 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 NOx emission and OH concentration. Both methods result in higher NOx emissions and shorter lifetimes in summer; lower NOx 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 emission in summer than in winter (Lange et al., 2021). During the summer, EMG and the WRF optimization method both increase the NOx emission and OH 465 concentration compared with the prior. The size of the NOx emission and OH concentration increase, obtained using the WRF optimization method is higher than the EMG method by 10% to 29 %. However, the difference between the EMG method and the component optimization method are smaller compared to the uncertainty 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 dissimilarity between the EMG method and the prior reduces after 470 optimization. The NOx emission after optimization differs from the EMG method by 33 %. Optimized OH concentration and NOx lifetime differs by <10 % compared to 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.

475
In contrast to EMG method, the optimization method can be used for a single TROPOMI overpass (see Section 3.6) and does not require yearly averaged NO2 data, allowing analysis of day-by-day OH, NOx and CO emission (see Section 3.3).
Segregation and averaging of NO2 urban plume by wind sector is not required in the optimization method. The effect of transport cancels out in taking the NO2 /CO ratio and loss of NO2 is mostly governed by OH during the mid-day. In this study, NOx emission and OH concentration is estimated iteratively whereas the EMG method arrives at the solution in a single step. 480 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 other NOx removal pathways. 485

WRF optimized emissions and emission trends
It should be realized that the a priori EDGAR emissions and TROPOMI optimized estimates represent different years ( Figure S23). For summer mid-day NOx emissions, the EDGAR emissions increased by 16.7 % from 2012 to 2018, which is lower than our optimization results. For winter, mid-day NOx emissions increase in EDGAR by 15.2 % from 2012 to 2018, whereas the WRF optimization yields reductions by 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 S4). Borsdorff et al., (2018b)  shows that for summer and winter NOx emission increases by 26 % from 2012 to 2018, which is higher by a factor 1.7 than EDGAR. CAMS-GLOB based summer and winter CO emission increases by 20 % from 2012 to 2018 which differs by ~40 % compared to EDGAR. In general, the relative increase in CO and NOx emission from EDGAR and CAMS-GLOB is much smaller compared to the difference with our optimization method. 505

Discussion
The TROPOMI retrieved XNO2/XCO ratio is useful for estimating mid-day 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 NOx and VOC emission, meteorological conditions, etc. (Sillman et al., 1990), which vary substantially between cities that are 510 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 NO2 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 NO2 this method cannot be applied. 515 We realise that our method only considers the first order loss of NO2 by OH forming HNO3. In reality, the NO2 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 NO2 via the formation of alkyl and multifunctional nitrates (RONO2) are also important reactions influencing the lifetime of NO2 (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 520 chemistry that includes all the sources and losses of NO2 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 NOx 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 XNOx,Bg is lost by OH. In this case the optimized NOx emission and OH for 525 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 S24 S25 shows that power plants and manufacturing industries are the largest pollutant emitter over Riyadh (Beirle et al., 2019).In this study, NOx and CO anthropogenic emissions 530 are introduced at the surface, whereas the emission height of different sources is expected to vary in reality. The different emission heights for NOx 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) (Guevara et al., 2020) suggests that temporal emission 535 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, spatial variations that is included. Road transport, CO emission has the largest contribution ~75 % to the total emission over Riyadh, whereas NOx emission from road contributes by 24 % to the total NOx emission. Residential combustion has the smallest contribution of ~0.3 to 0.4 % to total NOx and CO emissions (see Fig S24 S25 ). In the future, 540 the application of accurate diurnal emission factors for road transport (see Fig S25S26) can further improve the accuracy of urban OH concentrations estimated using TROPOMI derived XNO2/CO ratios. In addition, the seasonality for NOx and CO emissions is different in Riyadh than in Europe, which should be accounted for in future studies also.

Conclusions
In this study, a new method is presented for estimating OH concentrations in urban plumes using TROPOMI observed 545 XNO2/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 (Nov, 2018to March, 2019.
WRF is well able to reproduce the spatial distribution of TROPOMI retrieved XNO2 and XCO plumes over Riyadh during the summer and winter seasons. However, before the optimization, WRF overestimates XNO2 by 15 % to 30 % in summer and 550 25 derived XNO2/XCO ratio is higher by 15 % to 30 % in summer and 40% to 60 % in winter compared to TROPOMI, explained mostly by differences in XNO2.
The differences between WRF and TROPOMI observations have been used to optimize emissions and the NO2 lifetime. To this end, scaling factors for the city emissions, OH and the background level have been optimized iteratively using a least 555 squares method. Ratio and component wise optimizations have been compared to test the overall consistency of the method.
In summer, the ratio and XNO2 optimization for XNO2 suggest that the OH prior from CAMS is underestimated by 32 ± 5.3 %. The OH estimates obtained from the ratio and NO2-only optimization differs by <10 %, demonstrating the robustness of the method. Summertime emissions of NOx 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 560 of winter data, NOx 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 NO2.
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 two methods. These results 565 confirm that urban emissions and OH concentrations can robustly be estimated from TROPOMI data. With our method, single TROPOMI overpasses can be used to estimate OH whereas EMG method requires averaging of NO2 urban NO2 plume by wind sector. The iterative approach allows to test the factors i.e. femis, foh and fbg obtained from optimization method, whereas EMG method does not allows such flexibility.
An important remaining uncertainty is the bias correction of the TROPOMI XNO2 retrieval. Following the recommended 570 procedure, the air mass factor AMF is recalculated by replacing the tropospheric AMF based on TM5, that is provided with the data, with WRF-chem. The TROPOMI XNO2 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 XNO2 bias correction, the uncertainty in scaling factor for OH can vary up to 20 % and NOx emission to 60 % over Riyadh.

Appendix A: AMF recalculation 575
The air mass factor (AMF) used in the retrieval of TROPOMI XNO2 has been re-calculated by replacing the tropospheric AMF, calculated from the NO2 column simulated by TM5, with its WRF-chem equivalent, as described by Lamsal et al. (2010) and Boersma et al. (2016) where, M trop,WRF and M trop,TM5 are the tropospheric air mass factors derived from WRF and TM5, respectively. A trop,l is 580 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 NO2 column density in model layer l, based on WRF. A trop in Eq. (16)  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 595 is the change in XCO due to emission and ∆XCO Bg is the change in the XCO background level.