Evaluating Sensitivity to Different Options and Parameterizations of a Coupled Air Quality Modelling System over Bogotá , Colombia . Part I : WRF Model Configuration

Meteorological inputs are of great importance when implementing an air quality prediction system. In this contribution, the Weather Research and Forecast (WRF-ARW) model was used to compare the performance of the different cumulus, microphysics and Planet Boundary Layer parameterizations over Bogotá, Colombia. Surface observations were used for comparison and the evaluated meteorological variables include temperature, wind speed and direction and relative humidity. Differences between parameterizations were observed in meteorological variables and Betts-Miller-Janjic, Morrison 2-moment and BouLac schemes proved to be the best parameterizations for cumulus, microphysics and PBL, respectively. As a complement to this study, a WRF-Large Eddy Simulation was conducted in order to evaluate model results with finer horizontal resolution for air quality purposes.


Introduction
Air quality is one of the main issues that are concerned by current atmospheric research.Global air pollution has an impact on human health, climate change and on the physics and chemistry of the atmosphere [1].Air pollution has become one of the most important interests of the local authorities in Latin America and represents the greatest social and economic costs of environmental damage after water pollution and natural disasters in Colombia [2].Urban agglomerations as Bogotá are major sources of regional and global atmospheric pollution with the pertinent environmental impact [3].Bogotá (4.6˚N, 74.1˚W) is the 5 th most populated city with around 7.4 million inhabitants in Latin America and one of most polluted cities [4], emissions from traffic linked to the increasing numbers of vehicles contribute to this concern [5] [6].
Air quality modelling has become a useful tool for administrations since it provides them a method to deal with human resources, production, emergency proceedings or to improve existing air quality plans and test abatement strategies [7].There are several air pollution modelling studies in South America [8] [9] but none of them are developed in Colombia or nearby countries.There are a few works focused on Colombia [10] [11] which analyze sensitivity of a mesoscale meteorological model to couple with an emission model and with a photochemical model.Together, these three models compose an air quality modelling system [12].Accordingly, implementing an air quality system in a particular area starts with setting up the meteorological model (the final aim of this study) which provides inputs for emission and photochemical models.The main interest of this work is to evaluate how the Weather Research and Forecasting (WRF) mesoscale meteorological model responses to different parameterizations during high air pollution episodes, and more specifically during days of high ozone concentrations in Bogotá.
Mesoscale meteorological models allow us to study and simulate meteorological variables.These models have a wide range of physical options to set up.It is a fundamental factor when configuring a model the selection of the physical parameterizations that are used to simplify somehow unresolved processes applying diverse approximations, the determination of the suitable model setup is one of the challenges when establishing a mesoscale model in a new region.Apart from the existence of a large array of available options, the best combination for one region is not necessarily applicable to another [13].
In this paper, we focus our attention on the meteorological modelling system.Exploring its sensitivity to variation in its configuration options, it is an important model evaluation exercise [14].In terms of air quality applications, the simulated concentration depends on the accuracy of this meteorological model and the importance of meteorological inputs on air quality modelling has been clearly stated [15] [16] so this analysis allow us to reduce the total uncertainty associated to the air quality modelling system since meteorological outputs are inputs both in the emission and photochemical models.Few studies of WRF sensitivity to diverse parameterizations exist over tropical regions, and most of them are related to PBL parameterization schemes [17]- [19].ARW (Advanced Research WRF) core has been used to obtain meteorological fields.Meteorological outputs were evaluated by means of statistical techniques.Numerical deterministic evaluation has been realized to compare modelling results with measurements.
Description of the studied area is presented in Section 2.1, as well as simulation domains and selected episodes.A characterization of the model and the methodology to evaluate it is presented in Sections 2.2 and 2.3, respectively.Detailed analysis of the experiments developed is presented in Section 2.4 and results obtained are presented in Section 3. Finally, some conclusions are reported in Section 4.

Methodology
In the following sections we show a more detailed description of the studied area features, the simulation domains and periods analyzed as well as a more comprehensive explanation of the modelling approach and model evaluation.

Studied Area, Simulation Domains and Episodes Selected
Following the aim of implementing an air quality modelling system in Colombia, Bogotá was chosen to perform WRF model sensitivity.
Bogotá is the capital of Colombia, the fourth biggest country in South America.It is divided into 32 departments and one capital district (Figure 1), Bogotá, the capital of the department of Cundinamarca and also treated as a department itself.Bogotá ranks fourth in the list of national capitals ordered by altitude with 2625 m above sea level.It lies in a 40 km wide and 100 km long plateau placed in one of the three Andean mountain ranges which cross the country.Mountainous complex terrain borders the high plateau.Its longest river is Bogotá River which has shown high pollution levels in recent years.Bogotá registers average yearly rainfall of 1013 mm and average yearly temperatures of 15˚C.

Modelling Approach
The Advanced Research WRF (WRF-ARWv3.5.1) mesoscale model developed by the National Center for Atmospheric Research (NCAR), USA, was the model chosen to conduct the simulations.It is a universally used community mesoscale model and a state-of-the-art atmospheric modelling system that is applicable for both meteorological research and numerical weather prediction.Different physical options that WRF offers can be combined in many different ways.Further details and description on this model appear in [20].WRF has different parameterizations for microphysics, radiation (long and short wave), cumulus, surface layer, planetary boundary layer and land surface.The initial and boundary conditions for domain D01 were supplied by the National Centers for Environmental Prediction and National Center for Atmospheric Research (NCEP/NCAR) Climate Forecast System Reanalysis (v1), with 0.5˚ (~55 km × 55 km) of spatial resolution and 6h of temporal sampling.Numerical simulations are executed for 48 hours corresponding on every day selected, taking the first 24 hours as spin-up time to minimize the effects of initial conditions and in order to represent a complete diurnal cycle.This is a common practice in meteorological modelling for air quality applications [21].
Two-way nesting was used for the three external domains (D01, D02 and D03) and one-way nesting for D04 and D05.The vertical structure of the model includes 32 vertical layers covering the whole troposphere and a resolution decreasing slowly with height in order to allow low-level flow details to be captured.The first 20 levels are inside atmospheric boundary layer (below 1500 m), with the first level at approximately 16 meters, and the domain top is about 100 hPa.The higher resolution close to the surface is a common practice in air quality studies in order to better represent the physical-chemical processes within de Atmospheric Boundary Layer [9] [12] [22]- [24].A total of 224 simulations have been run during the project development-14 configurations × 16 simulations/configuration.Meteorological modelling system works operationally in a computing cluster owned by Meteosim S.L. with 25 nodes and more than 212 cores.

Datasets and Model Evaluation
The evaluation performed is focused on the innermost domains, D04 and D05, since the final aim of this study is to find the best model setup for high resolution simulations.Meteorological observations were provided by 10 air quality stations that belong to the Red de Monitoreo y Calidad del Aire de Bogotá (RMCAB).Figure 3 shows the location of these stations and There are several methodologies for model evaluation that all together complement themselves [25] [26].The approach of comparing measurements with model results through different statistics (statistical deterministic approach) has been applied.The evaluations include the speed and wind direction at 10 m and air temperature and relative humidity at 2 m.Temperature (K) is calculated from WRF T2 predictions, wind speed (m•s −1 ) and wind direction (˚) are computed from U10 and V10 and relative humidity (%) is obtained from Q2 (water mixing ratio at 2 m), T2 and PSFC (surface pressure) using Magnus formula and specific humidity definition.The statistics have been calculated from hourly data of the model and observations, obtaining a global statistical value for the total period.These statistics provide information on how uncertain a model is in regard to the observations [27] and according to them a benchmark is given following Emery and Tai [28] suggestions.Table 2 shows the statistics used for model evaluation: the Mean Bias (MB), the Mean Absolute Gross Error (MAGE), the Root-Mean-Square Error (RMSE) and the Index of Agreement (IOA) and its benchmarks.
The circular nature of wind direction makes that statistical parameters should be carefully considered.Then, for the wind direction evaluation: D represents the minimum difference between modelled values and observed ones and it is always between -180˚ and +180˚ and N is the total number of measurements for all the days considered.if

Extent of the Sensitivity Analysis: Experiments
Many different physics options in WRF are available for microphysics, radiation, surface layer, land surface, Planet Boundary Layer (PBL) and cumulus.Physics options (schemes) considered in our study are listed in Table 3.We focus our attention on the study of cumulus, microphysics and PBL schemes; and radiation and land surface schemes have been fixed for all configurations: Rapid Radiative Transfer Model (RRTM) as a longwave radiation scheme [31] and the Dudhia scheme as a shortwave radiation scheme [32].One only option was tested as land-surface model (LSM): Noah LSM [33].RRTM, Dudhia and Noah LSM schemes correspond to the default WRF physical options.A total of 14 experiments have been evaluated progressively, as Table 4 shows.Three of them by varying cumulus parameterizations, two experiments by varying microphysics and a total number of eight by varying PBL schemes.We have focus most part of the configurations on PBL parameterizations due to the relevance of these schemes on air quality modelling [16].Additionally, an experiment has been undertaken at a higher resolution to find out the effects on predictions when increasing horizontal resolution.
Cumulus parameterization is used to predict the collective effects of convective clouds at smaller scales as a function of larger-scale processes and conditions.First, three configurations, i.e.Default, C1 and C2, were analyzed to take out the best cumulus parameterization between Kain-Fritsch (KF) scheme [34] that has a deep and shallow convection sub-grid scheme, Betts-Miller-Janjic (BMJ) scheme [35] [36] that is the most popular for tropical systems and Grell-Freitas (GF) scheme that is a stochastic convective parameterization for air quality modelling [37].Once cumulus option was selected, experiments M1 and M2 were evaluated together with the previous "best cumulus case" and with different Microphysics options.Microphysics parameterizations resolve water vapour, cloud and precipitation processes and that is the reason why they play such a significant role on air pollution levels [15].The three microphysics schemes considered have been the WRF Single-Moment 3-class scheme (WSM3) [39], the Stony Book University (Y.Lin) scheme [39] and the Morrison double-moment scheme (Morrison 2-mom) described in [40].
Several authors have recently shown the impact of PBL parameterizations on air quality modelling applications.Some of these examples would be the [15] or [41] sensitivity analysis.Consequently, taking into consideration the future air quality applications of this contribution, more experiments were tested by varying PBL parameterizations.A total of nine PBL schemes are evaluated in this study.Once cumulus and microphysics options were selected, experiments P1, P2, P3, P4, P5, P6, P7 and P8 were tested together with the previous "best cumulus and microphysics case" and with different PBL options.The schemes to describe vertical sub-gridscale PBL fluxes due to eddy transport in the atmosphere are the Yonsei University (YU) PBL [42], the Mellor-Yamada-Janjic (MYJ) PBL [35], the Assymetric Convective Model (ACM2) PBL [43], the Quasi-Normal Scale Elimination (QNSE) PBL [44], the Mellor-Yamada Nakanishi and Niino Level 3 (MYNN3) PBL [45], the Grenier-Bretherton-McCaa (GBM) PBL [46] that is a TKE scheme new in the WRF version used for conduct these simulations, the Bougeault-Lacarrère (BouLac) PBL [47] that is a parameterization of orography-induced turbulence, the UW [48] and the Total Energy-Mass-Flux (TEMF) scheme [49].The surface layer schemes calculate friction velocities and exchange coefficients that enable the calculation of surface heat and moisture fluxes by the land-surface models and surface stress in the planetary boundary layer scheme.These coefficients are computed by the similarity theory (MM5 similarity) surface layer scheme (described in [20]) for YSU, ACM2, GBM, BouLac and UW PBL schemes; similarity theory (Eta) surface layer scheme [36] for the MYJ PBL scheme and QNSE, MYNN and TEMF surface layer schemes for QNSE, MYNN3 and TEMF PBL schemes, respectively.
As a result of the experiments evaluation and comparison, a model setup was chosen for prospective air quality applications in Bogotá.Additionally, we have included into the analysis, a modelling experiment with finer horizontal resolution (333 m) over Bogotá centre (D05).meteorological maximum horizontal resolution places a restriction on the maximum horizontal of coupled air quality modelling systems.In order to couple the different meteorological scales and to deal with the step from regional to local scale are a state-of-art topic in the atmospheric modelling science [50] and several approaches have been evaluated during the last years to solve this problem.Every approach uses different frameworks to characterize sub-grid features.WRF model includes several urban parameterizations as the Urban Canopy Model [51] or the Building Effect Parameterization [52].Both of them present a major disadvantage because they need the use of detailed urban database.Moreover, WRF includes the possibility to use WRF with a large-eddy-simulation (LES) module that replaces the use of a traditional planetary boundary layer scheme.Other approaches are based on the coupling between air quality models indicated for different meteorological scales [53]- [55], or on a detailed monitoring of air quality levels to analyze sub-grid variability [56].To complement this work, we have focus our attention in one of these approaches and a Large Eddy Simulation configuration has been run at a finer resolution (D05).

Results and Discussion
Results of the comparison of every configuration are presented below using the proposed statistics.They have been compared for each meteorological parameter; temperature, wind speed, wind direction and relative humidity, and the one that showed best results for the maximum meteorological parameters was selected as "best case".It is necessary to clarify that in the event of a "tie" or not conclusive differences, wind direction will carry the most sway when selecting "best case" due to the importance of this variable in air quality modelling.

Cumulus Schemes
The first schemes analyzed have been cumulus.Wind direction errors are not within the benchmark for any of the simulations ran.Terrain complexity has a considerable influence on wind direction errors and the values found are substantially above the MB and MAGE benchmarks.However, these values were found in similar studies [12] [15] [16] [29].For the rest of the parameters, all of them follow the recommendation value (except wind speed RMSE for C1 (2.17 m•s −1 ) and C2 (2.15 m•s −1 ) configurations).
The three schemes produced similar results for temperature, with all values within the benchmarks and slightly overpredicting it.The C2 configuration produced the lowest MB for temperature (0.07 K) while the lowest MAGE (1.67 K) and highest IOA (0.91) corresponded to Default configuration, even though no significant differences are observed between them, as can be seen in Table 5.As for wind speed, C1 and C2 produced similar MB and RMSE values, it is the Default configuration which minimized wind speed MB (0.16 m•s −1 ) and  5 and wind statistics for wind direction, the cumulus parameterization of C1 (BMJ cumulus scheme) configuration provides the optimum results.For this reason BMJ was selected for next simulations to come as cumulus scheme.Graphics in Figure 4 [left] reflect the mean daily temperature evolution (a), the mean daily wind speed evolu-

Microphysics Schemes
Once BMJ cumulus parameterization was selected, three configurations were compared with this setting and by varying microphysics schemes: previous C1 "cumulus best case" using WSM3 microphysics scheme, M1 configuration with SBU-YLin and M2 using Morrison 2-moment.
Results for the three configurations with different microphysics schemes tested are shown in Table 6.The C1 configuration produced the lowest MB for temperature (0.13 K) while the lowest MAGE (1.64 K) corresponded to M2 configuration, while no conclusive differences where found for IOA for this parameter.M1 minimized wind speed MB (0.08 m•s −1 ) and wind speed RMSE (1.84 m•s −1 ).If we focus on wind direction, it is also C1 which produced the lowest MB (−9.30˚) but not the lowest MAGE (66.32˚) which is given by M2 configuration.Likewise, even though no significant differences were found for MAGE for relative humidity, M2 presented the lowest MAGE (10.34%) and highest IOA (0.80) together with C1.Although microphysics parameterization is considered to be highly influential for precipitation outputs and therefore wet deposition predictions [57], results for relative humidity are quite similar in three configurations (Table 6 and Figure 4(f)).According to these results, the better overall description was given by the Morrisson 2-moment microphysics parameterization that belongs to M2 configuration.
Graphics in Figure 4 [right] reflect the mean daily temperature evolution (d), the mean daily wind speed evolution (e) and the mean daily relative humidity evolution (f) for C1, M1 and M2 configurations comparing with the same observed parameters.Graphic for temperature (Figure 4(d)) shows that microphysics does not affect temperature profile significantly because similar results are observed for C1, M1 and M2. Figure 4(e) shows that wind speed tends to be overestimated for all configurations.

PBL Schemes
The last evaluation of WRF physics options involves PBL parameterizations.Once Morrison-2moment microphysics parameterization was set for the next configurations as a result of the C1, M1 and M2 experiments, other nine configurations were compared with this microphysics scheme and by varying PBL parameterizations as Table 4 summarize.Results are shown in Table 7. P6 produced the lowest MB (0.02 K) for temperature while M2 did the same with MAGE (1.64 K).PBL is also influential for wind speed, a parameter lightly overpredicted under all the PBL configurations tested.P4 reduced the MB (0.07 ms −1 ) and RMSE (1.73 ms −1 ) for wind speed.It is quite clear that P6 is the scheme that showed the best MAGE results for wind direction (57.24˚) reducing by up to 12% the average MB for all configurations (65.38˚).P6 also minimized relative humidity MB (0.26%) and relative humidity MAGE (9.30%) and improved the results of relative humidity IOA (0.82) with values of the three metrics that did not show important differences with P5.P8 produced the worst results for all the metrics calculated for temperature, wind direction MAGE and both relative humidity MAGE and IOA with remarkable variation between configurations (up to 18.73˚ difference in terms of wind direction MAGE and 0.17 difference in terms of wind direction IOA if we compare both with P6).According to this, P6 proved to be the best configuration improving the results for wind direction and relative humidity.Graphics in Figure 5 [left] show the mean daily temperature evolution (a), the mean daily wind speed evolution (b) and the mean daily relative humidity evolution (c) for M2, P1, P2, P3, P4, P5, P6, P7 and P8 configurations comparing with the same observed parameters.P6 is the best configuration in forecasting maximum wind speed and P8 the worst one.Almost all configurations accurately predict temperature, with the exception of P8, and the same conclusion can be drawn for relative humidity, for which P8 continues to show the worst results with the TEMF Planet Boundary Layer scheme.
P2 (ACM2 PBL scheme), P4 (MYNN3 PBL scheme), P6 (BouLac PBL scheme) and P8 (TEMF scheme) configurations turned to be computationally more expensive than the others (about 30% -40%) and P7 (UW scheme) up to 120%.In the later case, this can be explained by a reduction of the time step from 60 s to 40 s due to computational errors.

LES (Bogotá-333 m Resolution)
A finer-grid LES covering a smaller horizontal domain (D05) is nested inside a coarser-grid covering a larger horizontal domain (D04).This last contribution aims to validate the model results by increasing the resolution so that a future coupling of the meteorological model and the photochemichal model would be interesting in terms of air quality applications.M2 configuration was selected to run LES simulation because cumulus and microphysics parameterizations were already evaluated obtaining the best results.M2-LES is compared with M2 for a smaller domain (D05) so that validation is consistent including the same stations.
Comparisons between M2-LES configuration and M2 configuration within the D05 are shown in Table 8 and Figure 5 [right].Even though M2 (D05) improved most metrics for all the meteorological parameters, there are not conclusive differences between them and this is an interesting outcome as it would allow us to apply WRF-LES approach to forecast air quality at an urban scale without deteriorating the quality of results.Figure 6 displays a wind flow comparison between M2 (D05) and M2-LES.This figure shows that by increasing reso-  lution with LES approach, we find similar wind direction patterns and lower wind speed values for the same area at a higher resolution.

Conclusions
A total of thirteen WRF sensitivity experiments were conducted over Bogotá by varying cumulus, microphysics and Planet Boundary layer schemes during high air pollution episodes of 2010 and aiming to find the optimal setup of the model over this region.This work has focused most part of the configurations on PBL parameterizations due to its relevance on air quality modelling.We evaluate the differences in meteorological parameters of temperature, wind and relative humidity compared with observations in the innermost domain following a statistical analysis and the results show that no significant differences were found for temperature and relative humidity predictions depending on microphysics and cumulus parameterizations and no configuration perfectly works for all the variables.Among all the configurations analyzed, the best for the maximum meteorological parameters and selected as "best case" for cumulus, microphysics and PBL, proved to be P6, which improves the results for wind direction MAGE (57.24˚) and relative humidity MB (0.26%), MAGE (9.30%) and IOA (0.82).P6 has Betts-Miller-Janjic as cumulus scheme, the popular cumulus parameterization for tropical systems, Morrison 2-moment as microphysics scheme and Bougeault-Lacarrère (BouLac) as PBL scheme, a parameter- The model replicated temperature observations with a global index of agreement of 0.90.Not so precisely wind direction was predicted, but uncertainty of the prediction associated to this variable plays an important role.
Finally, a WRF-Large Eddy Simulation was included into the analysis, a modelling experiment with finer ho-rizontal resolution (333 m) over Bogotá centre (D05).This experiment was compared with M2 configuration and meteorological evaluation found that although the latter improved most metrics for all the meteorological parameters, there were not conclusive differences between them.These findings will allow us to couple WRF-LES with the emission and photochemical models at a higher resolution as an area of work for the future.However, default WRF physiographic data sets (topography and land uses) were used for 333 m resolution simulations.This analysis may be extended in the future by including higher resolution data sets so that we can accurately evaluate model performance of the LES approach.To achieve conclusive results, both in WRF and WRF-LES simulations, it will be useful to extend this study to a large period.

Figure 1 .
Figure 1.Main administrative divisions and topography map of Colombia [Image generated with ArcGIS].In Figure 2, we show modelling domains used for simulations.The WRF model is built over a mother domain (D01) with 27 km spatial resolution, centered at 4.6˚N, 74.1˚W.It comprises Central America, northern South America and part of Brazil and Peru, Pacific and Atlantic Oceans and Caribbean Sea and it is intended to capture synoptic features and general circulation patterns.The first nested domain (D02), with a spatial resolution of 9 km, covers northwestern South America and part of the Caribbean Sea and Pacific Ocean.The third nested domain (D03), with a spatial resolution of 3 km, comprises the Cundinamarca department and the fourth nested domain (D04) covers Bogotá.A fifth domain was included to take further the sensitivity analysis of WRF model at a higher resolution (WRF-Large Eddy Simulation): it is the innermost domain (D05), with a 333 m resolution.Table 1 shows the main characteristics of the simulation domains.Simulations were conducted in 16 specific days of the year 2010 (1-2/01/2010; 5-6/01/2010; 13-14/02/2010; 27-28/02/2010; 21-22/08/2010; 11-12/09/2010; 1-2/04/2010, 11-12/12/2010).These days present ozone concentrations above 60 ppb as a maximum running average over eight hours according to air pollution records supplied by the Red de Monitoreo y Calidad del Aire de Bogotá (RMCAB).

Figure 2 .
Figure 2. Modelling domains for simulations.[Image on the right generated using Google Earth].

Figure 3 (Figure 3 .
Figure 3. Meteorological stations used to conduct the analysis.[Images generated using Google Earth (a) and ArcGIS (b)].

Figure 4 .
Figure 4. Daily evolution of the mean hourly temperature (a); wind speed (b) and relative humidity (c) for cumulus experiments [left] and hourly evolution of temperature (d); wind speed (e) and relative humidity (f) for microphysics experiments [right].

Figure 6 .
Figure 6.Daily evolution of the mean hourly temperature (a); wind speed (b) and relative humidity (c) for PBL experiments [left] and hourly evolution of temperature (d); wind speed (e) and relative humidity (f) for finer resolution experiments [right].zation of orography-induced turbulence.The model replicated temperature observations with a global index of agreement of 0.90.Not so precisely wind direction was predicted, but uncertainty of the prediction associated to this variable plays an important role.Finally, a WRF-Large Eddy Simulation was included into the analysis, a modelling experiment with finer ho-

Table 2 .
Statistics used for model evaluation.

Table 4 .
Experiments developed and physics parameterizations.[Default configuration for default schemes in WRF].

Table 5 .
Statistical evaluation.Results for configurations by varying cumulus schemes.Results within the benchmark are highlighted in bold, and the best for each statistic is shaded in gray.
wind speed RMSE (1.90 m•s −1 ).Nevertheless, it is C1 configuration which produced the lowest MB (−9.30˚) and MAGE (66.43˚) for wind direction, and the lowest MAGE (10.45%) and highest IOA (0.80) for relative humidity.According to the results shown in Table

Table 6 .
Statistical evaluation.Results for configurations by varying microphysics schemes.Results within the benchmark are highlighted in bold, and the best for each statistic is shaded in gray.

Table 7 .
Statistical evaluation.Results for configurations by varying PBL schemes.Results within the benchmark are highlighted in bold, and the best for each statistic is shaded in gray.

Table 8 .
Statistical evaluation for M2 configuration (M2-LES and M2 (D05) comparison).Results within the benchmark are highlighted in bold, and the best for each statistic is shaded in gray.