Intercomparison of the CALMET/CALPUFF Modeling System for Selected Horizontal Grid Resolutions at a Local Scale: A Case Study of the MSWI Plant in Krakow, Poland

Featured Application

Yielded research results indicate that the problems of the computational grid resolution in the CALMET/CALPUFF modeling system cannot be neglected during air quality impact assessment at a local scale for regulatory purposes. Decrease in grid resolution results in the underestimation of the highest maximum concentrations, which may lead to inadequate administrative decisions.

Abstract

Increase in grid resolution in atmospheric non-steady-state dispersion models induces a more faithful reflection of the area surface, and thus contributes to more detailed and diversified calculation results but also significantly prolongs the calculation time. This paper presents the influence of horizontal grid resolution in the CALMET/CALPUFF modeling system on the results of air quality impact assessment in a local scale carried out for the Municipal Solid Waste Incineration (MSWI) Plant in Krakow using the maximum permissible emission of NO_x. Subject to comparative analysis were four grids of the following resolutions: 100, 250, 500 and 1000 m. A direct intercomparison of air concentrations was made for 676 discrete receptors with the use of statistical indicators. On the basis of the calculations and analyses, it has been stated that, depending on the regular grid spacing, some differences in calculated concentrations can occur affecting the results of the air quality impact assessment. The highest concentrations in all computational receptors present in the given case were obtained for 100 m grid spacing. When compared to a grid of 100 m, the relatively smallest discrepancies were obtained for a grid of 250 m, with an already significantly shortened calculation time.

1. Introduction

The intense development in the remote sensing measurement techniques has contributed to establishing many spatial databases containing information on the terrain elevation and land cover, among others. Many resources of this type are made available for non-commercial applications by various organizations, e.g., by the United States Geological Survey (USGS) and the European Environment Agency (EEA). Special attention should be paid to the Digital Elevation Model data—Shuttle Radar Topography Mission (SRTM) [1] and the Digital Terrain’s Cover data—Corine Land Cover raster data [2], which are characterized by high resolution (approximately 100 m). More precise models reflecting the topography and land cover are also available, however they cover only the territory of North America. The above-mentioned sets of spatial data constitute the necessary source of information for the advanced dispersion modeling system of air pollution. Furthermore, they enable calculations at a very high spatial resolution of 100 m. Increase in the number of computational receptors significantly prolongs the calculation time, especially in more advanced models. Thus, in the case of more complex simulations and for larger areas, there is a tendency to limit the number of receptors and thereby increase the grid spacing.

The CALMET/CALPUFF modeling system (Exponent, Inc., Dallas, TX, USA) is more and more commonly used around the world in evaluations of the impact on air quality, including determination of shares of particular emission sources in shaping the level of air pollutant concentrations, especially in areas characterized by highly diversified landscape. A significant factor in this type of evaluations can be the method of description of the topography and land cover, as well as the resolution of the computational grid. The impact of topography on the modeling results of atmospheric dispersion using this model system was already noted (e.g., see [3,4,5]). In turn, the impact of the computational grid resolution on the air pollution modeling has been widely discussed thus far only for the mesoscale models with an extended description of chemical transformations in the atmosphere [6,7,8,9,10,11,12,13,14,15].

The review of the literature indicates that, depending on the size of the calculation domain, a different resolution of grids is used in the CALMET/CALPUFF modeling system. In the case of the mesoscale calculations, a resolution of several or even a dozen or so kilometers is used [16,17,18,19]; for the needs of air quality assessment in urban areas, a grid with the resolution of ca. 1 km is used [20,21,22,23,24,25,26,27,28]; and an analysis of the impact of single emission sources on air quality is usually conducted using finer computational grids with a spacing within the range of 100–500 m [4,29,30,31,32,33,34]. Sirithian and Thepanondh [35] tested four grid resolutions (0.75, 1, 2 and 3 km) in the modeling domain of 100 km × 100 km for open burning sources. In this case, it was concluded that the grid resolution of 1 km was the most appropriate resolution with regard to both accuracy of predicted data and acceptable computational time for the model simulation.

Some experiments were also conducted, with the use of a marker, to assess the accuracy of the CALMET/CALPUFF model for high grid resolution (e.g., at 200 m and 400 m) [36,37,38]. The results of the accuracy assessment of the CALMET model in high-resolution grids indicate that this model correctly predicts the wind speed and direction at near field, but the 1-h averaging may be too long to reflect the actual conditions. The CALPUFF model is characterized by general underprediction of the results for a near field; in addition, it is characterized by good results for short and long range at a prolonged averaging time [37]. MacIntosh et al. [36] stated that the model returns good results in a near field and for a long averaging period in the case of calculations of deposition of the particulates. The assessment of many dispersion models, conducted by Rood [39], at distances of 8 and 16 km from the emission source, showed that the CALPUFF model is characterized by a smaller error and a better correlation of the results than the stationary Gaussian slug models in the near field range area (<50 km). However, it requires a considerable amount of work associated with the preparation of data, adjusting the wind field model, and use of highly efficient calculation units. The research conducted by Dresser and Huizer [39] confirms that the CALPUFF model reflects high values of concentrations in the air well and predicts the results obtained for 3-h and 24-h averaging periods better than the AERMOD model. Its accuracy for neutral and permanent states of atmospheric stability is excellent, while in the case of unstable classes, some overprediction of calculation results can be expected.

Because the CALMET/CALPUFF modeling system is characterized by a good correlation between results and observations when using high-resolution grids, it is advised to recognize the impact of the changes in grid spacing on the results of calculations, as well as identify the resolution recommended for the purposes of the air quality impact assessment especially in the case of more complex terrain and land cover. The use of grids with too low resolution may in the case of such area create large discrepancies when compared to the results of calculations obtained for a dense grid. However, this type of action significantly shortens the time spent on preparation of the simulations and on the calculations themselves.

This paper presents an assessment of the impact of the horizontal grid resolution in the CALMET/CALPUFF models on the results of atmospheric dispersion modeling in a local scale characterized by a terrain of medium complexity and highly varying land use (urban, industrial and partly rural areas). The calculations were made for point emission sources related to the Municipal Solid Waste Incineration (MSWI) Plant in Krakow. To simplify the description of the problem, the calculations were limited only to one group of gaseous pollutants—nitrogen oxides (NO_x)—assuming their emissions at a constant level of the daily average emission limit value established for waste incineration plants [40]. NO_x emissions from MSWI plants cause the most significant air pollution from all emitted pollutants [41,42]. More detailed results of air quality impact assessment carried out for MSWI Plant in Krakow using CALMET/CALPUFF model with grid resolution of 100 m are presented in the literature [43]. The results of preliminary studies on the effect of horizontal grid resolution on modeling results for this facility are summarized in the literature [44].

2. Materials and Methods

2.1. The Study Area and Geophysical Data

A study area was adopted with the dimensions of 26 km × 26 km, covering almost the entire city of Krakow, with the emitters of the MSWI Plant located in the middle of this area (Figure 1). Krakow is one of the most urbanized areas in Poland (760,000 inhabitants), and the adopted test area is characterized by a particular topography (location of the city in low terrain impeding the undisturbed air flow, relatively high hills surrounding the city from two sides) and a high diversity of land cover classes (Figure 2).

The CALMET/CALPUFF modeling system requires preparation of the so-called geophysical data, which include the following information: elevation and land use, aerodynamic roughness length, albedo, Bowen ratio, anthropogenic heat flux parameter, soil heat flux parameter and leaf area index. For the purpose of this analysis, four rectangular computational grids of geophysical parameters of different resolutions (100, 250, 500, and 1000 m) were prepared. The grids, reflecting the land height, were prepared using the TERREL pre-processor, while the grid concerning land coverage, using the ArcMAP application (Esri, Redlands, CA, USA), in accordance with the methodology described by Oleniacz and Rzeszutek [45]. The source of the input data constituted the Digital Evaluation Model SRTM3 obtained from the resources of the United Stated Geological Survey [1] and the raster data of Corine Land Cover available at the European Environment Agency (EEA) service [2]. The MAKEGEO pre-processor was used to integrate the two abovementioned sets of data and to assign geophysical factors in respect of the land use category.

Figure 2 and Figure 3 present four grids of land use and digital elevation model prepared for the CALMET model for the particular examined resolutions. The graphic visualization of these grids indicates that an increase in the calculation resolution did not result in a greater diversity of land cover classes, however, it significantly influenced the spatial diversity of the occurrence of classes such as bare lands, rivers and forests. The grids with higher resolutions also enable a more accurate representation of the land orography.

2.2. Meteorological Grid

The meteorological data reflecting the conditions prevailing at area surface were obtained from the surface stations located in Krakow and four neighboring cities: Katowice, Tarnów, Nowy Sącz, and Bielsko-Biała. The following meteorological parameters obtained from these stations were used, recorded with the resolution of one hour: wind speed, wind direction, ceiling height, opaque sky cover, air temperature, relative humidity, atmospheric pressure, and type (code) of precipitation were acquired from the resources of the National Oceanic and Atmospheric Administration (NOAA) Earth System Research Laboratory [46], and then accordingly processed and formatted by means of the SMERGE pre-processor.

Radiosonde surveys were also used, obtained from three aerological stations, located in Poprad (Slovakia), Wrocław and Legionowo (Poland). The measurement data of this type, depending on the time interval and location of the station, are characterized by a varying time resolution. The intervals between various measurements can amount to 6, 8 and 12 h. During each measurement, the following meteorological parameters are recorded: pressure, temperature, wind direction, and wind speed. The data from the aerological stations were obtained from the database of the NOAA National Climatic Data Center [47]. They were subsequently processed and formatted using the READ62 pre-processor. This program also extrapolated the missing observational data to a pressure level specified by the user.

On the basis of the prepared meteorological and geophysical data, the meteorological and microclimatic parameter grids were established, with the use of the CALMET diagnostic model. The calculations were made for the meteorological data for 2012. During the calculation process, the following additional factors affecting the wind direction and speed were taken into account: a kinematic effect of terrain, slope flows and blocking effects. After considering the aforementioned effects, the model additionally carries out the procedure of minimizing the discrepancies in three dimensions to eliminate, e.g., the opposite wind flow.

According to the previous assumptions in the analyzed calculation area, four meteorological parameter grids were created, with the resolutions of 100, 250, 500 and 1000 m. Eight vertical layers were defined at the following heights: 20, 40, 80, 160, 300, 600, 1000 and 1500 m. Because the aerological stations are located more than 200 km from the calculation area, it was decided to include the vertical function of the wind field extrapolation with the use of the similarity theory on the basis of the data from the surface stations. Additionally, the silence extrapolation function was turned off and for each vertical layer the following BIAS variable values were adopted: −1, −0.8, −0.6, −0.5, −0.4, −0.3, −0.2, and −0.1; in this case, negative BIAS variable values mean a reduction in weight of the data of the wind field, originating from the aerological station.

2.3. The Emission Source Parameters

The source of the air pollutant emissions constituted the newly built MSWI Plant in Krakow, located in the western part of the city. The plant will have two independent lines of municipal solid waste incineration based on mechanical grate technology, with the total processing capacity of 220,000 tons/year. Flue gases from both incineration lines, after appropriate cleaning, will be released to the atmosphere with two identical emitters with a height of 80 m and inner diameter of 2.6 m. The calculations of atmospheric dispersion of pollutants include emission of sample gas substances—NO and NO₂ (NO_x) expressed as NO₂ (full conversion of NO to NO₂ was assumed). However, the assessed impact of the grid resolution on the calculation results for the adopted calculation area can refer to all produced gaseous pollutants (with the same dispersion parameters). To limit the impact of variability of pollutant emissions and flue gas parameters on the results of the calculations, the emission level of NO₂ and the velocity of flue gases at the outlet of both emitters were set at the same permanent level during the period of a year. The emission of NO₂ was set as the amount of the average daily emission limit values for waste incineration in the EU countries [40]. It amounted to 13.81 kg/h for each incineration line. The outlet velocity of flue gases, amounting to 13.84 m/s, was set as the function of the emitter’s diameter and the average flue gas stream volume, calculated theoretically for the waste combustion process in the MSWI Plant in Krakow.

2.4. Atmospheric Dispersion Modeling

The calculations of the pollutant spread in the ambient air were made using the CALMET/CALPUFF modeling system. CALMET is a diagnostic meteorological model, which generates variables, in terms of time and location, a three-dimensional field representing the wind and temperature. Additionally, on a two-dimensional scale, it calculates the microclimatic parameters determining the description of air pollution dispersion. On the basis of spatial data, this model determines the initial wind field (including such factors as kinematic effect of terrain, slope flows and blocking effects) and carries out a three-dimensional procedure minimizing the discrepancies, and then, using the observational data, enables an objective analysis of the final wind field [48]. CALPUFF is a multi-layer, multi-species non-steady-state puff dispersion model that simulates the concentration of air pollution. It has the ability to calculate dispersion of contaminants for point, surface and volumetric emitters, taking account of the variation in time of the level of air pollutant emissions and other considered parameters, such as: temperature, outlet speed of flue gases, the initial dispersion rates, etc. This model has the possibility to calculate in the area of a few dozen meters to several hundred kilometers [49]. It enables taking topography into consideration, along with dry and wet deposition of dust and gas substances and chemical changes of pollutants in the atmosphere. The CALMET/CALPUFF modeling system is recommended by the U.S. Environmental Protection Agency (U.S. EPA) for the calculations of air pollution dispersion to great distances and in a complex area [50,51]. It is also possible to use it for near field calculation, when a diverse wind field is of critical importance for the expected concentrations [52].

The modeling of atmospheric dispersion of NO₂ was performed for four variants characterized previously, differing in the resolution of the computational grid. The completed numeric simulations used the method of Simple CALPUFF of Terrain Adjustment, prepared on the basis of the assumptions of the CTDM model (Complex Terrain Dispersion Model). This method enables determination of the highest values of pollution concentrations in the air, as compared to other available options. For reasons of simplification, the calculations did not take the wet scrubbing mechanisms and chemical transformations into account [53]. In each grid receptor, the calculations of NO₂ concentrations in the air were made at area surface each hour. On the basis of the results of these calculations, the following values, among others, were defined: maximum 1-h (C_{1_max}) and 24-h (C_{24_max}) values, and average annual values (C_a).

Table 1. A general description of the adopted calculation variants of the air pollutant dispersion.

Variant	Grid Resolution (m)	Number of Receptors in the Grid	The Total Computational Time for the Grid 1
V1	100	67,600	453 h 43 min
V2	250	10,816	28 h 56 min
V3	500	2704	3 h 48 min
V4	1000	676	17 min

¹ with regard to a PC with an Intel Core i7-4930 K processor (12 M Cache, up to 3.90 GHz) BOX.

compares the adopted computational variants (differing in the resolution of the grid) and the corresponding duration of calculations. To ensure comparability of the obtained calculations results from various calculation variants, non-gridded discrete receptors were additionally defined in grids, in the number of 676, regularly distributed in the adopted calculation area (26 km × 26 km), precisely coinciding with the grid with the resolution of 1000 m.

2.5. Statistical Indicators

The comparative analysis of the obtained calculation results of air pollutant concentrations was conducted on the basis of numeric values in the form of statistical indicators included in the application BOOT Statistical Model Evaluation Software Package Version 2.0. For the evaluation of the models of air pollution dispersion, the U.S. EPA recommends six basic indicators: Fractional Bias (FB), Geometric Mean bias (MG), Normalized Mean Square Error (NMSE), Geometric Variance (VG), correlation coefficient (R) and FAC2 (fraction of predictions within a factor of two of observations) [54,55,56]. These indicators were used in this study in the form of following dependencies:

$$\mathrm{FB}=\frac{\left(\overline{C_1}-\overline{C_2}\right)}{0.5·\left(\overline{C_1}+\overline{C_2}\right)},$$

(1)

$$\mathrm{MG}=\exp \left(\overline{\ln C_1}-\overline{\ln C_2}\right)$$

(2)

$$\mathrm{NMSE}=\frac{\overline{{\left(C_1-C_2\right)}^2}}{\overline{C_1}·\overline{C_2}}$$

(3)

$$\mathrm{VG}=\exp \left[\overline{{\left(\ln C_1-\ln C_2\right)}^2}\right]$$

(4)

$$\mathrm{R}=\frac{\overline{\left(C_1-\overline{C_1}\right)\left(C_2-\overline{C_2}\right)}}{{\sigma }_{C_1}·{\sigma }_{C_2}}$$

(5)

$$\mathrm{FAC}2=\mathrm{fraction}\mathrm{of}\mathrm{data}\mathrm{satisfy}0.5\le \frac{C_2}{C_1}\le 2$$

(6)

where C₁ denotes results for Variant V1 (100 m); C₂ denotes results for Variants V2, V3 or V4; “ˉ” (top line) denotes the average over the dataset; and σ_C denotes the standard deviation over the dataset.

To conduct a detailed analysis, the values of false-negative (FN) and false-positive (FP) components of a fractional bias (FB) and of a geometric mean bias (MG) were determined. False-negative and false-positive components of the FB and MG indicators are determined according to the basic equation with the following assumptions:

• False-negative component is set only for those pairs (C₁, C₂), which meet the condition C₁ > C₂.
• False-positive component is set only for those pairs (C₁, C₂), which meet the condition C₁ < C₂.

The above parameters were calculated using spatially paired discrete receptors.

3. Results and Discussion

Table 2. Summary of minimum, maximum, average and standard deviation values for the calculation results for different types of NO₂ air concentrations obtained in the examined computational grids (variants).

The Type of Concentration	Variant	Grid Resolution (m)	NO2 Concentration (µg/m3)
			Min.	Max.	Average	Standard Deviation
C1_max	V1	100	0.61	167.0	9.8	9.1
	V2	250	0.60	140.0	10.3	8.4
	V3	500	0.84	133.1	11.1	9.3
	V4	1000	0.80	60.4	10.4	7.0
C24_max	V1	100	0.046	11.33	1.04	0.81
	V2	250	0.046	9.08	1.10	0.79
	V3	500	0.066	8.62	1.29	0.93
	V4	1000	0.070	7.83	1.19	0.76
Ca	V1	100	0.0030	1.712	0.106	0.112
	V2	250	0.0035	1.643	0.110	0.111
	V3	500	0.0054	1.263	0.125	0.123
	V4	1000	0.0052	0.896	0.111	0.100

presents the minimum, maximum and average values, as well as standard deviations, for different types of NO₂ concentrations obtained in the adopted calculation area for all receptors present under different examined variants of grids of different resolutions.

Table 3. Differences in minimum, maximum and average values of the analyzed NO₂ air concentrations obtained in the examined computational grids between Variant V1 and Variants V2, V3 and V4 (μg/m³).

Type of Concentration	Differences between Minimum Values			Differences between Maximum Values			Differences between Average Values
	V1–V2	V1–V3	V1–V4	V1–V2	V1–V3	V1–V4	V1–V2	V1–V3	V1–V4
C1_max	0.01	−0.23	−0.19	27.0	33.9	106.6	−0.5	−1.3	−0.6
C24_max	0.000	−0.020	−0.024	2.25	2.71	3.50	−0.06	−0.25	−0.15
Ca	−0.0005	−0.0024	−0.0022	0.069	0.449	0.816	−0.004	−0.019	−0.005

summarizes the differences in minimum, maximum and average values of these concentrations obtained in the examined computational grids between Variant V1 and other variants (Variants V2, V3 and V4). According to the calculations, using grids of lower resolution in the whole calculation area results in a decrease in the highest values (local maxima) among the maximum 1 h and 24 h and average annual concentrations, as compared to a grid with the resolution of 100 m (V1). In the case of grids with the resolutions of 250 m (Variant V2), 500 m (Variant V3) and 1000 m (Variant V4), this reduction amounted to 16%, 20%, and 64%, respectively, in the case of the highest values of maximum 1-h concentrations; 20%, 24%, and 31%, respectively, in the case of maxima of 24-h concentrations; and 4%, 26%, and 38%, respectively, in the case of maxima of average annual concentrations. Therefore, an increase in the computational grid resolution results in an increase in the accuracy of the area surface representation and enables reflecting the extremely high values of concentration that cannot be obtained during calculations performed in grids of lower density.

On the other hand, the computational grid with the highest examined resolution (Variant V1) also enabled mapping more points with low concentration values, as the average values of different types of concentrations from the entire examined calculation area are the lowest for this grid. The average area values of the maximum 1-h and 24-h and average annual concentrations for a grid with the resolution of 250 m were higher by, respectively, 4.2%, 6.1%, and 3.2%; for a grid with the resolution of 500 m by, respectively, 12.8%, 23.9%, and 17.5%; and for a grid with the resolution of 1000 m by, respectively, 5.8%, 14.3%, and 4.2% (as compared to the results obtained for a grid with the resolution of 100 m). The maximum overestimated average values (but also the minimal ones), as compared to Variant V1, were therefore obtained for a grid with the resolution of 500 m (Variant V3). In addition, in the case of Variant V3, the greatest spread of the calculation results of the maximum 1-h concentrations occurred, as compared to the average value, while the lowest possible spread was observed in the case of Variant V4 (a grid with the resolution of 1000 m).

It should be emphasized that, in the case of grids with the resolution of 100 and 250 m (Variants V1 and V2), the highest values of the maximum 1-h concentrations occurred at the same time, under almost identical weather conditions (wind speed ca. 0.5–0.6 m/s, atmospheric stability class 3, air temperature 259 K, and height of the mixing layer ca. 134–135 m). The grid receptors, for which the aforementioned maximum was obtained, in the case of Variants V1 and V2, are located, respectively, at a distance of ca. 650, and 500 m northeast of the emission source. On the other hand, in the case of Variants V3 and V4, the highest values of maximum 1-h concentrations occurred at different times and with different weather conditions. These conditions for Variants V1 and V2 were, respectively, as follows: wind speed 0.70 and 0.68 m/s; atmospheric stability classes 3 and 2; temperature 269 and 280 K; and height of the mixing layer 149 and 134 m. They were observed at the receptors located, respectively, ca. 400 and 600 m to the west of the emission source, that is, in a completely different place than in the case of Variants V1 and V2. Therefore, the use of grids with the resolution of 100–250 m, enabled obtaining the maximum 1-h concentrations in a similar time and place, while the less dense computational grids (in the range of 500–1000 m) did not make it possible, because the maximum values in these grids occurred in other places and under different weather conditions.

The highest values of the maximum daily concentrations in the case of Variants V1 and V2 also had similar characteristics with regard to the location (ca. 500 m east of the emission source), however, these maximum values occurred at different times and under different weather conditions. On the other hand, in the case of a grid with the resolution of 500 m (Variant V3), the maximum was obtained at a distance of ca. 1300 m, and, for a grid with the resolution of 1000 m (Variant V4), at a distance of ca. 1000 m from the source of emission (towards the north-east).

Figure 4, Figure 5 and Figure 6 present spatial distributions of maximum 1-h, maximum 24-h and average annual concentrations of NO₂ in the air obtained in the calculation area. Horizontal resolution of the calculation grid has some influence on the location of maximum values and the ability to identify local maxima. For example, the highest average annual values of concentration in the case of grids with the resolution of 100, 250, 500, and 1000 m occurred at distances of, respectively, ca. 700, 500, 1300 and 1700 m northeast of the emission source (Figure 6). In the case of grids with the resolution of 100 and 250 m, these maximum values were obtained near the emission source, while in the case of grids which resolution is 500 and 1000 m—at a distance up to three times larger. Significant differences in spatial distribution occurred also in the case of maximum 1-h (Figure 5) and 24-h (Figure 6) concentrations. Therefore, the use of coarse grids in the assessment of the impact on the air quality carried out using the CALPUFF model not only results in lowering the highest concentration values in the analyzed average periods, but also indicates different locations of these values in the adopted calculation area.

A more detailed analysis of this problem is possible by comparing the air concentrations calculated in a given spot to the reference concentration, which was the result of the calculation obtained for Variant V1 (a grid with the resolution of 100 m).

Table 4. Comparison of statistical indicator values of the model evaluation for the examined cases with regard to particular kinds of concentrations obtained in 676 discrete receptors.

Indicator	C1_max			C24_max			Ca
	V1–V2	V1–V3	V1–V3	V1–V2	V1–V3	V1–V3	V1–V2	V1–V3	V1–V3
FB	−0.05	−0.14	−0.06	−0.06	−0.22	−0.13	−0.04	−0.17	−0.05
FBFN	0.08	0.08	0.11	0.04	0.020	0.05	0.02	0.01	0.04
FBFP	0.12	0.22	0.17	0.10	0.24	0.18	0.06	0.18	0.09
MG	0.95	0.88	0.92	0.93	0.80	0.85	0.94	0.82	0.90
MGFN	1.08	1.08	1.10	1.04	1.02	1.04	1.01	1.00	1.02
MGFP	1.13	1.23	1.20	1.11	1.27	1.23	1.07	1.23	1.14
NMSE	0.10	0.32	0.24	0.04	0.16	0.13	0.04	0.12	0.06
VG	1.07	1.14	1.13	1.04	1.12	1.10	1.01	1.06	1.04
R	0.92	0.86	0.85	0.97	0.94	0.93	1.00	0.99	0.99
FAC2	0.99	0.95	0.95	1.00	0.96	0.95	1.00	1.00	1.00

compares the values of the particular statistical indicators obtained in the adopted calculation area for three cases reflecting a relative comparison of the results of concentration calculations obtained for the following pairs of variants: V1 and V2 (V1–V2), V1 and V3 (V1–V3), and V1 and V4 (V1–V4). Only the concentrations calculated in 676 discrete receptors were subject to this analysis, accurately coinciding with the points of the grid with the spacing of 1000 m.

The results of the calculations of the average fractional bias FB, presented in Table 4, clearly indicate that Variants V2, V3 and V4 are characterized by, on average, an overestimation of calculation results during all analyzed concentration averaging periods, as compared to Variant V1. The aforementioned overestimation does not exceed 30%, which means that the analyzed computational variants are in the range of the so-called good models, which in turn means that the obtained results did not differ significantly from the reference values (Variant V1). The fact that the analyzed cases show mostly average overestimation of the results in Variants V2, V3 and V4, as compared to Variant V1, is additionally proven by the fact that the values of the positive components of a systematic deviation FBFP are higher than the values of the negative components FBFN. It is worth emphasizing that the highest values of FBFP and the lowest values of FBFN, as well as the highest average deviation of the results in all examined concentration averaging periods were recorded for Case V1–V3, comparing Variants V1 and V3.

In the case of a geometric mean bias (MG), the good models include those which fit in the range of 0.7 ≤ MG ≤ 1.3 [54]. MG value = 1 means that the compared results are perfectly convergent. The MG values, presented in Table 4, show that each analyzed case deals with overestimation of the results, as compared to the V1 variant. The underestimation does not exceed 25% of the reference value. The analysis of the negative and positive MG components provides similar observations, as in the case of the values obtained for FBFP and FBFN. In this comparison, Case V1–V2 definitely returns the least diverse results, as it shows the lowest MG values for all average times, and a maximum underestimation of the calculation results at 7.5% of the reference value for the maximum daily concentrations.

The values of normalized mean square error (NMSE) and geometrical variance (VG) presented in Table 4, measuring both systematic and random errors, are within the so-called good models (NMSE <1.5; VG < 4.0), which in this case means good compliance. They show that, along with prolongation of the concentrations averaging time, the compared variants have smaller dispersion in relation to the average value. The largest deviation from the average value was obtained for Case V1–V3, while the lowest for Case V1–V2. However, for none of the analyzed cases, the twofold deviation of the results was reached in relation to the average value (NMSE > 0.5), therefore the deviations are very small. The above statements are reflected in the values obtained for the correlation coefficient (R) and FAC2. The values of FAC2 indicate that, in the case of average annual concentrations, two times the standard overestimations and underestimations of the calculation results, as compared to Variant V1, are very rarely observed. The highest number of double overestimations and underestimations of the results were recorded for Case V1–V4, with regard to the maximum 24-h concentrations. On the other hand, in Case V1–V2, only 1% of discrete receptors showed twofold overestimations or underestimations of the calculation results of the maximum 1-h concentrations, while there were no observed double overestimations or underestimations of the calculation results of the maximum 24-h and average annual concentrations.

Since the FB ratio is based on the average deviation, it is possible to obtain FB values = 0.0 for the compared variants (perfect convergence), even though, in reality, the modeling results are not convergent. Therefore, it has been decided to additionally present the FB values for each discrete calculation receptor to ensure better presentation of the differences between particular calculation variants. The fractional bias variability (FB), presented in Figure 7 for 676 discrete calculation receptors, shows a high diversity of the modeling results for different types of concentrations in the air obtained for grids of lower resolution, as compared to Variant V1 (the grid with the resolution of 100 m). An analysis of the examined cases indicates that Variants V2, V3 and V4 are characterized by a similar deviation, as compared to the reference variant. Furthermore, the standard deviations decrease along with the prolongation of the averaging period. It should be emphasized that the highest FB value ratios were recorded in each concerned averaging period for Case V1–V3, and the lowest—for Case V1–V2. This fact indicates that the mapping of Variant V4 is better than in the case of Variant V3, as compared to Variant V1.

The indicator value FB = 0 means a perfect convergence, and negative values mean an overestimation of the concentrations, as compared to the reference value (V1). The greatest spread of the calculation results occurred in the case of the maximum 1-h concentrations. In most receptors there is an overestimation of the modelled concentrations in the relation to Variant V1. For the examined Cases V1–V2, V1–V3 and V1–V4, respectively, 57%, 64% and 61% of discrete calculation receptors showed overestimation of the results, and due to the highest values of maximum 1-h concentrations, as much as a two times the standard overestimation of these results (FB < −2/3) was observed for, respectively, 1%, 4% and 3% of the analyzed receptors. It should be pointed out that the area of good compliance (|FB| < 0.3) includes, respectively, 75%, 68% and 58% of the aforementioned calculation results obtained in the discrete receptors for those cases. The lowest deviation of the obtained values was observed in Case V1–V2, the highest deviation of the results for FB > 0 was observed in Case V1–V4, and the highest deviation with respect to FB < 0—in Case V1–V3.

As it results from the fractional bias (FB) variability presented in Figure 7 for the calculations of maximum 24-h air concentrations, in every considered case for most of the analyzed calculation receptors, a value overestimation of these concentrations was obtained, in the case of Variants V2, V3 and V4, as compared to Variant V1. Depending on the case, 63–82% of the analyzed results were subject to overestimation. The highest and the lowest FB values were recorded for Cases V1–V3 and V1–V4, for which the values amounted to, respectively, −1.07 and 0.78. It should be pointed out that twofold underestimation of the results was recorded only in Case V1–V4, and twofold overestimation in Cases V1–V3 and V1–V4. The FB values, presented in Figure 7, indicate that for Variant V4 (Case V1–V4), for as much as 5.5% of the analyzed receptors, the twofold underestimation or overestimation was recorded as compared to the reference Variant V1. The range of the good compliance for Variant V2 (Case V1–V2, |FB| < 0.3) includes as much as 88% of the receptors subject to analysis.

The FB ratio variability for the calculation results of average annual concentrations presented in Figure 7, has similar characteristics to the variability of this ratio obtained for the maximum 24-h concentrations. The average annual concentrations received for Variants V2, V3 and V4 are characterized mainly by overestimation as compared to Variant V1. The most distinctive is Case V1–V3, in which 98% of the analyzed receptors showed overestimation of the calculation results. On the other hand, the range of |FB| < 0.3 in Case V1–V2, includes as much as 99% of calculation receptors. The results obtained for the average annual concentrations have smaller deviation with respect to the reference value, as compared to other averaging periods. It is also worth mentioning that twofold overestimation of the concentrations results was not recorded, and only one receptor in the case of V1–V4 showed at least a twofold underestimation.

Figure 8 presents the spatial differentiation of the correlation coefficients (R) for the analyzed cases (compared variants) with the discrete receptors considered and grouped depending on the distance from the emission source. The lowest values of the correlation coefficients were obtained between the calculation results of the maximum 1-h concentrations. Moreover, along with increase of the computational grid resolution, the correlation coefficient values decrease on general. Especially noticeable are the differences of the correlation coefficients between the groups of receptors located at different distance ranges from the emission source, which represent distinct terrain characteristics as well. This situation is also illustrated in Figure 4 and Figure 5, where depending on the computational grid resolution, the area depicting the occurrence of a certain concentration level changes. This stems from a fact that, in the case of low resolutions (e.g., 1000 m), a given area is represented solely by one terrain characteristic when actually, in its surroundings, there may be a greater diversification of the height of the terrain and land cover classes, which becomes visible only at higher grid resolutions. Consequently, the algorithms of the kinematic terrain effect, the flow of air masses from the slope or the blocking effects of the terrain in the wind field diagnostic model [48] result in obtaining different diffusion conditions for different horizontal resolutions.

On the other hand, according to the methodology described in Reference [49], establishing the actual receptor height resulting from the varying landform determines the air concentration level. Increase in the computational resolution causes the receptors located below or above the receptor representing a given area at lower grid resolution (greater grid spacing) to be included in the calculations. For example, the higher the elevation of the computational point, the shorter the vertical distances between the receptor and emitted plume, so the concentration obtained at the surface area is higher. A certain solution to this problem is the introduction of additional computational receptors, which represent the highest elevation points located within the impact of the emitter, as it is, for example, recommended for the AERMOD model [57]. It should also be noted that, at the long averaging periods, the resolution is not as significant, which is shown by high correlation coefficients between the results of calculations between Variant V1 and other analyzed variants for the annual average concentrations (Table 4 and Figure 8). In this case, the spatial distribution of the concentrations is determined primarily by the wind rose and not only by the terrain. However, it should be kept in mind that the wind rose characteristic for this terrain is dependent to a great extent on the flow direction of the main ventilation channel in Krakow, which runs along the Vistula River valley [58] and on the presence of elevations located mainly in the southern and northern part of the computational area (Figure 3).

4. Conclusions

As the results from the CALMET/CALPUFF modeling system assessment conducted by Cui et al. [37] using high-resolution grids (200 m) show, this system is characterized by some underestimation of the calculation results in the vicinity of the emitter. On the other hand, the research performed by Dresser and Huizer [39] for gaseous substances (SO₂) indicated that, with the grid resolution of 200 m, the CALPUFF model underestimates the maximum 1-h concentration values, at the same time allowing for obtaining very good results in a short range, for the averaging period of 24 h. On the basis of the calculations and studies presented in this paper, it was possible to obtain a partial answer to the question regarding to what extent the results of the gaseous pollutant dispersion in the CALMET/CALPUFF modeling system can vary for point sources with the height of 80 m, depending on the adopted computational grid resolution in the research area of 26 km × 26 km (the region of the Krakow Urban Area).

Increase in the grid spacing, from 100 to 250, 500 or 1000 m, resulted in a general (area average) overestimation of calculation results during all analyzed averaging periods, and thus their overestimation in the calculation area. The relatively smallest overestimation, as compared to a grid with the resolution of 100 m, occurred in the case of calculations made for a grid with the resolution of 250 m. For medium area values of maximum 1-h and 24-h and average annual concentrations, it amounted to, respectively, 4.2%, 6.1%, and 3.2%. On the other hand, in the case of a grid with the resolution of 500 m, the overestimation was the highest and reached, for these averaging periods, respectively, 12.8%, 23.9%, and 17.5%.

From the point of view of assessing the impact on the air quality, performed for point emission sources, obtaining the most reliable information about the level and location of the highest concentration values in the air requires the use of a computational grid with possibly the highest resolution. For this purpose, it is recommended to use a grid with the resolution of 100 m. Only in the case of a grid with such spacing is it possible to reflect the potentially least favorable situations, including coinciding of the local maximum with the calculation point. Unfortunately, along with the increase in the number of calculation points, a growth in calculation duration occurs (calculation costs), which, in the case of grids with high density and considering many point emitters in the analysis, may prove to be too great of a barrier. Application of a grid with the resolution of 100 m is characterized by a ca. 16-times longer calculation time, as compared to a grid with the resolution of 250 m, and as much as a 1600-times longer as compared to a grid with the resolution of 1000 m. Making these calculations for a possibly greater number of emitters will also contribute to an additional extension of the calculation time. The size of the calculation area is also significant here, as an increase in it leads to an increase in the number of receptors in a grid of a given density. High labor consumption, related to the preparation of the model input data, creates an additional barrier.

The use of a grid with the resolution of approximately 250 m seems to be good solution in the assessment of the impact on the air quality, as it is characterized by a slight general overestimation of the results as compared to a grid with the resolution of 100 m, and a relatively short duration of the calculations. Furthermore, the grid should reflect the close-range maxima results very well for long averaging periods. It is confirmed by the conducted statistical analysis of the modeling results of concentrations in the air obtained with various computational grid resolutions, especially in the analysis of average values of the examined statistical indicators and the variability of the fractional bias (FB) in the adopted discrete receptors. However, it should be pointed out that the calculations made for a grid with the resolution of 250 m do not guarantee achieving the highest values of maximum 1-h concentrations possible. To determine these maxima, it is therefore recommended to conduct additional calculations for a denser grid (e.g., of a 100 m spacing), which may be limited only to short episodes related to their occurrence. It is possible, since the highest values of maximum 1-h concentrations, in the case of grids with the resolution of 100 and 250 m, are likely to occur at the same time and at the same weather conditions.

The conducted calculations and analyses also enable taking a position on the problem of the use of computational grids with the resolution of 500 or 1000 m in the CALMET/CALPUFF modeling system for near-field applications. Considering the better proportion of obtainable overestimated values to the underestimated ones, the results obtained for the negative and positive FB components, as well as such statistical indicator values as NMSE, VG, FB, and MG, using a computational grid with the resolution of 1000 m is a slightly better solution than using one with the resolution of 500 m. However, an important argument, indicating the disadvantage of a grid with the resolution of 1000 m is the occurrence of high underestimations of the highest values of maximum 1-h and 24-h and average annual concentrations, respectively, at the level of ca. 64%, 31% and 48%, as compared to the results obtained for a grid with the resolution of 100 m. In addition, in the case of a grid with the resolution of 1000 m, the standard underestimations of the highest values of these concentrations occur several times. For comparison, in the case of a grid with the resolution of 500 m, the underestimation of the highest values of maximum 1-h and 24-h and average annual concentrations, as compared to the variant with the resolution of 100 m, amounted to ca. 20–26%.

The results yielded in the study correspond well with general guidelines regarding the selection of the grid resolution introduced in the manual describing the optimal settings of the analyzed modeling system [59]. According to these guidelines, it is necessary to find the optimum balance between the desire to make the grid size as large as feasible to reduce the run time and files size, and the need to make the grid size small enough to optimize the terrain effects on the wind field. Smaller domains for near-field applications may require a grid spacing of about 250 m, but if the dominant terrain features are not resolved it is recommended to change the grid resolution to about 150 m [59].

The present study provides an area for discussion with regard to the regulations concerning the use of the appropriate resolution for close range calculations with the use of the CALMET/CALPUFF modeling system. However, it does not exhaust the subject, and requires additional assessment and validation, based on the measurement data. It is recommended to take into account the problem of selecting the optimal horizontal grid resolution for a given case in further efficiency assessments of this modeling system at a local scale, especially in the case of terrain characterized by more complex topography and land use.