Evaluation of nitrogen oxides sources and sinks and ozone production in Colombia and surrounding areas

In Colombia, industrialization and a shift towards intensified agriculture have led to increased emissions of air pollutants. However, the baseline state of air quality in Colombia is relatively unknown. In this study we aim to assess the baseline state of air quality in Colombia with a focus on the spatial and temporal variability in emissions and atmospheric burden of nitrogen oxides (NOx = NO + NO2) and evaluate surface NOx, ozone (O3) and carbon monoxide (CO) mixing ratios. We quantify the magnitude and spatial distribution of the four major NOx sources (lightning, anthropogenic activities, soil bio5 genic emissions and biomass burning), by integrating global NOx emission inventories into the mesoscale meteorology and atmospheric chemistry model WRF-Chem. The comparison with in situ measurements is bound to urban areas whereas the use of remote sensing data allows to also evaluate air quality in remote regions. WRF-Chem was set up for a domain centered over Colombia with a similar resolution as OMI observed NO2 vertical columns as well as the EDGAR anthropogenic emission inventory, both providing information on a ∼20 km resolution. However, this apparently poses a challenge regard10 ing comparison with these urban observations. Air mass factors were recalculated based on the vertical distribution of NO2 within WRF-Chem, with respect to the coarse (1◦x1◦) a priori profiles because WRF-Chem is expected to better resolve spatial contrasts in NO2 profiles. The main reason for recalculation is a more consistent satellite-model comparison but it also reduced the mean bias. WRF-Chem was, on average, able to provide good estimates for tropospheric NO2 columns with an averaged difference of 0.02·1015 molecules cm-2, which is <5% of the mean column. However, the simulated NO2 columns 15 are overestimated in regions with abundant modeled lightning emissions and underestimated in regions where biomass burning emissions dominate in the model. This result reflects the high contribution by lightning emissions (1258 Gg N yr-1), even after already significantly reducing the emissions, and the low contribution by biomass burning emissions (104 Gg N yr-1) to total NOx emissions within the WRF-Chem domain. WRF-Chem was unable to capture NOx and CO urban pollutant mixing ratios, both in timing and magnitude. Yet, WRF-Chem was able to simulate the urban diurnal cycle of O3 satisfactory but with a 20 systematic overestimation of 10 ppb due to the equally large underestimation of NO mixing ratios and, consequently, titration. This indicates that these city environments are in the NOx-saturated regime with frequent O3 titration. We also applied an online meteorology-chemistry single column model (SCM) to evaluate how enhanced emissions and different representation of advection and mixing conditions could explain an improved representation of the observed O3 and NOx diurnal cycles. The 1 https://doi.org/10.5194/acp-2019-781 Preprint. Discussion started: 4 November 2019 c © Author(s) 2019. CC BY 4.0 License.

and seasonal variability during the last 1½ decade. It serves as a base to assess scenarios of future air quality in Colombia, or 25 similar regions with contrasting emission regimes, complex terrain and a limited air quality monitoring network.

Introduction
Nitrogen oxides (NO x = NO + NO 2 ) are one of the main precursors of lower atmospheric ozone (O 3 ). Exposure to NO x has an adverse effect on human health on acute and long-term basis (Panella et al., 2000;Wolfe and Patz, 2002). Likewise,O 3 is toxic to humans (WHO, 2003) and can also reduce agricultural yields (Ashmore and Marshall, 1998). Therefore, accurate 30 monitoring and predictions of surface concentrations of these air pollutants are key. Especially in densely populated regions air pollution has been a major concern and is expected to even have larger impacts in the future due to the continuous urbanization and increasing emissions from for example traffic.
Anthropogenic NO x is produced in combustion processes and is an indicator of industrial activity and transportation as well as other anthropogenic activities like biomass burning and agricultural activities. Anthropogenic sources add up to ∼70% (∼50% 35 industrial activity/transportation, ∼20% biomass burning) of the total global annual NO x emissions (Lamarque et al., 2010).
In addition to anthropogenic sources, natural sources contribute to total nitrogen budgets. NO emissions from soils add up to ∼12-20% of the global NO x emissions on a yearly basis (Bradshaw et al., 2000;Ganzeveld et al., 2002a;Jaeglé et al., 2005;Vinken et al., 2014). Lightning emissions are estimated to attribute on average 10-18% to the global yearly NO x emissions (Pickering et al., 2016). In the tropics (35 • N -35 • S), anthropogenic activities (7.81 Tg N yr -1 ), biomass burning (8.28 Tg N 40 yr -1 ), soil emissions (5.44 Tg N yr -1 ) and lightning discharges (6.33 Tg N yr -1 ) all contribute an approximately equal fraction to the total NO x emission budget (Bond et al., 2002). A modeling study in the tropics must therefore provide accurate estimates of all these source categories.
In Colombia, where economy is thriving after a period of civil war (Vargas et al., 2015), further industrialization and intensified agriculture have already resulted in-and are expected to further increase-NO x emissions (Ganzeveld et al., 2010). Previously, often limited by the presence of clouds.
During the last decades, computational advances have increased the possibility to conduct more detailed meteorology and air quality studies (Bauer et al., 2015). The recognition of the effects of chemical composition of the atmosphere on meteorology 60 have stimulated the development of online coupled meteorology/chemistry models (Baklanov et al., 2014). Nowadays, these models can be run for a large range of temporal and spatial scales. Not only the models, but also global emission inventories have considerably improved in spatial resolution during the last decades (González et al., 2018). Even though they may not provide enough spatial detail and heterogeneity for local scale (< 1 km) studies, e.g. to compare with in situ observations, they have provided essential information regarding emissions for regional scale (∼20 km) studies (Saide et al., 2012;Ghude 65 et al., 2013). In this study, rather than using high resolution urban emission inventories (e.g. González et al., 2018), we will demonstrate the importance of boundary layer mixing and advection in the comparison of simulated and observed in situ measurements.
The primary objective of this study is to assess the current baseline state of air quality in Colombia, diagnosed with a focus on NO x , using global emission inventories in a regional atmospheric chemistry model resolving the atmospheric chemistry 70 and meteorology at a resolution comparable to that of the emission inventories as well as the remote sensing observations. Furthermore, we evaluate surface NO x , O 3 and CO mixing ratios in urban regions. We are aware that concerns about air quality in Colombia are generally not limited to smog photochemistry mainly involving O 3 -NO x -VOC chemistry. Actually high concentrations of particulate matter might pose the largest risk to public health in many Colombian urban areas (Kumar et al., 2016). However, in this study we focus on NO x as an insightful metric to assess the spatial and temporal patterns in 75 air quality in this region given its role in O 3 photochemistry as well as the availability of remote sensing observations to be integrated with a bottom-up model analysis. In this study we use the Weather Research and Forecasting model coupled with Chemistry (WRF-Chem) (Grell et al., 2005). The model outcomes will be compared to in situ measurements and satellite retrievals to address the performance of the model both at the surface and integrated over the troposphere. This evaluation of surface and total column -using a highly resolving coupled meteorology-air quality model including the identification of regions that we refer to in this research.
The simulation was set up for one domain with a spatial resolution of 27 km centered at 4.89 • N, 71.07 • W. The entire domain consists of 100 grid points in both the North-South and the East-West direction with 60 vertical levels -in a sigma coordinate system-up to 50 hPa. The simulation length is one month (also given technical constrains on conducting much longer integrations with WRF-Chem), with a spin-up time of 24h, covering the whole month of January 2014. Selection of this study 95 period is motivated by the fact that January is the dry season in Colombia where loss of remote sensing data due to presence of clouds is minimized (see Sect. 3.1.1). In addition, in Sect. 5 we show how the selected study period can be put into context of the baseline state of air quality in Colombia using the interannual and seasonal variability in emission sources inferred from the remote sensing observations.  (Gery et al., 1989;Zaveri and Peters, 1999) is used here because it has been successfully implemented and tested in similar 105 studies (Gupta and Mohan, 2015). Additional parametrization schemes used in this research are listed in Table 1.
The Global Fire Emissions Database version 4 (GFEDv4) dataset (Randerson et al., 2015) provides us with the biomass burning emissions. GFED is available on a spatial resolution of 0.25 • x0.25 • , approximately the same size as the WRF-Chem grid cells.
The lightning-NO x parametrization scheme (Price and Rind, 1992), embedded in WRF-Chem, is used to account for NO x emissions by lightning. For this study we used an IC:CG (intracloud:cloud-to-ground) ratio of 2:1 constant over the whole domain with a flashrate factor of 0.1. Per lightning flash (both for IC and CG strikes), it is assumed that 250 moles of NO are emitted (Miyazaki et al., 2014). It has to be noted that in an initial simulation, using standard WRF-Chem settings (flashrate factor = 1.0 & 500 moles of NO per strike), resulted in a significant overestimation of the lightning emissions (see Sect. 4.1)

130
( Bradshaw et al., 2000;Miyazaki et al., 2014;Murray, 2016). The settings we used resulted in a twentyfold decrease of lighting emissions compared to standard WRF-Chem settings.
3 Observations of atmospheric composition

Satellite retrievals
Observational data on the distribution of NO 2 is retrieved from the Ozone Monitoring Instrument (OMI) onboard the National

135
Aeronautics and Space Administration (NASA) Aura satellite (Levelt et al., 2006). OMI measures, among other pollutants, NO 2 column densities (Boersma et al., 2007) with daily, global coverage. The pixel size of 24x13 km 2 may be coarse for particular applications, such as assessing urban pollution, but is suitable to assess contrasts in regional-scale air quality with apparent contrasting emission regimes. In addition, the resolution of the OMI observations is also comparable to the resolution of the anthropogenic emission inventory.

140
In this research we use the Quality Assurance for Essential Climate Variables (QA4ECV) NO 2 data product (Boersma et al., 2018). The measured slant columns -the tilted path directly from sun through the atmosphere to surface back to satelliteare converted to vertical columns using Air Mass Factors (AMFs) [-] by where VCD and SCD are the Vertical Column Density and the Slant Column Density [molecules cm -2 ], respectively. The

145
AMFs define the relation between slant column and the vertical column above a pixel based on external information on e.g.
surface albedo, scattering, clouds and the vertical distribution of NO 2 (Boersma et al., 2011). The vertical distributions of NO 2 in the QA4ECV product, which are used to calculate the AMFs, are simulated by the TM5-MP global chemistry transport model at a resolution of 1 • x1 • (Williams et al., 2017).

150
We follow the data filtering recommendations by the QA4ECV consortium. Presence of clouds (cloud radiance fraction > 0.5) led to omission of 63% of OMI NO 2 data. Figure 2 shows the amount of OMI data per WRF-Chem grid cell after filtering the observations of January 2014. Especially above mountainous regions, where we also find the main urban areas of Bogotá and Medellín, there is a lack of available data due to the continuous presence of clouds. This limits the quality of the measurements and which increases the uncertainty in the averaged tropospheric NO 2 column (Boersma et al., 2018). On average 9 data points 155 Figure 2. Spatial distribution of the available OMI measurements in January 2014 after filtering has been applied.
per grid cell are available for this specific domain in January 2014, but with a large spatial heterogeneity. Some areas have >20 data points and other only two valid observations in this month.

AMF recalculation
The AMF depends on assumptions of the state of the atmosphere and surface (e.g. surface albedo, cloud fraction, vertical distribution of NO 2 ) at the specific moment and location of a satellite observation (Lorente et al., 2017). This vertical sensitivity 160 is described by an averaging kernel, which describes the relationship between the true column and the estimated, or retrieved column (Boersma et al., 2016). High-resolution models such as WRF-Chem are expected to better represent spatial gradients in NO 2 profiles compared to coarse-scale global models such as GEOS-Chem or TM5-MP. Consqeuently, we can expect WRF-Chem to better resolve strong enhancements in tropospheric NO 2 VCDs in densely populated areas. Using grid sizes comparable to the size of such large urban areas is a major advantage of this procedure (Krotkov et al., 2017). The application 165 of the averaging kernel is shown to reduce systematic representativeness errors for a satellite-model comparison (Boersma et al., 2016). We can recalculate the AMF based on the a priori concentration profile x a (from the TM5-MP model) and the concentration profile in the high-resolution model x m , in this study WRF-Chem (Boersma et al., 2016): where M (x a ) is the tropospheric AMF used in the retrieval, A l are the elements of the averaging kernel for each l th vertical 170 layer and M (x m ) is the recalculated AMF. In a next step, the new VCDs can be calculated by dividing the SCDs (retrieved by the satellite) with the recalculated AMFs (Eq. (1)). On average, we find a mean decrease in AMF of 0.05 with a standard deviation of 0.15. Regarding inferred changes in the VCD due to this recalculation of AMF, we find a mean increase in the VCD of 0.02·10 15 molecules cm -2 (∼3% of the average 175 VCD) and a standard deviation of 0.07·10 15 molecules cm -2 . Above cities (e.g. Caracas, Bogotá, Medellín), we find mostly a decrese in AMF (Fig. 3a), and, consequently an increase in VCD (Fig. 3b). This indicates that there is more NO 2 present near the surface in WRF-Chem compared to TM5. This is consistent with our expectation that WRF-Chem better captures the sub-1 • x1 • processes that are not resolved by TM5, such as the localized urban emissions. Furthermore, we find increases in VCD (0.5·10 15 molecules cm -2 ) above the Amazon region. Decreases in VCD ( Fig. 3b) are found mostly across the border 180 from Colombia to Venezuela, better known as the Orinoco region ( Fig. 1). This reflects a higher abundance of NO 2 higher up in the troposphere from lightning sources, which are potentially underestimated by the global TM5 model (Silvern et al., 2018), combined with less NO 2 near the surface. We also find two isolated hot spots of decreases in VCD in southern Venezuela which correlate well with topography within the WRF-Chem domain which is less well resolved in the coarser resolution of TM5.
Despite locally significant changes in VCDs, a domain average of 0.6·10 15 molecules cm -2 indicates that the difference in the 185 NO 2 a priori profiles of TM5 compared to those in WRF-Chem does not lead to domain-wide significant changes in VCDs.

Comparison of OMI and WRF-Chem
In WRF-Chem we calculate the tropospheric NO 2 column by integrating from the surface to the tropopause, determined to be approximately the 50 th model level (∼90 hPa, ∼17km). This level is determined based on the average temperature profile (from surface to 50 hPa) of the complete simulation. Furthermore, to assess the daily differences in total NO 2 columns from time). Grid points with none or only one measurement after filtering (see Fig. 2) will be completely discarded. In this way we aim to get a reliable comparison between WRF-Chem and OMI, which enables us to determine systematic biases in the regions dominated by different emission sources.

In situ data 195
To further evaluate the model, observational data from air quality monitoring stations in Colombia are used. These include 30 stations confined to four cities in Colombia: Bogotá, Bucaramanga, Cali and Medellín (see Fig. 1). Observational data consists of 1-hourly averaged CO, NO, NO 2 and O 3 concentrations. The complete availability and locations of the station within the WRF-Chem domain can be found in Table A1. The data and information on the measurements is publicly available at http://sisaire.ideam.gov.co/ideam-sisaire-web/ (last access: March 2020). In this paper we only show the results of Bogotá prominent sources of NO x across the Colombian-Venezuelan border (Fig. 4b), in the Orinoco region, in our model study. This region is dominated by savanna type grasslands which emit a relatively high amount of soil NO x but also have a high probability of catching fire. NO x emissions in these regions are up to one or two orders of magnitude smaller compared to anthropogenic emissions, but, on the other hand, cover a larger area.
Lightning NO x emissions is the most dominant emissions source over land in 63% of all grid cells, followed by anthropogenic-220 (22%), biomass burning-(9%) and biogenic (6%) emissions (Fig. 4b). Since we use four different emission inventories, all with their own estimates and uncertainties, the distinct contrasts in the spatial distribution of emission sources will be key to  determine spatially heterogeneous biases in satellite retrievals compared to WRF-Chem. From budget calculations integrating over the whole domain, and using these January emissions to infer a NO x emission budget expressed per year (see Table 2) we also find that lightning NO x is the largest source in absolute terms, followed by anthropogenic, biogenic and biomass burning 225 respectively.

WRF-Chem & OMI comparison
To assess whether WRF-Chem is able to reproduce filtered and recalculated NO 2 VCDs satisfactorily we check for the spatial and frequency distributions for both WRF-Chem and OMI (see Fig. 5). For WRF-Chem we find a wide range of column densities (Fig. 5a). We find very low VCDs (∼0.3·10 15 molecules cm -2 ) over the Caribbean sea and across the eastern border 230 of Colombia into Venezuela. High VCDs in WRF-Chem are simulated above the major Colombian cities and the northeastern part of the domain (∼5·10 15 molecules cm -2 ) while the highest VCDs are simulated above the city of Caracas with values up to 8·10 15 molecules cm -2 . Similar to WRF-Chem, we find the lowest VCDs over the Caribbean sea in OMI (Fig. 5b). Also, we find the highest VCDs above major cities -most pronounced for Caracas and Medellín-but the magnitude of the OMI observed VCD (∼2.5·10 15 molecules cm -2 ) is much smaller compared to WRF-Chem. In OMI we find low VCDs above the 235 Amazon rainforest.
The large WRF-Chem VCDs (∼5·10 15 molecules cm -2 ) we find in the northeastern part of the domain (Fig. 5c)  VCD above Caracas might be due to an overestimation of anthropogenic emissions but this is not supported by a systematic major overestimation above cities (for example, we find no overestimation above Panama City or Bucaramanga). However, the EDGAR emission inventories are based on the year of 2010, when Venezuelan economy was still at its maximum (Wang and Li, 2016). After 2010, Venezuelan economy and oil production have declined strongly (Wang and Li, 2016) and therefore also emissions of pollutants have likely also decreased substantially. Lastly, we find a systematic overestimation in the WRF-Chem 245 simulated VCD above the Amazon rainforest. Even though the model overestimation is small in absolute terms (∼0.5·10 15 molecules cm -2 ) it is quite substantial relative to the background mixing ratios. In this region, soil NO x release is small, anthropogenic activities are hardly present and there are no known sources of biomass burning during January 2014. Also the role of advection from outside the model domain by the prevailing easterly winds seems to be limited indicated by lower VCDs at the eastern most grid cells. Consequently, overestimation in the simulated VCDs is most likely due to the simulated major 250 influence of lightning NO x emissions in this region (Fig. 4b), even though they have already been significantly reduced relative to the standard settings (see Sect. 2.2). However, we have to take into account that the OMI retrievals used for this comparison are near the detection limit and reflect those conditions when cloud formation, and therefore lightning production, is less active resulting in very low VCDs. In contrast, the co-sampled WRF-Chem columns might reflect simulated cloud cover resulting in production of NO by lightning. Nonetheless, the question whether lightning production was actually present or that it could 255 not be picked up by OMI, being less sensitive to the presence of NO 2 below clouds, remains unanswered.
Remarkably, we find a region with systematic underestimations ranging from the center of Colombia to the northeastern border with Venezuela (the Orinoco region). In this region, there is no presence of major cities and lightning NO x emissions are small. The discrepancy we find might be due to missing agricultural-or biomass burning emissions. Localized enhancements observed in OMI (∼2.5·10 15 molecules cm -2 ) might also be caused by biomass burning emissions since enhanced soil NO x are 260 expected to result in a more homogeneous enhancement of VCDs over a larger area with smaller intensities. We find that this intensity of biomass burning is not picked up by the WRF-Chem simulation using the GFED biomass burning inventory. Figure 5d shows, for both WRF-Chem and OMI, the frequency distribution in the NO 2 VCD. We find that both model simulated and observed VCDs show similar distributions, peaking at approximately the same VCD. However, WRF-Chem shows more outliers especially regarding the simulation of high NO 2 VCDs. The 90% confidence interval of the WRF-Chem simulated 265 VCDs is (0.33,1.33)·10 15 molecules cm -2 while for OMI the 90% confidence interval is (0.32,1.06)·10 15 molecules cm -2 , with medians of 0.59·10 15 and 0.56·10 15 molecules cm -2 , respectively.
We find the median and mean of the absolute overestimation by WRF-Chem to be 0.02·10 15 and 0.09·10 15 molecules cm -2 respectively (Fig. 5e). The 90% confidence interval equals (-0.43,0.70)·10 15 molecules cm -2 . The distribution is approximately Gaussian with a standard deviation of 0.53·10 15 molecules cm -2 but somewhat left-skewed indicating an overestimation by 270 WRF-Chem. This confirms the finding that WRF-Chem is able to produce on average good estimations for vertical NO 2 columns above Colombia. However, over-and underestimations can be significant, e.g. larger than the ∼10% uncertainty in monthly averaged OMI VCD over polluted regions (Boersma et al., 2018), due to numerous factors in both the model setup and the characteristics of the retrievals.

Surface mixing ratios
275 Figure 6 shows the model simulated and observed temporal evolution in NO x and O 3 over the whole simulation period. WRF-Chem represents the lower limit (generally mid-day) of observed NO x mixing ratios, but is unable to simulate the observed maxima (generally morning rush-hour) up to 200 ppb for particular events. Regarding O 3 , WRF-Chem resembles the upper limit of observed mixing ratios (∼40 ppb) during daytime but is unable to reproduce the observed low (<5 ppb) nighttime mixing ratios. 280 We retrieve January averaged diurnal cycles in NO x , CO and O 3 surface mixing ratios (Fig. 7a-c) by removing the significant spread in observed surface mixing ratios and by averaging the day-to-day variation in both the model and observations. Regarding simulated NO x , the average nocturnal mixing ratios are 20 ppb, with some day-to-day variation (standard deviation = 10 ppb), and a minimum of 2 ppb during daytime but with less day-to-day variation (Fig. 7a). The observations reach peak mixing ratios around 7:00 local time where vehicular emissions during rush hour are mixed in a shallow boundary layer increasing NO x 285 mixing ratios to 85 ppb on average. After rush hour mixing ratios decrease due to decreasing emissions, increasing boundary layer height and decreasing NO x lifetime. It is interesting to note that there does not seem to be a clear signal of evening rush hour in the NO x measurements and simulation.
The averaged diurnal cycle of CO in WRF-Chem shows a similar pattern to that of NO x (Fig. 7b). WRF-Chem shows daytime mixing ratios of ∼150 ppb (well above rural background mixing ratios of 100 ppb) and ∼350 ppb during nighttime while 290 the surface measurements show a significantly larger variation. Averaged surface measurements during rush hour exceed CO mixing ratios of 1500 ppb. Some measurement stations even report mixing ratios above 3000 ppb. Even though nighttime emissions are mixed over a smaller boundary layer we find that WRF-Chem underestimates surface mixing ratios of CO by a factor of 4 during rush-hour and by a factor of 2 for nighttime conditions. These ratios are similar to the NO x ratios in Fig. 7a.
Since CO has a relatively long lifetime compared to that of NO x we argue that observed differences regarding simulated and   We find that for WRF-Chem most of the NO x is present as NO 2 with NO mixing ratios being very close to 0 ppb (not shown here  (Fig. 7c). We also find that in WRF-Chem, the formation of O 3 immediately starts at 6:00 local time (sunrise) while for the observations we find the lowest mixing ratios at 7:00 local time due to the extra NO titration caused by rush hour. Nonetheless, it seems that chemical production and destruction rates of O 3 , as well as other processes contributing to the overall magnitude and diurnal cycle in O 3 , e.g., entrainment and deposition, are well captured by WRF-Chem considering the similar shape and amplitude of the diurnal cycle.
To test the hypothesis that the model-data mismatch over Bogatá is caused by a too coarse model resolution, or a misrepresen-310 tation of emissions, we conducted additional experiments with a Single Column chemistry-meteorological Model (SCM), also being previously applied for an analysis of observations of the plume of pollution downwind of the city of Manaus (Brasil) (Kuhn et al., 2010). The SCM simulates online, similar to WRF-Chem, atmospheric chemistry processes, including anthropogenic and natural emissions, gas-phase chemistry, wet and dry deposition and turbulent and convective tracer transport as a function of meteorological and hydrological drivers, surface cover, and land use properties (Ganzeveld et al., 2002b(Ganzeveld et al., , 2008.

315
For these urban area simulations with the SCM we have modified the surface cover properties by prescribing surface roughness at 1 meter, assuming a reduced vegetation fraction of 0.6, using a city area albedo of 0.18 and nudged the SCM meteorology with wind speed, moisture, and temperature profiles from the WRF-Chem simulation. The SCM is also nudged with long-lived tracers such as O 3 , NO x and CO above the boundary layer using WRF-Chem mixing ratios. Finally, we also used the same emissions, including diurnal cycle, as in the WRF-Chem simulation.

320
Using these settings in the SCM results in simulated January average diurnal cycles in NO x and CO, quite different from WRF-Chem but, in better agreement with the observations in terms of 30-day average diurnal cycles, maximum early morning peak and daytime minimum mixing ratios of NO x and CO ( Fig. 7a-b). The skewed O 3 diurnal cycle (Fig. 7c) is also better reproduced compared to WRF-Chem although the overestimation of the maximum afternoon mixing ratios is larger. Figure   7d shows a comparison of the SCM, WRF-Chem simulated-as well as the ERA5 reanalysis boundary layer height for the 325 grid point resembling the location of Bogotá. The SCM is showing a substantially deeper daytime maximum boundary layer with more day-to-day variation compared to WRF-Chem and ERA5 reanalysis data. The SCM also simulates a relatively fast afternoon transition to suppressed nocturnal mixing conditions reflected by a nocturnal inversion layer which agrees well with the ERA5 boundary layer height beging shallower than that simulated by WRF-Chem. Interestingly, the SCM simulation results in a better representation of the observed diurnal cycle of urban pollutant mixing ratios, especially regarding the observed 330 early morning maximum CO and NO x and minimum O 3 concentrations, without requiring the hypothesized enhancement in emissions. This stresses that, besides application of higher-resolution emission inventories and model experiments, the diurnal cycle in boundary layer dynamics (and advection) should be critically evaluated in models such as WRF-Chem which, however, would then also require urban boundary layer structure measurements.

335
The integration of global emission inventories in a highly resolved coupled meteorology-air quality model (WRF-Chem), with roughly the same spatial scale, allowed us to assess the state of-and contribution by different sources to the air quality in Colombia and neighbouring countries, diagnosed with a focus on NO x . We identified four major sources of NO x in Colombia which were implemented in WRF-Chem partly through emission inventories (anthropogenic and biomass burning) and partly through emission models (soil NO and lightning). Using January emissions to infer a NO x emission budget expressed per 340 year we found that lightning NO x emissions are the main source for the domain applied in this study, with 1258 Gg N yr -1 .
These are followed by respectively anthropogenic (933 Gg N yr -1 ), soil biogenic (187 Gg N yr -1 ) and biomass burning (104 Gg N yr -1 ) emissions. Figure 8 shows the averaged VCDs over the regions dominated by one of the four emissions classes (Fig. 4b). Figure 8a shows The top-down validation approach, using satellite retrievals, is a valuable tool to evaluate air quality in remote regions (Bailey 360 et al., 2006;Webley et al., 2012) with a missing network of air quality monitoring in both urban and rural sites. The daily global coverage and retrievals of NO 2 by OMI (Levelt et al., 2006) were used to assess the quality of all emission inventories over the whole domain. However, 63% of the data is lost for this specific model setup mainly due to the continuous presence of clouds.
Thus, longer simulation times have to be considered in the tropics compared to mid-latitudes. The vertical distribution of NO x within a modeling environment is key to identify discrepancies for a top-down validation approach using satellite retrievals.

365
It has to be recognized that the satellite sensitivity is reduced towards the surface (Boersma et al., 2016), inducing enhanced differences between observed and modeled profiles. However, this can be overcome by replacing a priori TM5 profiles with those from the applied model (Boersma et al., 2016).
In contrast to the bottom-up validation approach, where WRF-Chem showed a significant underestimation of NO x compared to the in situ measurements, we found that WRF-Chem does not systematically underestimate urban VCDs. This suggests that biomass burning emissions. Biogenic emissions are expected to show a more homogeneous distribution over a larger area with less pronounced peak emissions.. Therefore, they are also not expected to explain VCDs over 2·10 15 molecules cm -2 we found 375 in OMI retrievals. This connects to the findings of Grajales and Baquero-Bernal (2014) who concluded that high VCDs in this region are most likely related to biomass burning, which is apparently underestimated by the emission inventory we applied in this study. Soil NO emissions might be underestimated due to the missing anthropogenic term (fertilizer and manure application) (Visser et al., 2019). However, enhanced soil NO emissions due to the use of fertilizer is estimated to only contribute ∼1.3% to the total soil NO flux in the global chemistry-climate model EMAC (∼2.8 • ) for this domain (Ganzeveld et al., 2010).

380
Furthermore, the simulated soil biogenic NO flux (187 Gg N yr -1 ) from MEGAN is in the range between the total soil NO flux (230 Gg N yr -1 ) and the NO flux at the top of the canopy (105 Gg N yr -1 ) estimated by EMAC.
We find a large area with overestimations of modeled VCDs in the region dominated by lightning NO x emissions. These findings are in contrast with Grajales and Baquero-Bernal (2014) who found in their study with the GEOS-Chem modeling system that in remote regions without biomass burning there is an underestimation of modeled VCDs. Our study indicates 385 that lightning NO x emissions are a major source of NO x which might explain the discrepancy in the study by Grajales and Baquero-Bernal (2014) in which this source was not considered. Also, the use of WRF-Chem, having a spatial resolution approximately the same size as the OMI observations, can be advantageous over coarser models such as GEOS-Chem used by Grajales and Baquero-Bernal (2014). Further attention is required not only regarding the lightning NO x parametrization scheme, but also the model representation of convection and clouds, in follow-up studies on atmospheric NO x over Colombia, 390 or other regions where lightning is a dominant source of NO x . This study does not aim to provide comprehensive estimates of any of the emission sources using OMI data. Rather, we show the potential use of satellite data in a region with a limited air quality monitoring network in determining the regional scale air quality and NO x source regions. The use of cloud covered OMI observations to get a more comprehensive estimate of lightning NO x emissions (Beirle et al., 2010;Pickering et al., 2016) would make a very interesting follow up study.

395
The air-quality monitoring network in Colombia is limited to four major cities. This implies that the validation is limited to urban areas where anthropogenic emissions are the dominant source of pollution. A comparison with in situ data showed that WRF-Chem systematically underestimates urban surface mixing ratios of NO x and CO. The surface observations showed a clear signal of morning rush-hour emissions with average observed NO x mixing ratios up to 90 ppb and single observations not rarely exceeding 150 ppb. Similar to González et al. (2018), who focused on O 3 dynamics in Manizales (medium sized grazing shift rapidly towards a more intensified production of food, biofuels and rubber (Lavelle et al., 2014). Especially oil palm, which is one of the world's most rapidly expanding crops (Fitzherbert et al., 2008), is becoming more and more dominant in the Orinoco region (Vargas et al., 2015). Also, urbanization in Colombia is continuously increasing (Samad et al., 2012).

415
Ongoing and anticipated future transformation of both rural and urban areas, in combination with expected increases in temperature and changes in the hydrological cycle, imply changes in emission budgets affecting air quality in the future. Gg N yr -1 ), soil biogenic (187 Gg N yr -1 ) and biomass burning emissions (104 Gg N yr -1 ) all contribute to the total nitrogen emission budget. Especially the spatial distribution, clearly identifying regions with different dominating NO x sources, shows the importance of providing good estimates of each individual NO x source.
The top-down validation approach, using OMI retrievals, indicated a mean bias of NO 2 Vertical Column Densities (VCDs) of 0.02·10 15 molecules cm -2 , which is <5% of the mean column, with a 90% confidence interval of (-0.43, 0.70)·10 15 molecules 430 cm -2 . The VCDs in the Amazon region are overestimated in WRF-Chem, even after an already strongly reduced production efficiency, with respect to the low cloud free VCDs in OMI which is operating near the detection limit. This is a region where lightning NO x emissions are the only significant source of NO 2 . Additionally, the comparison indicates that GFED biomass burning emissions are potentially underestimated for January 2014 since OMI showed some strong enhancements in NO 2 not being reproduced by WRF-Chem. The biomass burning emission inventory shows some presence of wildfires in that region but 435 the model only produces estimates of VCDs of ∼1·10 15 molecules cm -2 , compared to OMI VCDs up to 2·10 15 molecules cm -2 , in regions where it is known to have significant biomass burning sources. Air Mass Factors (AMFs) were recalculated based on the vertical distribution of NO 2 within WRF-Chem with respect to the coarse (1 • x1 • ) a priori profiles for a more consistent model satellite comparison. An analysis of the past one and a half decade of OMI NO 2 VCD data showed that the selected simulation period is representative for the baseline state of air quality in Colombia but that interannual and seasonal variability 440 has to be taken into account in interpreting the OMI data and WRF-Chem simulations. The interannual variability in NO 2 columns over the different source regions can be attributed to specific events such as ENSO whereas the seasonal variability shows a strong enhancement of NO 2 VCDs above biogenic and biomass burning regions at the end of the dry season.
The bottom-up validation approach using air quality monitoring stations in urban areas showed that WRF-Chem, at the relative coarse resolution, does not reproduce these observations given the role of large heterogeneity in the emissions and other 445 processes determining pollution levels. Application of the anthropogenic EDGAR emission inventory (0.1 • x0.1 • resolution) in WRF-Chem resulted in a simulated underestimation of NO x and CO mixing ratios with respect to the local urban surface measurements. However, WRF-Chem was able to simulate the diurnal amplitude in O 3 reasonably well for all urban locations.
It seems that the underestimation of ∼10 ppb O 3 both during day-and nighttime can be attributed to the underestimation of NO by ∼10 ppb. Additional sensitivity simulations were performed with a Single Column Model (SCM) showed especially a 450 shallower nocturnal inversion layer compared to that simulated in WRF-Chem. This resulted in a better representation of the observed diurnal cycle of urban pollutant mixing ratio without the hypothesized enhancement in emissions. This indicated that besides the use of local emissions inventories in highly resolved modeling systems, it is also essential to carefully assess the role of boundary layer dynamics, in particular the representation of nocturnal mixing conditions, to evaluate simulations of pollutant concentrations.

455
In this study we presented a concise method, integrating both in situ and remote sensing observations with a mesoscale modeling system, to arrive at a quantification of air quality in regions with a limited measurement network to cover the large spatial heterogeneity in air pollution source distribution. Results obtained in this study provide insight in the baseline state of air quality in Colombia and which is essential to apply the presented combined modeling and measurement approach also to assess how air quality will further change due to future industrialization and land use changes.
Janić, Z. I.: Nonsingular implementation of the Mellor-Yamada level 2.5 scheme in the NCEP Meso model, US Department of Commerce,