Benzo[a]pyrene in the Ambient Air in the Czech Republic: Emission Sources, Current and Long-Term Monitoring Analysis and Human Exposure

: This paper provides a detailed, thorough analysis of air pollution by benzo[a]pyrene (BaP) in the Czech Republic. The Czech residential sector is responsible for more than 98.8% of BaP, based on the national emission inventory. According to the data from 48 sites of the National Air Quality Monitoring Network, the range of annual average concentration of BaP ranges from 0.4 ng · m − 3 at a rural regional station to 7.7 ng · m − 3 at an industrial station. Additionally, short-term campaign measurements in small settlements have recorded high values of daily benzo[a]pyrene concentrations (0.1–13.6 ng · m − 3 ) in winter months linked to local heating of household heating. The transboundary contribution to the annual average concentrations of BaP was estimated by the CAMx model to range from 46% to 70% over most of the country. However, the contribution of Czech sources can exceed 80% in residential heating hot spots. It is likely that the transboundary contribution to BaP concentrations was overestimated by a factor of 1.5 due to limitations of the modeling approach used. During the period of 2012–2018, 35–58% of the urban population in the Czech Republic were exposed to BaP concentrations above target. A signiﬁcant decreasing trend, estimated by the Mann-Kendall test, was found for annual and winter BaP concentrations between 2008 and 2018.


Introduction
Polycyclic aromatic hydrocarbons (PAHs) are ubiquitously distributed in the environment [1]. They are common by-products of combustion processes of fossil fuels and wood. PAHs represent a group of substances, many of which have toxic teratogenic, mutagenic or carcinogenic properties [2,3]. They affect fetal growth. Prenatal exposure to PAHs is related to markedly lower birth weight [2] and probably also has negative effects on the cognitive development of young children [4]. Due to their physical and chemical properties, all these substances can be transported over long distances and deposited in remote areas [5][6][7]. PAHs can bioaccumulate, enter the food chain [1] and be toxic to the environment.
Benzo(a)pyrene, occurring in the atmosphere primarily bound to particulate matter, has been set as a suitable marker of ambient air pollution caused by PAHs. A European directive has set a target value of 1 ng·m −3 for the total content of BaP in the PM 10 fraction, averaged over a calendar year [8], with the aim of avoiding, preventing or reducing harmful effects on human health and/or the environment as a whole. The WHO has not drafted a guideline for BaP, which is a potent carcinogen. The reference level of 0.12 ng·m −3 was estimated assuming the WHO unit risk for lung cancer for PAH mixtures and an acceptable risk of additional lifetime cancer of approximately 1 × 10 −5 [9,10].
Residential combustion as an important source of PAHs and other air pollutants are responsible for the majority of anthropogenic emissions of BaP in Europe. Such emissions are linked to adverse health effects, especially in urban and suburban areas where emissions and population densities are higher [11]. A modeling study of Europe [11] stated that it is necessary to assess concentrations of BaP in Europe, as an indicator for PAHs, and quantify their health-related effects. The European Environmental Agency estimates that in 2017, 17% of the EU−28 urban population was exposed to above-target annual mean BaP concentrations; this is the lowest value since 2008. As in previous years, values above 1 ng·m −3 are predominantly found in Central and Eastern Europe. The highest concentrations were recorded mainly at stations in Poland and the Czech Republic [12].
Air pollution by BaP is one of the main problems associated with ensuring air quality in the Czech Republic. BaP concentrations exhibit significant intra-annual variation with maxima in winter that are related to emissions from seasonal anthropogenic sources (local heating units) and generally worsened dispersion conditions.
The aim of this study is to assess the current levels and long-term trends of air pollution by benzo(a)pyrene in the Czech Republic, together with their causes. A transboundary contribution to the annual mean concentrations of BaP is also quantified via a modeling-based approach.

Area Description
The Czech Republic is located in Central Europe. The topography of the Czech Republic consists predominantly of hills and highlands. More than two thirds of the territory is located below 500 m altitude, with the large majority of the rest between 500 m and 1000 m and only 1% above 1000 m.
The climate in the Czech Republic is mild; it could be classified as somewhere between continental and maritime. It has 4 seasons. Local variations of weather and climate are influenced by ruggedness and altitude. The average annual air temperature in the Czech Republic usually varies from 5.5 to 9.5 • C. The coldest months of the year are December, January and February. The hottest are July and August. Usually, precipitation in the Czech Republic is at a maximum in July and a minimum in February. Currently, about 60% of the population lives in cities with more than 5000 inhabitants [13,14].

Sampling and Analytical Methods
In the Czech Republic, concentrations of BaP in the PM 10 fraction as measured at manual monitoring stations form the basis for the evaluation of air quality. The monitoring stations are placed mainly in cities and in areas with known high BaP concentrations ( Figure 1). In 2018, BaP concentrations were measured at 48 monitoring sites (Table A1). The majority of stations are located in cities, with 28 urban and suburban stations. Transport and industrial contributions to BaP are monitored at six traffic and nine industrial stations. Background levels of BaP concentrations are monitored at 5 rural monitoring stations. PM 10 samples are taken by low or high volume samplers on a quartz filter. The samples are processed in certified chemical laboratories and analyzed by high pressure liquid chromatography (HPCL) or gas chromatography with mass detection (GC/MS). The measured concentrations of BaP are daily averaged value and are collected with a minimum three-day frequency. The concentrations measured at the pollution monitoring stations are stored in the Czech national Air Quality Information System (AQIS) database. The lower detection limit is 0.04 ng·m −3 for GC/MS and 0.10 ng·m −3 for HPLC. The measurement uncertainty of BaP is up to 25%. Short-term monitoring of BaP concentrations in small settlements has been made in campaigns within the Czech national project TITSMZP704-Measurement and Analysis of Air Pollution with Emphasis on the Evaluation of the Share of Individual Groups of Sources. These are case studies that monitor the variability of short-term BaP concentrations during the heating season under the specific conditions of particular small settlements. As this is an ongoing project finishing in 2021, here we will present only a sample of data from 3 small settlements ( Figure 1)-Bochovice,Černíny and Hřivice-to show the level of BaP concentrations in villages where they are not regularly monitored and where solid fuel heating is predominant. Bochovice is a small village with 143 inhabitants with 58 houses, which are heated by solid fuels. InČerníny there are 370 permanent residents and 140 houses out of 167 are heated by coal or wood. Hřivice is a village with 631 inhabitants, with 251 houses heated by coal or wood. We also present the Kladno-Švermov locality as an example of an area surrounding a current monitoring station measuring very high levels of air pollution caused by local heating. Kladno-Švermov is a district to the north of Kladno city situated in a shallow valley, where almost 5000 people live. It has a high building density with both central and local heating.

Emission Calculation
BaP is mainly a product of incomplete combustion of organic substances at temperatures between 300 and 600 • C. Unsurprisingly, the main contribution to total BaP emission in the Czech Republic is from combustion of solid fuels in low-capacity boilers. This is mainly household fuel combustion for heating, cooking and water heating (residential sector).
BaP emissions from the residential sector are calculated on the basis of emission factors for various combinations of fuel type and installed combustion plant (see Tables A2-A5 in Appendix B). The combinations used are set from the annually updated distribution of solid fuel and type of heating equipment (Table A5 in Appendix B), resulting in country specific emission factors which reflect the particular circumstances for a given year. For the purpose of national emission inventory, the total fuel consumption in households is determined by the Czech statistical office [13]. The emissions in this model are calculated as a sum for the whole country and comprise local heating, cooking and water heating. This national model uses emission factors at nominal heat output.
For territorial distribution of residential heating emissions, a bottom-up up approach is used. The emissions are calculated for each basic territorial unit-municipalities and city districts-and comprise only the local heating from permanently occupied households. The fuel consumption at a basic territorial unit is calculated based on the average annual heating amount and specific fuel type consumption per average housing type. The base data are obtained from the 2011 Population Census (number of households, their type of heating and average floor area) and the results of the ENERGO 2015 statistical survey (share of given fuel burned in a particular installation type, share of insulated/noninsulated flats). The year-by-year changes in fuel consumption are mostly influenced by the characteristics of the heating period, which is expressed as the number of heating degree days (the sum of the differences between the reference indoor temperature and the average daily outdoor temperature on heating days). Other regional annually updated parameters are the number of households and their type of heating. Solid fuel parameters and the share of solid fuel consumption according to the installation type of the combustion plant are annually updated at a national level, based on the statistics of boiler sales ascertained by the Ministry of Industry and Trade and data from the subsidy program for the boiler replacements. Solid fuel parameters are updated according to the results of the survey on supplies and quality of solid fuels in the Czech Republic carried out by the TEKO company [15]. The regional calculation model uses a 15/85 boiler operating mode, i.e., 15% of time at nominal heat output and 85% at lower heat output. This assumption is in accordance with the Ecodesign Directive 2009/125/EC. Calculation model for the territorial distribution of emissions is more sensitive to climatic conditions in a given year than model for the national emission inventory. The highest emission difference between these two models is in 2014 and 2018, when the heating period was mild and short (see Table A6 in Appendix B).
Emission factors of solid fuels for local heating were obtained from measurement results of the most common fuel and boiler type combinations used for household heating in the Czech Republic. These measurements were carried out at the Energy Research Center of VSB-TUO [16] during the years 2008-2016. Emission factors for liquid and gas fuel were taken from the Emission inventory guidebook [17].
The transport sector comprises emissions from road transport, railways, air and water transport, off-road transport used in agriculture, forestry, building construction and area transport within large industrial enterprises. Emissions are calculated at the Transport Research Center [18] in an up-to-date version of the COPERT program [19] based on nationwide fuel consumption.
BaP emissions from combustion in nonresidential stationary sources such as power and heat generation, combustion processes in industry and manufacturing, institutions and services and waste incineration are calculated from activity data and given emission factors; national emission factors are estimated for a particular sector or taken from the Guidebook [17]. Emissions from industrial sources are either reported by the operator or calculated from activity data and emission factors.

Spatial Mapping and Population Exposure
The methodology used for the creation of the BaP concentration maps is a linear regression model followed by an interpolation of its residuals; rural and urban areas are mapped separately and then merged by population density [20]. The methodology is referred to as the Regression-Interpolation-Merging Mapping (RIMM) method and is used for air quality mapping in the Czech Republic and elsewhere in Europe as well [11,21,22]. The estimate of concentrations is calculated using the relationship: whereẐ is the estimated concentration value at point s 0 , X i are the various supplementary data, c and a i are the parameters of the linear regression model andR is the spatial interpolation of the residuals of the linear regression model at point s 0 , calculated on the basis of the residuals at the points of measurement. The primary data for creating air pollution maps of BaP are concentrations measured at individual monitoring stations. Since there are only a limited number of monitoring stations and their spatial representativeness is variable, various supplementary (secondary) data are also used. These secondary data both provide comprehensive information about the entire territory and also exhibit regression dependence on the measured data. The main secondary sources of information are outputs of dispersion models, which combine data from emission inventories and meteorological data. In the Czech Republic, the secondary data mainly used are annual mean concentrations provided by EMEP/MSC-E [23] together with annual mean concentrations from the Czech Gaussian model SYMOS. Other supplementary data can be provided by maps of annual mean PM 10 and PM 2.5 concentrations.
The kriging and Inverse Distance Weighting (IDW) techniques are used as interpolation methods [24]. Interpolation of residuals using IDW is calculated using the relationship: whereR is the estimate of the field of residuals at point s 0 , R(s i ) is the residual of the linear regression model at the measuring site s i , N is the number of surrounding stations used in the interpolation, d 0i is the distance between points s 0 and s i , and β is the weight. In case of ordinary kriging, the interpolation of the residuals is calculated using the relationship: where R(s i ) is the residual of the linear regression model at the measuring site s i and λ i are estimated weights based on the theory of spatial statistics [24] derived from a variogram fitted to an empirical variogram 2γ ν (h) of the field of residuals. The variogram expresses the dependence of the variability between points on the mutual distance between the points and the empirical variogram is calculated as follows: where R are the residuals at measuring points s i and s j , d ij is the distance between points s i and s j , n is the number of pairs of stations s i and s j whose mutual distance is h ± δ, and δ is the tolerance. The calculated urban and rural map layers are subsequently merged by a layer of population density α:Ẑ whereẐ is the final estimate of the concentration at point s 0 ,Ẑ r ,Ẑ u is the concentration for the rural or urban map layer, and α 1 , α 2 are the classification intervals corresponding to the population density. For the BaP concentration maps α 1 was set to 200 inhabitants per km 2 and α 2 was set to 1000 inhabitants per km 2 . The entire concept of separate mapping of rural and urban pollution is based on the assumption thatẐ r (s 0 ) ≤Ẑ u (s 0 ) for BaP. For areas where this assumption is not fulfilled, a third layer created in a similar fashion to the urban and rural layers is used; this third layer is created using all the background stations without distinguishing between urban and rural stations.
The maps are constructed with a spatial resolution of 1 × 1 km 2 . The uncertainty of the map was assessed using the cross-validation method: concentration at the location of a measuring site is always estimated from other station data only, thus providing an objective estimate of the map quality away from measurement site locations. In this article, the uncertainty of the maps is expressed by the relative root-mean-square error (RRMSE): where Z is the measured value of the concentration at point s i ,Ẑ is its estimate using cross-validation and N is the number of measuring stations. The spatial distribution of the uncertainty was not estimated. It should be noted that the cross-validation is applied only during the interpolation of residuals; parameters of linear regression are always estimated using all the stations. Therefore, the overall uncertainty of the maps is somewhat underestimated. The uncertainties (RRMSE) were calculated for each map layer separately and were up to 30% for urban and over 60% for rural areas.
The higher uncertainty of rural areas is due to lack of measurements at rural regional stations and the absence of more extensive measurements in smaller settlements in the Czech Republic. The annual mean BaP concentration maps 2012-2018 were prepared at CHMI during the annual air quality assessments. Estimation of population exposure to above-target BaP concentrations were calculated based on maps of BaP and population density data with resolution 1 × 1 km [13].

Trend Analysis
Trends of annual average BaP concentrations were analyzed at six selected monitoring stations in the Czech Republic ( Figure 1). Five stations were classified as urban or suburban, the remaining station was classified as a rural regional station. Station selection was based on their classification and the quantity of data availability for trend analysis. We focus on data from urban and suburban monitoring sites since one of the aims of this study is to assess human exposure to BaP concentrations. For a comprehensive overview, the data from the Košetice rural regional site are also presented.
Trends for annual, winter (October-March) and summer (April-September) average concentrations are analyzed. The authors annually prepare average monthly concentrations of BaP for the "Air Pollution in the Czech Republic" yearbooks [22] (the newest report) and have partitioned this data into winter and summer periods. From April to September, average monthly BaP concentrations are usually below or just above the target value while for the rest of the year they very often exceed the target value.
Temporal trends, i.e., annual averages of BaP concentrations and emission, were analyzed using the nonparametric Mann-Kendall trend test with a level of significance of 0.05 [25,26]. This test is among the most widely used statistical methods for this kind of data [27][28][29][30][31] and is particularly useful since it tolerates missing values and the data need not conform to any particular distribution. Moreover, as only relative rather than absolute magnitudes of the data are used, this test is less sensitive towards incomplete data capture and special meteorological conditions leading to extreme values [32] that often affect air quality data.
If a linear trend is significant, the slope and severity of the trend is estimated by Sen's test [32]. For all stations, annual average concentration and emission were analyzed for the same 2008−2018 time period. There were no missing averaged annual data.
We used R-Studio software for statistical analyses [33]. The maps were created using the geographic information system ArcGIS by ESRI [34].

Source Apportionment
During the update of the National air quality plans (NAQP) in 2018, a transboundary contribution to annual mean concentrations of BaP in 2015 had to be quantified. It is known that polycyclic aromatic hydrocarbons including BaP undergo gas-particle partitioning and degradation in the atmosphere, but these processes are not fully understood and other processes, e.g., secondary organic aerosol coating can protect PAHs from ozone degradation during long range transport [35]. Due to the limited time for the update of NAQP and unavailability of ready-to-use models for PAHs long-range transport, a chemical transport model CAMx v5.41 [36] was adopted to account for BaP as a passive tracer. The limitations of this approach are discussed further in the text.
The CAMx model was run in two nested domains d01 and d02 with 14.1-km and 4.7-km resolution respectively ( Figure 2). The transboundary contribution was estimated with a brute-force method. Sources outside the Czech Republic were set to zero. The spatial distribution of concentrations originating from Czech sources within the 4.7-km CAMx grid was determined by the Gaussian model SYMOS [37] at 0.5-km resolution: where C CZ is the contribution of Czech sources calculated by the CAMx model in a 4.7-km grid, S(i) is total contribution of Czech sources in subgrid point i calculated by the SYMOS model and n is number of subgrid points. Next, a relative contribution of sector C of Czech sources was determined: where C nonCZ is the contribution of sources outside of the Czech Republic calculated by the CAMx model in a 4.7-km grid, and S c (i) is the contribution of the Czech sources (sector C only) in subgrid point i calculated by the SYMOS model. Meteorological inputs with a 1-h time step were derived from the assimilation cycle of the numerical weather prediction model ALADIN/CE version ALARO [38] operated at the CHMI at 4.7-km resolution and with 87 vertical levels. In the assimilation cycle, the analysis at 0, 6, 12 and 18 UTC was followed by a 6-h forecast. Analysis of upper-air parameters combines the driving model ARPEGE with mesoscale structures of the ALADIN model through DFI blending complemented by 3DVAR assimilation of observations [39]. Analysis of surface temperature and relative humidity is based on optimal interpolation and serves as an input to the Interaction Soil Biosphere Atmosphere (ISBA) scheme describing exchanges between the atmosphere and the land surface [40]. The 68 lowest ALADIN levels were aggregated into 26 CAMx levels, with the top of the lowest level at approximately 50 m and the highest level at approx. 10 km above ground.
High-resolution BaP emissions for the Czech Republic were taken from the calculation model for territorial distribution. For the Polish Silesian and Lesser Poland Vovoidships, BaP emissions were estimated by the ATMOTERM company within the LIFE-IP MAŁOPOLSKA project (LIFE14 IPE/PL/000021). For Slovakia, emissions from the SNAP 2 sector were taken from the national emission inventory. For other sectors in Slovakia and the rest of the modeling domain, a top-down emission inventory based on the TNO emission inventory [41] was used: emissions for the year 2015 were estimated by linear interpolation and then distributed to ten sectors following the SNAP nomenclature. Industrial emissions were allocated to sources registered in the European point source emission register E-PRTR [42] and remaining industrial emissions as well as emissions from area sources were distributed using the "industrial area" land cover class from the Corine Land Classification (CLC) database [43]. For the spatial distribution of residential heating, a combination of population density and urban area was used, assuming different fuel mixes in metropolitan and rural areas. Typically, rural areas exhibit a larger per capita emission of BaP due to the increased usage of coal and wood. A map of annual emissions used is shown in Figure 2.
The time distribution of the residential heating emissions in the Czech Republic and Silesian and Lesser Poland Vovoidships was based on temperature profiles, otherwise factors for month, day of the week and hour of the day were used [44,45]. For the vertical distribution of point source emissions from top-down inventory, typical point source parameters based on the analysis of the data from the Czech database were used for each SNAP category.

Emissions of Benzo(a)pyrene
In 2018, the residential sector accounted for more than 98.8% of total Czech BaP emissions (15.56 t out of total 15.74 t). The remaining percentage was produced mainly by the transport sector (0.78%, 122 kg), especially passenger cars (0.48%, 76 kg). The industrial share of total BaP emissions was 0.27% (43 kg) and the share of institutions and services 0.10% (16 kg). The emissions from industrial sources are year-round, in contrast with local heating, especially in the Moravian-Silesian region. There is also a higher proportion of coal combustion in households in this region, which is reflected in a higher BaP emission load. The contributions from the different sources of BaP emissions have not changed significantly in the last 10 years. Almost 58% of BaP emissions from the residential sector are produced by combustion of fuel wood, and the consumption of wood is still increasing (Table A2 in Appendix B). In contrast, coal consumption is decreasing in recent years, which is reflected in its decreasing share of emissions ( Figure 3). In 2018, approximately 15% of households were using solid fuels for local heating, and these households are responsible for more than 98% of total BaP emissions in the Czech Republic. Figure 4 presents BaP emissions from combustion of solid fuel in the residential sector in 2018 sorted by combustion installation and fuel type. More than 50% of total BaP emissions from the residential sector are estimated to be from fuel combustion in over-fire boilers and about 30% from fuel combustion in fireplaces and stoves.

Ambient Air Concentrations of BaP and Population Exposure
The map of annual average concentration of BaP in 2018 is shown in Figure 5. Areas where BaP concentrations were higher than the target value of 1 ng·m −3 (above-target) are indicated in red or brown in the figure. The thresholds correspond to the upper and lower assessment thresholds of 0.6 ng m −3 and 0.4 ng m −3 , the target value of 1.0 ng m −3 set by EU legislation [8], and to 2.0 ng m −3 to distinguish the most polluted areas in the Czech Republic. High values of BaP concentrations were estimated in the North East of the Czech Republic, referred to as the Ostrava region. Other contaminated areas include the Kladno district (West of Prague), areas of Prague and a number of smaller towns. Areas with above-target concentrations comprised 12.6% of the Czech territory in 2018. The lowest annual average concentrations of BaP were estimated to be in locations distant from emission sources and therefore free from direct exposure (i.e., natural mountain areas). The range of measured concentration of BaP in 2018 was from 0.4 ng·m −3 at the Košetice rural station to 7.7 ng·m −3 at the Ostrava-Radvanice industrial station.  The highest daily average concentration of BaP over the sampling campaign-24.5 ng·m −3 -was observed in the Kladno-Švermov station whereas the lowest daily average concentration of BaP (0.1 ng·m −3 ) was recorded inČerníny. The limited amount data obtained by these campaign measurements, which were only obtained in winter, does not allow for calculation of annual average concentrations. Nevertheless, in Bochovice andČerníny the target limit value (1 ng·m −3 ) established by European legislation was exceeded on 59% and 54% of measurement days, respectively. In contrast, the average daily BaP concentrations monitored in Hřivice and Kladno-Švermov were below the target limit in only one case.    In contrast, the highest annual average value of 0.7 ng·m −3 at the Košetice rural regional site, representing the wider area without the local emission, was measured in 2013. The lowest annual average value and the median of 0.4 ng·m −3 was recorded for six years (2008,2011,(2014)(2015)(2016)2018). The median value for the Košetice site is 0.5 ng·m −3 .

BaP Concentration and Emission Trends
With respect to the main emission source of BaP (Figure 3), i.e., local heating causing higher levels of BaP in ambient air, we assessed concentration trends specified for the winter (October-March) and summer (April-September) period and for the year as a whole (Figures 8-10). More detailed graphs presenting the course of annual, winter and summer concentrations using boxplots at each measuring station can be found in Appendix D, Figure A2.
Around 65% of annual average concentrations from suburban and urban stations during the study period were higher than the target value of 1 ng·m −3 [8]. For the sake of completeness, this figure was around 96% for the winter periods and 7% for the summer period, respectively. Only one rural regional site registered no above-target concentration during the entire study period.
Comparing the annual average values for all stations between 2008 and 2018, there is a decrease of 27% in the BaP annual average concentration and a 24% decrease for winter and 50% decrease for summer average concentrations (Figures 8-10). For urban and suburban stations, the situation is very similar with a 28%, 25% and 49% decrease for the annual, winter and summer periods respectively (Figures 8-10).   The Mann-Kendall test was used to assess the monotonic trend of BaP at selected monitoring sites ( Table 2). Concerning annual averages, significant decreasing trends were detected at two stations (Kladno-Švermov with p value = 0.01 and Sen's slope of −0.15 ng·m −3 ·year −1 and Ostrava-Poruba with p value < 0.005 and Sen's slope of −0.12 ng·m −3 ·year −1 ). Concerning winter averages, the significant decreasing trend was detected at the same stations-Kladno-Švermov with p-Value < 0.005 and Sen's slope of −0.28 ng·m −3 ·year −1 , Ostrava-Poruba with p-Value = 0.03 and Sen's slope of −0.21 ng·m −3 ·year −1 .
The average annual and average winter concentrations overall for all stations exhibited a significant decreasing trend. The decrease per year in the BaP concentration was equal to 0.05 ng·m −3 for annual averages and 0.09 ng·m −3 for winter averages, respectively. For urban and suburban stations, the annual and average winter concentrations showed a similar significant decrease is similar-0.06 ng·m −3 for annual averages and 0.11 ng·m −3 for winter averages, respectively. No clear trend was found for BaP summer concentrations.
In terms of total BaP emission for the Czech Republic, 2013 was the year with the highest annual mean value, which was 17.5 t. BaP emissions increased by more than 11% between 2009 and 2010. This increase was due to the implementation of new statistical data for hard coal consumption from 2010. In addition, the years 2010 and 2013 were characterized by long and cold heating season compared to other years. After 2013, BaP emissions had a decreasing trend supported by milder winter seasons (especially 2014 and 2018) but also by decreasing coal consumption and especially the replacement of high-emission boilers. Consequently, no significant trend (p-Value = 0.88) was found for BaP emission development (for more details, see Section 3.1). Moreover, no correlation (p-Value = 0.054) was found between emission and the concentrations from the Košetice rural regional site representing background levels in the Czech Republic.

Source Apportionment
The transboundary contribution to annual average concentration, as estimated by the CAMx model in its 4.7-km resolution grid mode, is 46-70% for most (80%) of the territory of the Czech Republic, with a median value of 57% ( Figure 11). The highest transboundary contribution is modeled in the relatively clean Western and Southern parts of the country and in the North-Eastern mountain regions (cf. Figure 5). When subgrid scaling is applied, the contribution of Czech sources can exceed 80% in residential heating hot spots ( Figure 12). Three categories of Czech sources with a relative contribution to annual average concentration exceeding 10% were identified: residential heating-the absolutely dominant Czech source on most of the Czech territory; road transport-only in large cities Prague and Brno and in vicinity of major roads; and coke oven plants in the Ostrava agglomeration.
When mapped annual average BaP concentration is multiplied by the relative contribution of transboundary sources, regions in the North-East parts of the Czech Republic are still above 1 ng·m −3 , which indicates that the target value cannot be reached without significant reduction of BaP emissions in Poland in hand with measures to mitigate Czech sources.
The reliability of source apportionment results depends of course heavily on the emission inventory and dispersion model used. As stated above, BaP was treated as passive pollutant in CAMx, which can probably lead to overestimation of long range transport due to neglect of its degradation. Nevertheless, emission inventory plays probably the more important role. At the beginning of our modeling effort in 2017, there was only a limited amount of European-scale BaP emission inventory available that could be used in air quality models (authors were aware of [46]; gridded BaP emissions for 2015 were made available by EMEP in December 2017 and marked as unofficial data evaluation purposes only [47]). For this reason we used a top-down inventory based on [41].  Annual average BaP concentration modeled by the CAMx domain d02 was compared with measurements. Annual statistics at station locations were taken from Air Quality e-Reporting [48]. Statistics for the Czech stations were taken from the CHMI's Air Quality Information System, since there was an error in data provided to the AQ e-Reporting. Since annual statistic in AQ e-Reporting do not include information on station classification, all available data marked as valid and verified were used. As can be seen from Figures A3 and A4 and Table A7 in Appendix E, we get the best agreement with observations for the Czech stations (model/observation 0.3-3.6 with median 1.5). For Poland, values are generally underestimated (model/observation 0.1-1.5 with median 0.4), while for Germany and Austria the model largely overestimates measured values (model/observation 4-17.7 with median 5 for Germany and 1.9-26.1 with median 7.8 for Austria). This results in a clear South-West to North-East gradient in model bias. It seems reasonable to expect that the model bias will lead to overestimation of the transboundary contribution in the Southern part of the Czech Republic and to its underestimation in Northern parts (especially in the Ostrava agglomeration). To confirm this assumption, the transboundary contribution to annual average BaP concentrations was compared with data provided by EMEP/MSC-E [23] (transboundary contribution to annual mean concentration provided via personal communication with A. Gusev). For this purpose, CAMx results were aggregated to the EMEP 50-km grid. From Figure A5 we can see that over approximately one third of the Czech Republic the transboundary contribution estimated by this study and by EMEP does not differ by more than 5% (absolute difference). Compared to the EMEP results, the transboundary contribution in the South-central part of the Czech Republic is 10-20% higher, while in the Ostrava agglomeration it is 5-7% lower. Nevertheless it must be noted that EMEP results are based on EMEP/CEIP gridded emission for 2015 [47], where the total BaP emissions for the Czech Republic were estimated to be 8 t based on the Czech 2017 submission. This number was corrected to 16 t in resubmission in 2019. Therefore, the contribution of Czech sources must be underestimated approximately by a factor of two in EMEP results leading to overestimation of relative transboundary contributions by a factor of 1.1-1.6. From the text above it seems likely that modeling only the dispersion of BaP can lead to an overestimation of the transboundary contribution of BaP in the Czech Republic by factor of 1.5.

Discussion
Transport, industry and services combined do not contribute more than 2% to BaP emissions. The main source of BaP emissions in the Czech Republic is overwhelmingly local heating, especially the combustion of solid fuels in older types of boiler constructions with over-fire and under-fire type of burners (Figures 3 and 4). Compared to other European countries, the Czech Republic has the second highest average BaP emission per capita from the residential sector (1.5 g/(person·year)). The highest average emission is in Poland (1.6 g/(person·year)) whereas the EU average is 0.4 g/(person·year) [49].
The main reason for a high share of BaP emissions from local heating in the Czech Republic is specifically the combustion of solid fuel in older type of boilers (over-fire and under-fire). However, these types of boilers are being gradually replaced by low-emission boilers or by other types of heating. In 2018, the share of over-fire and under-fire boilers was estimated to be around 69%. These replacements are being accelerated by the legislative requirements of Act 201/2012 Coll. [50] that stipulates that after September 2022 only low-emission boilers meeting the parameters of at least 3rd boiler class (as defined in the EN 303-5:2012 [51]) can be in operation for solid fuel combustion in households. The replacement of old boilers was supported by a subsidy program between 2015 and 2019, with the aim of replacing up to 100,000 high-emission boilers. However, replacement of the boiler itself is not a guarantee of efficient emission reduction if the boiler is not operated properly in accordance with operating instructions.
The range of measured concentration of BaP was from 0.4 ng·m −3 at a rural regional station to 7.7 ng·m −3 at an industrial station in 2018. The area where BaP concentrations exceeded the target value was 12.6% in 2018 ( Figure 5). The highest annual average concentrations of BaP have long been recorded throughout the entire Ostrava Region. In this region, there is the highest emission load due to a combination of local heating and the largest share of heavy industry in the Czech Republic. The transboundary contribution in the Ostrava region was estimated to be 40-60% but may be in fact somewhat lower due to model limitations. High BaP concentration values due to the effect of local heating systems have also been monitored in Kladno for a long time. Concentrations exceeding the target value for BaP also occur in Central Bohemia and in a number of municipalities. The values of BaP concentration show that the ambient air concentrations of BaP are high in the Czech Republic in general, which is consistent with the model study of Europe [11], where the authors pointed out the ambient air concentrations of BaP to be substantially high in Central and Central Eastern Europe but also in some other European regions.
Relatively low levels of BaP have been recorded in large cities, like in the center of Helsinki [52] and in Porto [28]. In Porto, transport was identified as the main source of PAHs based on diagnostic ratios. In street canyons in Helsinki, the measured concentrations of BaP were at the same level as those in the urban background and clearly lower than those in suburban detached-house areas.
These results indicate that local traffic has only a minor effect on BaP concentrations, compared with the corresponding effect of small-scale combustion. In the Czech Republic, transport is also a minor emission source of BaP (0.8%); nevertheless, levels of BaP concentrations in Prague were higher than in Porto and in street canyons in Helsinki. The higher BaP concentrations in Prague (Table 1) were caused by regional and long-range transport, especially in the center of the city in areas with a high proportion of remote central heating, where the main emission source of BaP is traffic. In suburban Prague, local heating was as important as it was in suburban Helsinki [52].
The lowest average annual concentrations are estimated at places distant from direct exposure to emission sources (natural mountain areas). The lowest measured BaP concentrations, ranging from 0.3 to 0.7 ng·m −3 , have been recorded at the Košetice regional station. Nevertheless, these values are still relatively high (Table 1) and are above the WHO reference value (0.12 ng·m −3 ). These relative high values of BaP concentrations at a regional background station pointed out the important role of regional and long range transport in the Czech Republic. The importance of regional transport is related to the high contribution of coal combustion to BaP emissions (about 30-40%) and the long lifetime of BaP derived from coal combustion in the atmosphere as found in this study [35].
Moreover, no correlation between background BaP concentration at Košetice and emission was found. The same (lowest) value of background BaP concentration in 2008, 2011, 2014 and 2016−2018 supports our conclusion on the combined effects of emission, meteorological and dispersion conditions and transboundary contributions. Based on the CHMI data, the last four years 2015−2018 of the study period are also years with the highest ratio of good meteorological and dispersion conditions and with BaP emissions decreasing after their peak in 2013. The year 2014 can be characterized by a milder winter season and the shortest heating season for many years [22]. In contrast, the beginning of the study period could be characterized by the greater presence of moderately poor and poor dispersion conditions [53] and the lowest BaP emissions in the study period (Table 1).
Since the higher BaP concentrations are a problem in urban and suburban areas due to local heating, we chose five urban and suburban sites in the Czech Republic for the current concentration level and trend assessment. The sixth rural regional site represents the background ambient air concentration in the Czech Republic.
The highest annual average value in urban and suburban sites was detected in 2008, the lowest annual average value of 1.8 ng·m −3 was detected in 2016, 2017 and 2018 with the widest range, from 0.7 to 3.7 ng·m −3 in 2017. Nevertheless, the lack of correlation of BaP concentration in the five urban and suburban sites and emissions for the whole Czech Republic (mainly from local heating, Figure 8) highlights the possible domination of local influences and the influence of meteorological and dispersion conditions.
We found a significant decreasing trend for average concentrations and winter average concentrations. The development of concentrations from selected stations in the Czech Republic between 2008 and 2014 is comparable to the development of PAHs and BaP concentrations presented in the study of [28] who analyzed data from two suburban sites in Porto between 2004 and 2014. A significant decreasing trend in the framework of their study was also found. A downward trend for all types of stations and at two thirds of total stations over the period 2007−2014 was also presented by EEA [54]. A significant decreasing trend was found at 22% of European rural and urban stations [54].
The decrease in concentrations in the Czech Republic is especially noticeable since 2014 highlighting the influence of milder winter seasons in 2014 and 2018, the prevailing good dispersion conditions in 2015−2018 and decreasing coal consumption in the last few years. The significant decrease at two particular stations (Ostrava-Poruba and Kladno-Švermov) that are among those with the highest concentrations in the Czech Republic (and of course those above target value BaP concentration for the whole study period) point also to the effect of improvements in local heating.
Similarly to other studies [26,[55][56][57], the typical seasonal variation for the BaP concentration has been shown. The BaP concentrations for October-March are more than eight times higher than for April-September. For instance, Albuquerque et al. (2016) [28] found a December-January/June-August ratio of 5 for PAHs for the eleven year study in Porto. The reasons for this are generally known-i.e., seasonal sources as local heating, higher emissions from motor vehicles and less mixing in the atmosphere due to inversions [28,57,58]. During the warmer season, on the other hand, concentrations decrease due to unstable atmospheric conditions favorable towards dispersion, increased chemical and photochemical decomposition of PAHs due to higher levels of solar radiation and higher temperatures and of course also due to decreased emissions from anthropogenic sources [59][60][61].
High values of daily BaP concentrations in winter months associated with local household heating were also recorded during the 2017-2018 campaign measurements in the small settlements of Bochovice, Cerníny and Hřivice (Figure 7), where concentrations of BaP are not regularly monitored and where solid fuel heating predominates. The range of measured BaP concentrations was 0.1-8.0 ng·m −3 iň Cerníny, 0.2-9.8 ng·m −3 in Bochovice, and 1.0-13.6 ng·m −3 in Hřivice. In particular, measured BaP concentrations in the small settlement of Hřivice were as high as and on some days even higher than in Kladno-Švermov, where some of the highest concentrations of BaP in the Czech Republic were recorded. Every year, concentrations there reach high values and exceed the target value by almost four times. Such high levels of BaP concentrations are caused on the one hand by an extremely high density of buildings, which leads to a higher BaP emission density near the surface, and on the other hand by the fact that the town is located in a shallow valley, which leads to a reduction in dispersion of pollutants during cold days. The limited amount of winter season data obtained by the campaign measurements does not allow for a calculation of annual average concentrations. The target limit value of 1 ng·m −3 established by European legislation was exceeded in Bochovice,Černíny and Hřivice on 59%, 54% and almost 100% of measurement days, respectively. The measurements in these three small settlements with low-density population clearly indicates that emissions of BaP by local heating influence the short-term BaP concentration in the surroundings. Local meteorological conditions, orography of the populated area and regional and transboundary long-range transport are further factors which influence the ambient air concentration of BaP. In the Czech Republic, due to the rugged terrain, a number of settlements are located in valleys, where there may be a frequent deterioration of dispersion conditions and thus an increase in pollutant concentrations.
The reference level established for BaP by the WHO of 0.12 ng·m −3 was exceeded at all monitoring sites each year. The target value of 1 ng·m −3 is set by a European directive (EU, 2004) with the aim of avoiding, preventing, or reducing harmful effects on human health and/or the environment as a whole. During the period of 2012-2018, 35-58% of the urban population in the Czech Republic was found to be exposed to BaP concentrations exceeding the above mentioned target value. The lowest number of inhabitants living in the above-target areas was estimated to be in 2018. The highest number of inhabitants (58%) exposed to above-target concentrations of BaP was estimated for the years 2012 and 2017. On average, only 7% of the population lived in the areas with the lowest concentration of BaP between 2012 and 2018. BaP is carcinogenic to humans and has been considered a good indicator for the assessment of risk to human health associated with exposure from PAHs found in the environment. The individual carcinogenic potencies of PAH in relation to BaP can be expressed through the BaP equivalent concentrations (BaP eq.) and evaluation of BaP alone will probably underestimate the carcinogenic potential of the PAHs mixtures [28,62]. The uncertainty of the map is a result of the inadequate number of measurements at rural regional stations and the absence of more extensive measurements in smaller settlements in the Czech Republic, where the air pollution by BaP would demonstrate the fundamental effect of local heating units. In addition, the maps are prepared with a resolution of 1 × 1 km and therefore cannot take into account the local fragmentation of the terrain, which in the case of settlements located in valleys affects the levels of pollutants [63]. Thus, assessment of the interannual changes in the territory affected and population exposed to above-limit concentrations of BaP will also be accompanied by a greater margin of error.

Conclusions
A complex analysis of air pollution from BaP in the Czech Republic was carried out. Ambient air BaP concentrations and their long-term trends were studied to assess the level of BaP in the Czech Republic. The calculated emissions of BaP and modeling of the transboundary contribution to the annual mean concentrations of BaP were quantified to present the causes of BaP air pollution.
The measured concentrations of BaP are high due to high emission load from the combination of local heating and heavy industry in the Czech Republic. The residential sector is responsible for more than 98.8% of BaP emissions.
Many people (50% on average) in the Czech Republic live in the area where the BaP concentrations are above the target value set by a European directive with the aim of avoiding, preventing, or reducing harmful effects on human health and/or the environment as a whole.
On the basis of the observations in small settlements described above, where BaP concentrations are not regularly monitored and where solid fuel heating predominates, it can be assumed that in small settlements, carcinogenic BaP levels may reach high levels in the short-term. Above-target values where BaP is not routinely measured can also be expected in similar municipalities with a high proportion of local heating using solid fuels. Consequently, the fraction of people living in the above-target value areas will be higher than presented.
Due to the high number of small municipalities and particularly due to the high costs for laboratory analyses and limited capacity of the laboratories, the number of measurement locations will always be limited, and therefore it is desirable to specify in more detail data on emissions and to provide substantial support to modeling.
The transboundary contribution to annual average concentration was estimated to be between 46% and 70% for most (80%) of the territory of the Czech Republic. The contribution of Czech sources can exceed 80% in residential heating hot spots. Results are nevertheless subject to limitations of the modeling approach used-it is likely that the transboundary contribution to BaP concentrations in the Czech Republic was overestimated by a factor of 1.5. Apart from residential heating, two other categories of Czech sources with relative contribution to annual local average concentration exceeding 10% were identified: road transport (large cities and vicinity of major roads) and coke oven plants in the Ostrava agglomeration.
The typical seasonal variation for the BaP concentrations has been shown. The BaP concentrations for selected air quality monitoring stations for October-March were more than eight times higher than for April-September. This is in line with the composition of emission sources in the Czech Republic, with the dominance of the local heating sector and with the different influence of meteorological and dispersion conditions in the colder and warmer parts of the year connected with the atmospheric stability and chemical properties of BaP.
We found significant decreasing trends for average concentrations and winter average concentrations. No correlation between background BaP concentration and emission was found. BaP concentration development in the Czech Republic has been influenced by the combined effect of total emission, meteorological and dispersion conditions and transboundary contributions.
The significant decrease at two particular stations belonging to those with the highest concentrations in the Czech Republic also point out the effect of improvements in local heating infrastructure. To conclude, even assuming generally good dispersion conditions and milder winter seasons in future, a significant reduction of BaP emissions is needed to reach the target value for BaP in the Czech Republic.

Acknowledgments:
The authors want to thank all colleagues from CHMI who participate in the measurement and processing of air quality data. We would also like to thank Alexey Gusev and Victor Shatalov from EMEP/MSC-E for providing BaP modeling results and for commenting on the text.

Conflicts of Interest:
The authors declare that there is no conflict of interest regarding the publication of this paper.