Inﬂuence of the Grid Resolutions on the Computer-Simulated Surface Air Pollution Concentrations in Bulgaria

: The present study aims to demonstrate the effects of horizontal grid resolution on the simulated pollution concentration ﬁelds over Bulgaria. The computer simulations are performed with a set of models used worldwide—the Weather Research and Forecasting Model (WRF)—the meteorological preprocessor, the Community Multiscale Air Quality Modeling System (CMAQ)— chemical transport model, Sparse Matrix Operator Kernel Emissions (SMOKE)—emission model. The large-scale (background) meteorological data used in the study were taken from the ‘NCEP Global Analysis Data’ with a horizontal resolution of 1 ◦ × 1 ◦ . Using the ‘nesting’ capabilities of the WRF and CMAQ models, a resolution of 9 km was achieved for the territory of Bulgaria by sequentially solving the task in several consecutive nested areas. Three cases are considered in this paper: Case 1 : The computer simulations result from the domain with a horizontal resolution (both of the emission source description and the grid) of 27 km.; Case 2 : The computer simulations result from the domain with a horizontal resolution (both of the emission source description and the grid) of 9 km.; Case 3 : A hybrid case with the computer simulations performed with a grid resolution of 9 km, but with emissions such as in the 27 km × 27 km domain. The simulations were performed, for all the three cases, for the period 2007–2014 year, thus creating an ensemble large and comprehensive enough to reﬂect the most typical atmospheric conditions with their typical recurrence. The numerical experiments showed the signiﬁcant impact of the grid resolution not only in the pollution concentration pattern but also in the demonstrated generalized characteristics. Averaged over a large territory (Bulgaria); however, the performances for cases one and two are quite similar. Bulgaria is a country with a complex topography and with several considerably large point sources. Thus, some of the conclusions made, though based on Bulgarian-speciﬁc experiments, may be of general interest. Author Contributions: Conceptualization, K.G. and G.G.; methodology, K.G., G.G. and P.M.; software, K.G., G.G. and P.M.; formal analysis, K.G., G.G. and P.M.; investigation, K.G.; resources, G.G.; data curation, G.G.; writing—original draft preparation, K.G., G.G. and P.M.; writing—review and editing, K.G. and G.G.; visualization, G.G. and P.M.; supervision, K.G. and G.G.; project administra-tion, K.G. and G.G.; funding acquisition, G.G. All authors have read and agreed to the published version of the manuscript.


Introduction
Over the last years, it became obvious that our understanding of pollution and exposure processes at the regional-local-urban scales could be improved by combining multiscale models and creating new dedicated numerical approaches and that the representation of different scale interactions of dynamic, pollutant emission sources and pollutant transformations would be a critical element in obtaining more realistic simulation results.
Meteorology is one of the main uncertainties of air quality modeling and prediction. Many studies have investigated the role of meteorology on air quality [1][2][3][4][5][6][7][8][9][10][11][12]. The relationship between meteorology and air pollution is a result of a complex interaction between the atmospheric circulation and physical and chemical processes of the air pollutants in both gas and aerosol form. The improvement of atmospheric composition prediction capability is, therefore, tied to progress in both fields and their coupling. Recently, two-way coupling has been widely recommended as a more appropriate approach.

Modeling Tools
WRF v.3.2.1-Weather Research and Forecasting Model, [47], the meteorological preprocessor at CMAQ. The Weather Research and Forecasting Model (WRF) is a nextgeneration meso-scale numerical weather forecasting system designed to serve both operational forecasting and atmospheric research needs. It is the evolutionary successor to the MM5 model. The creation and further development of the WRF are due to the joint efforts of several US institutions such as NCAR, NOAA, NCEP, and others. WRF is a non-hydrostatic model with hydrostatic pressure coordinates following the terrain. The discretization is Arakawa-C type; CMAQ v.4.6-chemical transport model with proven qualities, applied worldwide and in European practices, [48][49][50]; SMOKE-emission model. It should be noted that SMOKE is very strongly oriented toward the American methodology for determining emissions with its categorizations and databases. For this reason, its application for modeling the levels of pollutants in Europe and Bulgaria is relatively limited. Most often, European scientists use emission models created by themselves. Intensive work is underway in a number of research groups to adapt SMOKE to European conditions. So far, the use of this processor is partly mainly for estimating biogenic emissions, emissions from large point sources and merging the various emission files, and saving them in the necessary formats.
It is important to note that the choice of this system of models is due not only to their great popularity and highly rated simulation qualities but also to the fact that in recent years these models have been well validated [51] and largely used for computer simulations of atmospheric composition of the Balkan Peninsula and Bulgaria [52][53][54][55][56][57].

Meteorological Data
The large-scale (background) meteorological data used in the study were taken from the 'NCEP Global Analysis Data' with a horizontal resolution of 1 • × 1 • . Using the 'nesting' capabilities of the WRF and CMAQ models, a resolution of 9 km was achieved for the territory of Bulgaria by sequentially solving the task in several consecutive, nested areas ( Figure 1a).

Modeling Tools
WRF v.3.2.1-Weather Research and Forecasting Model, [47], the meteorological preprocessor at CMAQ. The Weather Research and Forecasting Model (WRF) is a next-generation meso-scale numerical weather forecasting system designed to serve both operational forecasting and atmospheric research needs. It is the evolutionary successor to the MM5 model. The creation and further development of the WRF are due to the joint efforts of several US institutions such as NCAR, NOAA, NCEP, and others. WRF is a non-hydrostatic model with hydrostatic pressure coordinates following the terrain. The discretization is Arakawa-C type; CMAQ v.4.6-chemical transport model with proven qualities, applied worldwide and in European practices, [48][49][50]; SMOKE-emission model. It should be noted that SMOKE is very strongly oriented toward the American methodology for determining emissions with its categorizations and databases. For this reason, its application for modeling the levels of pollutants in Europe and Bulgaria is relatively limited. Most often, European scientists use emission models created by themselves. Intensive work is underway in a number of research groups to adapt SMOKE to European conditions. So far, the use of this processor is partly mainly for estimating biogenic emissions, emissions from large point sources and merging the various emission files, and saving them in the necessary formats.
It is important to note that the choice of this system of models is due not only to their great popularity and highly rated simulation qualities but also to the fact that in recent years these models have been well validated [51] and largely used for computer simulations of atmospheric composition of the Balkan Peninsula and Bulgaria [52][53][54][55][56][57].

Meteorological Data
The large-scale (background) meteorological data used in the study were taken from the 'NCEP Global Analysis Data' with a horizontal resolution of 1° × 1°. Using the 'nesting' capabilities of the WRF and CMAQ models, a resolution of 9 km was achieved for the territory of Bulgaria by sequentially solving the task in several consecutive, nested areas ( Figure 1a). (a)

Emission Data and Emission Modeling
The detailed emission inventory made by TNO, the Netherlands [58] was used for the domains outside Bulgaria. The inventory of emissions is made on an annual basis. The pollutants were calculated in groups such as CH 4 , CO, NH 3 , NMVOC (non-methane VOC, VOC-volatile organic compounds), NO x , SO x , PM 10, and PM 2.5 . The Bulgarian emissions are taken from the Bulgarian national emission inventory.
CMAQ, such as other chemical transport models, requires its input with emissions to be in a format that reflects the evolution over time of all pollutants involved in the chemical mechanism used. When preparing the input file for CMAQ, a number of additional procedures must be performed: • First, all primary information must be interpolated into the corresponding selected network/networks (gridded); • Second, time profiles should be imposed to modify the annual values so as to take into account seasonal, weekly, and daily variations in the work of the sources. • Finally, emissions from the "families" of organic gases and, to a lesser extent, SO x , NO x, and PM 2.5 must be split or "converted" into a larger number of components, according to the emission input requirements of CMAQ, which in turn depend on the chosen chemical mechanism-a procedure called "speciation".
In doing so, each of the different types of sources: surface (AS), large point (LPS), and biogenic (BgS), should be treated in a specific way (emissions from transport are also a separate category, but due to the way they are inventoried in our country, they are combined with area sources). Obviously, emission models are preprocessors needed for chemical transformation and pollutant transport models. One such component is SMOKE. Unfortunately, as already noted, it is very much adapted to the conditions in the United States-emission inventories, administrative division, categorization, combustion processes, etc.
For the purposes of the present study, time variations in emissions were calculated on the basis of daily, weekly, and monthly profiles provided in [59,60]. These time profiles are country, pollutant, and SNAP (Selected Nomenclature for Air Pollution) specific.
The "speciation" procedure depends on the chemical mechanism used. CMAQ supports various chemical mechanisms. For the purposes of the present study, Carbon Bond v.4-CB4 is used [61].
A specific approach to this speciation has been developed. It is proposed to follow the technology developed by the US EPA Emission Factor and Inventory Group [62]. More details about the speciation procedure can be seen in [63].
The inputs needed to calculate emissions are gridded data for area sources (AS), large point sources (LPS), and land use data needed to model biogenic sources (BgS). The latter emits organic matter, CO and NO, and their values depend heavily on weather conditions, including sunshine.
The area source data are processed by the specially designed AEmis program, which performs speciation and time profiles for each network cell for each SNAP for the respective Julian dates. The obtained hourly values of all 22 pollutants (CH4, CO, NH 3 , 10 types VOC, NO x , SO x , PMC, and 5 types PM 2.5 ) are saved in a file in NetCDF format.
The LPS database contains data for only 4 SNAP sectors-1, 3, 4, and 8. This information, together with the MCIP output, is fed to the SMOKE LPS processor, which produces the respective emission file. For this purpose, the inventory of powerful point sources is transformed into the requirements of SMOKE IDA format. This is a text file, but the order of the variables and their positions are fixed. This file, along with the inventory data, includes a number of parameters of the sources, such as geographical coordinates, height and diameter of the chimney, speed and temperature of discharge of pollutants, and more. The model not only performs speciation and time allocation but also calculates the so-called. "Plume-rise"-ejection of pollutants in height as a result of mechanical impulse and Archimedean forces. This increase in height also depends significantly on weather conditions-wind and atmospheric stability. As a result, SMOKE produces a 3D file-the pollutants are dumped at different levels (the levels coincide with the vertical structure of CMAQ).
SMOKE is also used to create a file with the third type of emissions-biogenic emissions. SMOKE currently supports the BEIS (Biogenic Emissions Inventory System) mechanism, versions 2 and 3 [64,65]. BEIS2 and BEIS3 are fed by the spatial distribution of the underlying surface type for the first step of the process-the calculation of normalized emissions for each network cell and for each underlying surface category (these are the emissions at fixed standard meteorological parameters). The final step is to bring normalized emissions up to date on the basis of grid and hourly meteorological information. The current version of SMOKE incorporates the BEIS3.13 mechanism [66].
Three cases are considered in this paper: Case 1: The computer simulations result from D2-the domain with a horizontal resolution (both of the emission source description and the grid) of 27 km, projected into the D3 grid. This case will be further referred to as C1; Case 2: The computer simulations result from D3-the domain with a horizontal resolution (both the emission source description and the grid) of 9 km. This case will be further referred to as C2; Case 3: The results of the computer simulations performed in D3 with a grid resolution of 9 km, but with emissions such as in D2-the emissions from each of the 27 × 27 km grid points are divided into 9 and allocated to each of the 9 points from the D3 grid, contained in the respective D2 grid cell. Thus the total emission amount of each of the D2 grid points is preserved, only distributed equally to the neighboring D3 grid points. This case will be further referred to as C3.
The simulations were performed, for all the three cases, for the period 2007-2014 year, thus creating an ensemble large and comprehensive enough to reflect the most typical atmospheric conditions with their typical recurrence.
The emission sources were constructed on the basis of 2005 emission inventory. This was performed on purpose-at that time, the emissions from the thermal power plants (TPPs) in Bulgaria were not reduced yet, so the experiment can also follow the effects of very large elevated point sources.

Results
Maps of the surface annually-averaged concentrations of NO 2 for cases C1, C2, and C3 are shown in Figure 2. The effect of the grid resolution is clearly manifested-the pattern of C2 is more detailed and displays the main road network and some of the big cities (large ground area sources). In the C3 fields, the local maximum near the southern border of Bulgaria can be noticed during the whole day. This probably is due to the emissions from Maritsa TPPs. This maximum is not observed in C1 and C2 fields, which suggests that it is perhaps created by the combined effects of coarse source description and detailed atmospheric dynamics and chemistry. The joint effect of coarse source description and detailed atmospheric dynamics and chemistry is also manifested by the fact that C3 fields are more diffused compared to C1 and C2 ones. The simulations were performed, for all the three cases, for the period 2007-2014 year, thus creating an ensemble large and comprehensive enough to reflect the most typical atmospheric conditions with their typical recurrence.
The emission sources were constructed on the basis of 2005 emission inventory. This was performed on purpose-at that time, the emissions from the thermal power plants (TPPs) in Bulgaria were not reduced yet, so the experiment can also follow the effects of very large elevated point sources.

Results
Maps of the surface annually-averaged concentrations of NO2 for cases C1, C2, and C3 are shown in Figure 2. The effect of the grid resolution is clearly manifested-the pattern of C2 is more detailed and displays the main road network and some of the big cities (large ground area sources). In the C3 fields, the local maximum near the southern border of Bulgaria can be noticed during the whole day. This probably is due to the emissions from Maritsa TPPs. This maximum is not observed in C1 and C2 fields, which suggests that it is perhaps created by the combined effects of coarse source description and detailed atmospheric dynamics and chemistry. The joint effect of coarse source description and detailed atmospheric dynamics and chemistry is also manifested by the fact that C3 fields are more diffused compared to C1 and C2 ones. The grid resolution effects are very well manifested in the maps of the annually averaged relative differences of the surface concentrations for different cases (Figures 3-6). The maps for NO2 virtually display the same effects as in Figure 2. It can be seen that the difference in the resolution of the road traffic emissions is also displayed in C2-C3 relative differences, but not so prominently as in C2-C1 relative differences. Due to the fact that the atmospheric dynamics in cases C2 and C3 are with the same resolution, the effect of the difference in the road traffic emission resolution is not so diffused and is concentrated The grid resolution effects are very well manifested in the maps of the annually averaged relative differences of the surface concentrations for different cases (Figures 3-6). The maps for NO 2 virtually display the same effects as in Figure 2. It can be seen that the difference in the resolution of the road traffic emissions is also displayed in C2-C3 relative differences, but not so prominently as in C2-C1 relative differences. Due to the fact that the atmospheric dynamics in cases C2 and C3 are with the same resolution, the effect of the difference in the road traffic emission resolution is not so diffused and is concentrated near the main roads only.
The relative differences between cases C1 and C2 are big for larger areas, which suggests a cumulative effect of emission and dynamics resolution. At noon, the relative differences between C1 and C3 show large values above mountain regions. This clearly is an effect of differences in dynamic resolution.
The relative differences for SO 2 are shown in Figure 4. The most remarkable feature in C2-C1 cases is the small spots with positive, fairly large relative differences located close to large SO 2 point sources-a clear effect of the emission source resolution. Significant local maximums and minimums can be seen for C1-C3 and C2-C3 cases, which probably is a result of the combined effects of emission source and atmospheric dynamic and chemistry resolution. Unlike the NO 2 case, the larger relative differences for SO 2 are between C1 and C3.
The biggest positive relative differences for PM 2.5 ( Figure 5) are between cases C2 and C1, and the biggest negative differences are between C1-C3 and C2-C3. The relative difference maps pattern shows the big positive differences are related to the road network and some area sources, while the big negative differences are related mostly to elevated large point sources.
As should be expected, the relative differences for O 3 ( Figure 6) are much smaller and not explicitly related to pollution sources.
Diurnal and seasonal course averages over the ensemble surface concentrations of different pollutants for cases C1, C2, and C3 are shown in Figures 7-9 for the city of Sofia, Maritsa TPPs, mountain point Rozhen (see Figure 1b) and averaged over Bulgaria.          The relative differences between cases C1 and C2 are big for larger areas, which suggests a cumulative effect of emission and dynamics resolution. At noon, the relative differences between C1 and C3 show large values above mountain regions. This clearly is an effect of differences in dynamic resolution.  The relative differences for SO2 are shown in Figure 4. The most remarkable feature in C2-C1 cases is the small spots with positive, fairly large relative differences located close to large SO2 point sources-a clear effect of the emission source resolution. Significant local maximums and minimums can be seen for C1-C3 and C2-C3 cases, which probably is a result of the combined effects of emission source and atmospheric dynamic and chemistry resolution. Unlike the NO2 case, the larger relative differences for SO2 are between C1 and C3.
The biggest positive relative differences for PM2.5 ( Figure 5) are between cases C2 and C1, and the biggest negative differences are between C1-C3 and C2-C3. The relative difference maps pattern shows the big positive differences are related to the road network and some area sources, while the big negative differences are related mostly to elevated large point sources.
As should be expected, the relative differences for O3 ( Figure 6) are much smaller and not explicitly related to pollution sources.
Diurnal and seasonal course averages over the ensemble surface concentrations of different pollutants for cases C1, C2, and C3 are shown in Figures 7-9 for the city of Sofia, Maritsa TPPs, mountain point Rozhen (see Figure 1b) and averaged over Bulgaria.   On all the graphs for NO2 (Figure 7), the diurnal and seasonal course is well displayed. For Bulgaria, the curves for C1 and C2 overlap almost completely throughout the day, as well as for C3 during the day when the concentrations have minimums. The concentrations of NO2 at night for C3 for Bulgaria have maximum and are the highest compared to the other cases.
For Sofia, the C1 and C3 curves overlap completely, with the exception of the maximum concentrations in the afternoon in summer. This can be explained by the fact that  On all the graphs for NO2 (Figure 7), the diurnal and seasonal course is well displayed. For Bulgaria, the curves for C1 and C2 overlap almost completely throughout the day, as well as for C3 during the day when the concentrations have minimums. The concentrations of NO2 at night for C3 for Bulgaria have maximum and are the highest compared to the other cases.
For Sofia, the C1 and C3 curves overlap completely, with the exception of the maxi- On all the graphs for NO 2 (Figure 7), the diurnal and seasonal course is well displayed. For Bulgaria, the curves for C1 and C2 overlap almost completely throughout the day, as well as for C3 during the day when the concentrations have minimums. The concentrations of NO 2 at night for C3 for Bulgaria have maximum and are the highest compared to the other cases.
For Sofia, the C1 and C3 curves overlap completely, with the exception of the maximum concentrations in the afternoon in summer. This can be explained by the fact that Sofia is a large and rather homogeneous NO 2 area source, so the source description resolution is not of much importance.
For Maritsa, TPPs C1 and C2 from November to February have the same course and close values. In the remaining months, there is a discrepancy in the phases of the minima and maxima in the average concentrations between these two simulations. From July to September, there was an overlap of the minimal values at noon of the average concentrations for all three cases. The average concentrations of C3 are the highest and in contra phase with the concentrations of C1 during the night. For Rozhen, it can be seen that during the winter C1 and C3 curves overlap almost all the time, but mostly in the afternoon when the average concentrations are again maximal, and in the other months, the overlap is at noon when the average concentrations are minimal and very close to those of C2. During the warm months, the concentrations of C3 are the highest, while those of C2 are the lowest throughout the year.
The average surface concentrations of SO 2 ( Figure 8) for all points also have a welldefined diurnal and seasonal course. As for Bulgaria, it can be seen that the average concentrations of the cases C1 and C2 are in contra phase with those of C3. For Sofia, there is no discrepancy between the minimums and the maximums for the different cases. It is interesting to note the period from November to May for Rozhen and Maritsa TPPs. For both points, the concentrations of the three cases have almost identical diurnal courses, while from July to October for Rozhen, the cases of C2 and C3 are in phase with those of C1. For Maritsa, TPPs C1 and C2 are in the contra phase with C3. On average for the country, C1 provides the highest average concentrations and C3 the lowest. For Sofia, from October to March, the average surface concentrations of C2 are the highest, and during the rest of the time, there is a slight dominance of the C1 with a more pronounced maximum in the evening during the summer. For Maritsa TPPs, the concentrations of C1 have the largest maximum values varying between 2-50 µg/m 3 , followed by C2 with values up to 35 µg/m 3 , while the variation of concentrations of C3 is in the range of 7-17 µg/m 3 . For Rozhen, the concentrations of C1 are the highest, and those of C2 and C3 are almost equal, those of C2 slightly higher.
The four graphs for PM 2.5 ( Figure 9) show that the concentrations of C3 are the highest, and for all cases, the concentrations have a well-defined diurnal and seasonal course. For Bulgaria, the concentrations have a maximum in the cold half of the year and a minimum in the warm. Generally, for the country, C1 and C2 are almost identical, as the concentrations of C1 are slightly higher. For Sofia, the concentrations of C1 or C2 are also almost identical; only in November the concentrations of C2 are larger than C1 and in the contra phase. The average concentrations for Sofia show that the minimum for all three cases is in April and May, which could be due to the "washing" of the spring rains.
For Maritsa TPPs, the concentrations of C2 are slightly higher than C1 as for the three cases, the maximums are in the cold half of the year and the minimums in the warm half of the year, which is probably due to lower TPPs production and more intense turbulent transport in summer. For Rozhen, differences between C1 and C2 are observed only at noon (lower values of average concentrations), as the concentrations of C2 are slightly lower than those obtained from C1. For this point, the seasonal course of C1 and C2 is weak, yet lower values of winter and spring concentrations are slightly noticeable.    Concentrations calculated on average for the whole country show the highest correlations and the lowest standard deviations, which is obviously due to the averaging over a sufficiently wide region. The highest correlations (higher than 0.7) are observed between scenarios C2-C3 and C2-C1.

Discussion
The NO 2 concentrations at the individual points ( Figure 10) show different regression and correlation dependencies. For Sofia, cases C1 and C3 show a significant decrease in values compared to C2, which can be seen from the regressions between them. Case C2 shows the highest scattering. In Maritsa TPPs, opposite results are obtained, the C2-C1 and C2-C3 regressions are close enough, and C3 shows the strongest scattering. The results for Rozhen are qualitatively and quantitatively close to the average for Bulgaria, which is obviously due to its remoteness from emission sources. For SO 2 in Sofia, where there are no significant sources of this pollutant, there are practically no diurnal changes in the correlation, which has a low value of about 0.5. For Maritsa TPPs, large differences were observed between the behaviors of the concentrations in the three cases.
The regression C2-C3 is very small, and the correlation between them reaches zero at some hours. This corresponds to the abnormally low C3 means shown in Figure 8. The scatter ratio in the three cases corresponds to the ratios between the averages.
The diurnal course of the scattering at this point is with a maximum of around noon, unlike the other points for Bulgaria as an average. The diurnal course of the scattering for Rozhen is similar to that obtained for Sofia, but the standard deviations in Rozhen are significantly higher than those for Sofia and especially from the average for Bulgaria. The average concentrations for these sites shown in Figure 8 do not differ significantly.
The fine particulate matters PM 2.5 have similar regression and correlation dependences both on average for Bulgaria and at specific points. The C2-C1 regression indicates an increase, C2-C3 is close to equivalence, and C2-C3 indicates a decrease. The correlation between C3 and C1 is the lowest for Bulgaria and all the selected points. Case C3 shows the highest scatter and also has the highest means, according to Figure 9.
For all the demonstrated compounds, the averaged over Bulgaria concentrations show a good fit and high correlation between C1 and C2 cases. The standard deviations for Bulgaria are also small. For the three selected points, the mutual behavior of cases C1, C2, and C3 is very different, largely depending on the emission source configuration.

Conclusions
The numerical experiments showed the significant impact of the grid resolution not only in the pollution concentration pattern but also in the demonstrated generalized characteristics. Averaged over a large territory (Bulgaria), however, the performance for cases C1 and C2 are quite similar.
Case C3, which is a kind of hybrid between cases C1 and C2, behaves strangely-it is not close to either of them. For Rozhen-a point distant from pollution sources, it correlates fairly well with C2. This probably is because, for the concentrations in points far from big sources, the resolution in the source description is not of much importance.
The above-demonstrated examples show that both the grid and the source description resolution play their role in the atmospheric composition formation-in different ways and to different extents, depending on the location and the specific pollutant. The grid resolution influences the atmospheric dynamics, the accuracy of the numerical solutions, and the chemical transformation rates. The experiments and the analysis performed do not suggest the idea of how these three factors interact to jointly form the atmospheric composition. Maybe applying the "Integrated process analysis" procedure of the CMAQ model can provide some clues in this direction. This could be a task for future work.