2.1 Study area

The Zonguldak province is located on the west coast of the Black Sea region in Turkey. The province is famous for coal mines and is called the land of “black diamond”. The study area, the Zonguldak metropolitan region, encompasses three counties (Zonguldak, Kozlu, and Kilimli) and three towns (Çatalağzı, Muslu, and Gelik). Çatalağzı is especially important because the coal-based thermal power plants are located there. Figure 1 shows the modeling domain, a 21 × 16 km area covering a total of 336 km². The region is surrounded by the Black Sea to the north and by hilly terrains to the south, east, and west. Meteorological parameters play an important role in the distribution of air pollutants in the area. Topographical and meteorological conditions prevent dispersion of air pollutants; therefore, inversion events occur frequently. The meteorological station in the modeling area is located at 41.4492 N and 31.7779 E. The Black Sea climate predominates, with an annual average temperature of 13.6 ºC. July and August are the warmest months, while January and February are the coolest. The annual average rainfall is 1217.8 mm. The yearly mean wind velocity is 2.3 m s^-1 and the predominant wind directions are ESE, SE, and SSE. The annual average pressure is 1000.1 hPa and the average relative humidity is 71.2% (^{Tecer et al., 2012}; ^{TSMS, 2019}). There was only one air quality monitoring station during the modeling period (2007-2011) and the station continuously monitored ambient concentrations of sulfur dioxide (SO₂) and PM (PM₁₀). The location of the Zonguldak AQMS is also shown in Figure 1 (41.4513 N and 31.7833 E). According to the Address Based Population System, 183 651 inhabitants lived in this region in 2011.

Fig. 1 Study area. 

2.2 Emission inventory

An emission inventory can be defined as a compilation of the amount of all air pollutant entering the atmosphere from all sources in a geographical area for a certain time (^{Elbir and Muezzinoglu, 2004}). The emission inventory of this study was prepared for PM₁₀ emissions of domestic heating, transportation, and industrial activities. Emissions of illegal fuel usage for heating were not included due to insufficient data. In this study, the PM₁₀ emission inventory was compiled for 2011. Since official emission factors for Turkey are quite limited, they were selected from two widely known databases published by the US-EPA and the European Environmental Agency (European Environmental Agency). PM₁₀ emissions for point and area sources were calculated by using Eq. (1).

where E is the amount of PM₁₀ emissions (t year^-1), Fueli is the activity data (total fuel consumption per year) for each fuel type, EF _i is the emission factor for each fuel type and ER is the emission reduction in percentage, providing there is an air pollution control equipment. Eq. (2) was used for line sources (motor vehicle emissions).

where E is the amount of PM₁₀ emissions (t year^-1), Veh _i is the number of vehicles per type, D _i is the distance traveled in a year for a specific type of vehicle, and EF _i,km is the emission factor per vehicle type (i) per km (^{Çetin et al., 2007}).

2.3 Air quality modeling

Air quality modeling is the mathematical simulation of the destination of air pollutants after release from source to receptor points including their movements (transport and dispersion), reactions, and removal (dry and wet deposition) processes (^{Holmes and Morawska, 2006}). The Gaussian dispersion model explains the dispersion of plume concentration in the vertical and horizontal axis perpendicular to the wind direction (^{Fraile et al., 2006}). On the other hand, puff models are used to model dispersion and movement of discrete emissions (series of puffs over time) from the stack. The Gaussian puff model successfully calculates pollutants concentrations in case of low wind speeds while the Gaussian dispersion model generally fails (^{Holmes and Morawska 2006}; ^{Gulia et al. 2015}). Dimensions of the puff are determined with dispersion parameters, which are functions of the puff’s travel duration (^{Korsakissok and Mallet, 2009}; ^{Abdul-Wahab et al., 2011}).

The AMS/EPA Regulatory Model (AERMOD) was developed by the American Meteorological Society (AMS) and the United States Environmental Protection Agency (US-EPA). AERMOD is one of the most widely used regulatory models. It is a near-field steady-state Gaussian plume model and can model point, area, and volume sources in flat and complex terrain (^{Holmes and Morawska, 2006}; ^{Yassin and al-Awadhi, 2011}). Line sources can be added to AERMOD as a series of area sources. This model can calculate pollutant concentrations up to 50 km from the source. The AERMOD model formulation was well defined by ^{Cimorelli et al. (2005)} and ^{Jeong (2011)}. To calculate pollutant concentrations, the Gaussian probability distribution function is used both in the vertical and horizontal axes in the stable boundary layer, and the horizontal axis in the convective boundary layer. However, vertical distribution in the convective boundary layer depends on the bi-Gaussian probability density function (^{Vijay Bhaskar et al., 2008}; ^{Yassin and al-Awadhi, 2011}). AERMOD uses two preprocessors: the meteorological preprocessor AERMET and the terrain preprocessor AERMAP (^{Jeong, 2011}; ^{Kakosimos et al., 2011}). Hourly meteorological data required by AERMET are cloud cover, temperature, relative humidity, pressure, wind direction, wind speed ceiling height, precipitation, and horizontal radiation. AERMET also needs some surface parameters like albedo, Bowen ratio, and surface roughness. Then, it uses the planetary boundary layer theory and calculates planetary boundary layer parameters like Monin-Obukhov length, friction velocity, sensible heat flux, and convective scaling velocity. With the help of these parameters, convective velocity scale and mixing height can be estimated (^{Cimorelli et al., 2005}; ^{Kakosimos et al., 2011}). The terrain preprocessor AERMAP is used for calculating terrain heights at receptor points (^{Heckel and Lemasters, 2011}; ^{Kakosimos et al., 2011}).

Another air quality model used in this study is The CALifornia PUFF (CALPUFF) model developed by Sigma Research Corporation. CALPUFF is a non-steady-state Lagrangian Gaussian puff model (^{Indumati et al., 2009}). Gaseous and particulate pollutants can be modeled by CALPUFF using space and time-varying meteorological conditions. CALPUFF can model point, area, and volume sources. Like AERMOD, line sources can be added as a series of area sources in this model. The CALPUFF formulation calculates partial penetration, buoyant and momentum plume rise, stack and building effects, dry and wet deposition (^{Holmes and Morawska, 2006}; ^{Abdul-Wahab et al., 2011}). CALPUFF has three components: CALMET preprocessor, CALPUFF, and CALPOST postprocessor. CALMET calculates hourly meteorological data and wind fields, while CALPOST is used to calculate averaging and reporting concentrations depending on CALPUFF results (^{Barna and Gimson, 2002}; ^{Oshan et al., 2006}). Details of the CALPUFF model formulation can be found in ^{Jeong (2011)} and ^{Scire et al. (2011)}.

In this study, AERMOD View and CALPUFF View software, developed by WebLakes, was used for modeling. PM₁₀ emissions of two thermal power plants, 11 main roads, 44 districts, and the harbor were entered as source inputs in both models. Area sources of residential heating were effective only in the heating season (from October 15 to May 15). Receptors were defined as a rectangular grid with 250 m between adjacent receptors. In total, 3576 receptor points were placed in the modeling domain. All sources and receptors were placed by using Universal Transverse Mercator (UTM) projection (Zone 36 North) on WGS-84 Datum. The same configuration of sources and receptors was applied for both models. Five years (2007-2011) of AERMET and CALMET ready meteorological data (hourly surface and upper air) were obtained from WebLakes. The reason for using 5 years of meteorological data was to increase the accuracy of forecasting. The wind rose plot of the given years is provided in Figure 2. Surface parameters of the modeling area, required by the AERMET preprocessor, are shown in Table V. Since the modeling region consisted of three different land-use types (sea, forest, and urban area), weighted averages of surface parameters were determined and used in the AERMET preprocessor. Although AERMOD View provides digital terrain data from United States Geological Survey (USGS) web server, global digital elevation model (GDEM) data obtained from the Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) sensor, onboard the Terra satellite was used due to GDEM’s better resolution (30 × 30 m). ASTER GDEM files can be freely obtained from the NASA Earthdata website (^{Earthdata, 2020}). In the CALPUFF model, a geophysical preprocessor is necessary to obtain the geophysical data file (geo.dat) required by CALMET. ASTER GDEM’s digital terrain data were used for topography. These files were converted into 7.5 min DEM format to run them with a geophysical processor. A 1 × 1 km grid resolution was applied t meteorological grids, producing 21 cells along the x-axis and 16 cells along y-axis. Ten vertical cells were applied, and default settings were used for cell face heights (0, 20, 40, 80, 160, 320, 640, 1200, 2000, and 3000 m). The Diagnostic Wind Module option was selected from the Wind Fields options. The diagnostic model allows certain features of the flow field such as sea breeze circulation with return flow aloft (^{Scire et al. 2011}). After that, the CALMET preprocessor and CALPUFF model were run. The final step in the CALPUFF modeling system was the CALPOST postprocessor which is used to calculate pollution concentrations for selected averaging times. Both AERMOD and CALPUFF models were run for 1 and 24 h, and annual averaging times.

Fig. 2 Wind rose plot for the modeling area (2007-2011). 

2.4 Determination of model performance

Fractional bias (FB) (Eq. 3), index of agreement (IA) (Eq. 4), geometric mean bias (MG) (Eq. 5), and geometric mean-variance (VG) (Eq. 6) are widely used statistical parameters for the determination of model performance in air quality modeling studies (^{Chang and Hanna, 2004}; ^{Donnelly et al., 2009}; ^{Paschalidou and Kassomenos, 2009}; ^{Behera et al., 2011}; ^{Amoatey et al., 2018}). A perfect model has IA, MG and VG = 1, and FB = 0 (^{Jeong, 2011}). The FB value falls in the interval of ENT#091;-2, 2ENT#093;. Negative FB represents the over-predicting model whereas positive FB means the under-predicting model (^{Behera et al., 2011}). ^{Hanna and Chang (2012)} stated that an acceptable model should have -0.3 ≤ FB ≤ 0.3 for rural areas and -0.67 ≤ FB ≤ 0.67 for urban areas. IA can be defined as the normalized measure of model prediction errors and varies between 0 and 1. Higher IA values (close to 1) represent an acceptable performance (^{Moriasi et al., 2007}). IA ≥ 0.5 is considered a good performance (^{Amoatey et al., 2018}). Acceptance criteria for geometric mean and geometric variance are 0.7 < MG < 1.3 and VG < 1.6, respectively (^{Chang and Hanna, 2004}; ^{Donnelly et al., 2009}).

where O is observed concentrations, P are predicted concentrations, and O̅ and P̅ are the means of observed and predicted concentrations. Finally, n is the amount of data. To determine the model performance, daily averages of PM₁₀ concentrations measured at Zonguldak AQMS between 2008 and 2011 were used in this study. For that period, this station’s records were published as raw and un-validated measurements. Therefore, outliers were removed from the dataset by examining the box plots. A discrete receptor was placed to obtain concentration values for the AQMS. After that, the only missing information was the background PM₁₀ concentration value (the pollution due to natural sources and far anthropogenic sources). For this, the histogram of PM₁₀ concentrations was drawn, where the peak point of the histogram represents the background concentration. Finally, model predictions plus background concentration were validated against measured concentrations.

3.1 Emission inventory results

Before employing air quality modeling it is necessary to represent the emission inventory results. The results of the PM₁₀ emission inventory are given in Figure 3, according to which a total of 2710.2 t of PM₁₀ emissions were released into the Zonguldak atmosphere annually, and CATES thermal power plant was found as the main source of particulate pollution in the study area. Total emissions of the power plants (CATES and ZETES II) were calculated as 1648.7 t year^-1 (60.8% of total emission). On the other hand, domestic heating (together with harbor emissions) accounted for 38.6% of the total emission. Domestic heating emissions of PM₁₀ for Zonguldak, Kilimli, Kozlu, Çatalağzı, Muslu, and Gelik were determined as 531.5, 151.9, 283.4, 43.8, 10.0, and 16.7 t year^-1, respectively. Motor vehicle emissions were only a small portion (0.6%) of the total emissions. There is only one previous study in the literature (^{Tecer et al., 2008}a) that provides a PM₁₀ emission inventory for Zonguldak, in which the total PM₁₀ emission amount was reported as 17170 t year^-1. According to this inventory, 11000 t of PM₁₀ resulted from the emission of the Erdemir steel factory, but this source was not included in our study since it is approximately 37 km outside the modeling domain. ^{Tecer et al. (2008a)} included all registered cars (with total emissions of 315 t year^-1) in their inventory. On the other hand, they only included 11 main roads regarding traffic emissions, together with the emissions of 55000 homes and workplaces (which were calculated as 355 t year^-1). In our inventory, 93568 homes and 156 public and private workplaces were included. Moreover, they did not address the details of the inventory (such as emission factors), therefore it is not possible to make an exact comparison between the two inventories.

Fig. 3. Emission inventory results for PM₁₀ (t year^-1) for point (CATES and ZETES II power plants), line (motor vehicles), and area (Zonguldak, Kilimli, and Kozlu counties, Çatalağzı, Gelik, and Muslu districts, and harbor) sources. 

3.2 Air quality modeling results

Having represented the emission inventory results, this section provides the results of air quality modeling. Maximum concentrations are obtained from AERMOD and CALPUFF models for 1 and 24 h, and 1-year averaging times. Table VI presents the maximum concentrations and their locations. As seen from this table, the CALPUFF model outputs are much higher than AERMOD. The maximum PM₁₀ concentrations predicted by CALPUFF are 743.2, 236.5, and 37.3 µg m^-3 for 1 h, 24 h, and 1-year averaging times, respectively. On the hand, AERMOD gives 364.2, 105.1, and 17.2 µg m^-3 values for the same averaging times. In literature, several studies compotherared AERMOD and CALPUFF models with respect to model performance. ^{Tartakovsky et al. (2016)} performed air quality modeling in complex terrain and suggested using AERMOD instead of CALPUFF for regulatory purposes to better protect the population. However, their modeling area was located to the west of Jerusalem and away from the shoreline. ^{Amoatey et al. (2018)} conducted a health risk assessment study in a coastal area with flat terrain and reported that AERMOD predictions were better than CALPUFF. On the other hand, according to the results of this study, CALPUFF predictions were higher than those of AERMOD. The higher concentrations of the CALPUFF model may result from recirculation of pollutants. Sea breezes (also called day breezes) and thermal internal boundary layers are the sources of the recirculation process; they occur because of the thermal contrast between dry land and sea (^{Efimov and Barabanov 2009}; ^{Shang et al. 2019}). ^{Levy et al. (2009)} mentioned that coastal recirculation adversely affects air pollution concentrations. ^{Indumati et al. (2009)} stated that sea breeze could increase ground-level concentrations near the coastal regions. Although CALPUFF does not include a specific sea breeze module, a diagnostic wind module option was used in this study. Therefore, the effects of sea breezes can be determined by the CALPUFF modeling system. Moreover, the CALPUFF model has a specific algorithm to handle overwater dispersion and coastal interactions (coastal fumigation and thermal internal boundary layer), so that the impacts of coastal fumigation can be revealed easily. The complex topography and coastal fumigations create complex meteorological conditions which can be handled more easily by the CALPUFF model (^{Abdul-Wahab et al., 2011}). ^{Bluett et al. (2004)}) suggested that using advanced models like CALPUFF for coastal areas can provide a better simulation of the fumigation process. As a result, it can be said that AERMOD underpredicts concentrations as compared to CALPUFF on shoreland with complex topography, while CALPUFF gives better results on complex topographical conditions for modeling emissions from gas refineries located near the shoreline (^{Atabi et al. 2016}). Therefore, advanced models like CALPUFF should be preferred over AERMOD for environmental impact studies of emission sources located in seashore and hilly regions.

The initial objective of this study was to determine polluted zones by air quality modeling. PM₁₀ distribution maps are displayed in Figure 4a-d for 1 and 24 h averaging times for both models. It is apparent from Figure 4 that CALPUFF predicted higher concentrations as compared to AERMOD, which reports the maximum concentration (364.2 µg m^-3) at Fatih district of Kozlu county for 1 h averaging time, while CALPUFF predicted the maximum concentration as 743.2 µg m^-3 at Subatan district of Kilimli county, west of Gelik, for the same averaging time. As seen from Figure 4a, the most polluted region is Kozlu, and polluted areas are close to the seashore in AERMOD. However, CALPUFF gives results (Fig. 4b) for different hotspot areas where PM₁₀ concentrations are higher than 400 µg m^-3. ^{Kallos et al. (2000)} investigated sea breezes on the Black Sea and emphasized that strong sea breezes are observed over the Turkish coast during the day. They also reported that urban plumes of Istanbul disperse to the SW direction. The location of Zonguldak is close to Istanbul and similar wind patterns can be expected. The thermal power plant emissions in the modeling domain may disperse to the central region, explaining the higher concentrations and hotspots of particulate matter emissions in CALPUFF results. In Figure 4c, 24-h averaging time concentrations of AERMOD show similar patterns as in the 1-h averaging time concentrations (Fig. 4a). The maximum concentration (105.1 µg m^-3) occurs at Zonguldak harbor and the surrounding region (city center). For the same averaging time, the CALPUFF model’s peak concentration is seen at the Terakki district of Zonguldak county (236.5 µg m^-3) (Fig. 4d). The daily PM₁₀ limit established by the World Health Organization (WHO) and the European Union (EU) is 50 µg m^-3 (^{Zeydan and Pekkaya, 2021}). As seen from Figure 4c, d, the highest concentrations are above this threshold value.

Fig. 4 PM₁₀ pollution distribution maps. (a) AERMOD model output for 1-hour averaging time. (b) CALPUFF model output for 1-hour averaging time. (c) AERMOD model output for 24-hours averaging time. (d) CALPUFF model output for 24-hours averaging time. 

The air quality modeling results of this study are especially important in terms of the Çatalağzı region. The CATES thermal power plant was installed in 1948 and renewed in 1991. However, due to its old technology, CATES could not meet air quality regulations and it was on the focus of environmentalists. In the last two decades, this region was considered as an energy basin by the Ministry of Energy and Natural Resources. ZETES I and II were installed in 2010 in the Çatalağzı region. It is obvious from the modeling results that PM₁₀ emissions of thermal power plants together with domestic heating emissions create unhealthy ambient air conditions for the people in Zonguldak. According to the recent report of the Health and Environment Alliance (HEAL), emissions of CATES and ZETES II thermal power plants may be related to 1100 and 1412 premature deaths, respectively (^{HEAL, 2021}), suggesting the need of taking urgent precautions for public health. It should be also noted that ZETES I was under maintenance and was not included in air quality modeling for the year 2011. However, in 2016, another thermal power plant called ZETES III (700 × 2 MW) was installed close to the other plants and the total capacities of the ZETES power plants reached 2790 MW, making them the largest in Turkey with respect to capacity. Currently ZETES (I, II, and III) and CATES produce 7% of the total electricity of Turkey (^{Cardoso and Turhan, 2018}). Perhaps the most important finding of this study is that the Çatalağzı region should be considered as polluted in terms of particulate pollution and no more thermal power plants should be installed in this region.

As stated earlier, this work is the first detailed air quality modeling study in the Zonguldak province that shows the spatial distribution of PM₁₀ concentrations. So, it is not possible to compare the results of this paper to previous studies in the same modeling area. Therefore, we examined the maximum PM₁₀ concentrations of air quality modeling studies of other cities in Turkey. ^{Elbir et al. (2010)} used the CALPUFF model and reported the maximum annual PM₁₀ concentration as 90 µg m-3 for Istanbul. In another study, the CALPUFF model predicted the maximum annual PM₁₀ concentration as 276 µg m^-3 in Izmir (^{Elbir, 2004}). Both cities are highly populated and heavily industrialized as compared to Zonguldak. Considering that CALPUFF estimated the maximum annual PM₁₀ concentration as 37.3 µg m^-3, it could be said that the annual maximum PM₁₀ levels in Zonguldak province are not as high as that of metropolitan cities like Istanbul and Izmir.

3.3 Model performance evaluation

The last but not least issue regarding air quality modeling studies is model performance evaluation. Model-predicted plus background concentrations are validated against daily measured PM₁₀ concentrations at Zonguldak AQMS between 2008 and 2011. The total dataset contains 1186 observations after the elimination of outliers. The background concentration was not given as an input for both models since it was aimed to predict PM₁₀ pollution distribution from known sources in the modeling area; nevertheless, as stated earlier, some PM₁₀ sources (illegal coal use, quarries, coal transport, storage facilities, and ash dams of thermal power plants) could not be included in the inventory study. Moreover, re-suspension of particulates from coal ash refuses by the effect of wind may have contributed to PM₁₀ concentrations. Also, some portions of particulates may be transported to the modeling domain via long-range transport. Therefore, background concentration should be included in the model performance evaluation. Unfortunately, there was no background AQMS in the modeling area. So, background concentration was obtained by drawing a histogram of PM₁₀ concentrations. The peak point of the histogram corresponds to background concentration. As seen from Figure 5, background concentration is determined as 49 µg m^-3.

Fig. 5 Histogram of PM₁₀ concentrations (µg m^-3). 

Table VII compares the model evaluation parameters for AERMOD and CALPUFF. FB values are found to be 0.30 for AERMOD and 0.28 for CALPUFF, both positive, which means that the two models underpredicted PM₁₀ concentrations. ^{Lee et al. (2014)} found FB values ranging from 0.043 to 0.821 for PM₁₀. ^{Barna and Gimson (2002)} reported a fractional bias ranging between -0.17 and 0.30 for PM₁₀ for their dispersion modeling in Christchurch, New Zealand. ^{Gibson et al. (2013)} found annual, monthly, and hourly FB values for PM_2.5 as 0.96, 0.88, and 0.89 for Halifax and 0.96, 0.95, and 0.90 for Pictou, respectively. ^{Otero-Pregigueiro et al. (2018)} found FB values of 0.44 in their PM₁₀ modeling study near an industrial plant. The FB values in this study are quite reasonable compared to the literature, falling into the acceptable range of ^{Hanna and Chang (2012)} definition for urban areas (-0.67 ≤ FB ≤ 0.67). In our study, IA values were found to be 0.50 for AERMOD and 0.51 for CALPUFF, passing a good performance criterion (IA ≥ 0.5) (^{Amoatey et al., 2018}). ^{Elbir et al. (2010)} reported lower IA values (between 0.33 and 0.50) for an air quality modeling study for PM₁₀ in Istanbul. On the other hand, similar IA findings (in the range of 0.37-0.59) were reported by ^{Elbir (2004)} for particulate pollution in Izmir. ^{Lee et al. (2014)} reported IA values for PM₁₀ between 0.50 and 0.85 in Ulsan, Korea. Therefore, it can be concluded that IA values in this study are consistent with the literature.

In this study, geometric mean bias (MG) values are determined as 1.21 for AERMOD and 1.20 for CALPUFF. ^{Song et al. (2006)} reported MG values in the range of 1.01-1.83 for PM₁₀ modeling in Beijing. MG findings in another PM₁₀ modeling in the greater Thessaloniki area, Greece was stated as 1.17 (^{Kakosimos et al., 2011}). MG values in this work satisfy the model acceptance rule (0.7 < MG < 1.3) and are in accordance with other studies. The last statistical parameter of performance evaluation in this study is geometric mean-variance (VG), whose values for both models were calculated as 1.24, also below the limit value (VG < 1.6). ^{Kakosimos et al. (2011)} reported VG values of 1.07 for PM₁₀ modeling in the greater Thessaloniki area. VG values for SO₂ and NO2 modeling in Tema metropolis, Ghana were found as 1.01 and 2.61 for CALPUFF, and 2.04 and 1.04 for AERMOD models (^{Amoatey et al., 2018}). Thus, FB, IA, MG, and VG values in this study are in line with the literature and all parameters fall in acceptable ranges. Therefore, it can be concluded that both models are satisfactory, but the performance of CALPUFF is slightly better than that of AERMOD due to lower FB, MG, and VG values, and higher IA value.