Sulfur and Nitrogen Depositions in BULGARIA—Model Results and Observations

: Atmospheric deposition processes are of primary importance for human health, forests, agricultural lands, aquatic bodies, and ecosystems. South-East Europe is still characterized by numerous hot spots of elevated sulfur deposition, despite the reduction in European emission sources. The purpose of this study is to discuss the results from two chemical transport models and observations for wet and dry depositions of sulfur (S), reduced nitrogen (RDN) and oxidized nitrogen (OXN) in Bulgaria in 2016–2017. The spatial distribution and the domain main deposition values by EMEP MSC-W (model of the MSC-W Centre of the Co-operative Programme for Monitoring and Evaluation of the Long-range Transmissions of Air Pollutants in Europe) and BgCWFS (Bulgarian Chemical Weather Forecast System) demonstrated S wet depositions to be higher than N depositions, and identiﬁed a rural area in south-east Bulgaria as a possible hot-spot. The chemical analysis of deposition samples at three sites showed a prevalence of sulfate in the western part of the country, and prevalence of Cl and Na at a coastal site. The comparison between modeled and observed depositions demonstrated that both models captured the prevalence of S wet depositions at all sites. Better performance of BgCWFS with an average absolute value of the normalized mean bias of 16% was found. in and underestimated S-DD in Ahtopol, while EMEP-CTM signiﬁcantly underestimated S-DD in Ahtopol. The underestimation of the S-DD in Ahtopol by both models was most probably linked to the resolution of the models, limiting the representation of local circulations in this complex coastal area. Another cause might be related to local emission sources—e.g., intensive road trafﬁc in summer and wood burning from residential heating in autumn—not accounted for in the model for the studied period.


Introduction
The atmospheric deposition has been intensively studied worldwide over the last decades, mainly due to concerns over acid rains, eutrophication, trace metal deposition, damages to forests and vegetation, ecosystem health, and global climate change [1][2][3][4]. The atmospheric deposition is part of complex air pollution and environmental processes that link emissions of pollutants, their chemical transformation and sinks, and their effects on the earth's surface. The naturally occurring process of removal of pollutants compounds (gaseous and aerosols) from the atmosphere is analyzed through two main pathways-wet and dry. Wet deposition refers to the scavenging of dissolved gases and particles through precipitation, while the dry deposition refers to processes of direct downward transport of pollutants from the atmosphere to the earth surface. Historically, wet deposition received more attention due to the negative impacts of acid rains on forests. Studies on total depositions (wet plus dry) were promoted in relation to preservation of ecosystems health, but also in support to estimating the effects of emission reductions.
In this work, we focus on sulfur (S), reduced nitrogen (RDN) and oxidized nitrogen (OXN) wet, dry and total depositions in Bulgaria, a country of South-East Europe. The significant anthropogenic emissions of sulfur dioxide (SO 2 ) come from combustions of fossil fuels, e.g., coal-fired thermal power plants (TPP), oil refineries and industrial facilities. Nitrogen oxides (NO x) originate mainly from the transport (area source), but large point BgCWFS is based on WRF (3.6.1), [35], and CMAQ (4.6), [36] models, recommended by the US Environmental Protection Agency. The system runs operationally over the following 5 nested domains: from European scale (81 km grid resolution), through Balkan Peninsula domain (27 km), Bulgaria domain (9 km), down to Sofia region (3 km) and Sofia city (1 km), [37,38]. BgCWFS was set up in 2016 for routine calculations of wet and dry depositions in Bulgaria. The results analyzed here are for domain Bulgaria.
The meteorological driver WRF was fed by analysis fields from the Global Forecast System of the US National Centers for Environmental Prediction (NCEP GFS) with horizontal resolution of 1 • × 1 • and temporal resolution of 6 h. Initial conditions for CMAQ were part of previous day calculations. A predefined set of vertical concentration profiles was used as chemical boundary condition for the European domain, while all other domains received their boundary conditions from the previous one in the hierarchy. The emissions were prepared using TNO inventory for 2011 [39] and Bulgarian national emission inventory for 2015. Extensively tested parameterizations schemes were selected for the simulations, e.g., in WRF-the Yonsei University scheme for the planetary boundary layer [40], the WSM6 scheme for microphysics [41], the Kain-Fritsch scheme for cumulus parameterization [42], and the Noah land surface model [43]. The chemical mechanism in CMAQ v.4.6 is the 4th generation module "cb4 ae4 aq", dry deposition of particles is calculated by turbulent air motion and by direct gravitational sedimentation of large particles [44], and wet deposition is calculated in the cloud module of CMAQ [45].
Preliminary analysis of BgCWFS wet depositions [25,29] suggested model overestimation due to overestimation of precipitation amounts. Thus, a precipitation bias adjustment (PBA) approach was applied for wet deposition calculations following [46]. The method involves linear correction of model wet depositions by the ratio of the observed to estimated precipitation. In previous studies, [25,29] it was shown that the PBA approach led to improvement of monthly wet depositions in Bulgaria when comparing the model to observed values at single sites. Sensitivity tests on the effect of PBA on monthly wet depositions in Bulgaria suggested a decrease in annual wet depositions by about 25-30%, [31]. In this work, the PBA correction is applied as post-processing to the gridded wet depositions on monthly and daily basis. The method requires gridded values for observed precipitation, [31]. The objective analysis of precipitation is a challenging task due to the heterogeneity in the spatial distribution. We used the precipitation analysis system developed at the Forecast Center at NIMH. It combines observations for accumulated precipitation in the country and gridded data from the analysis field of a weather forecasting model using Cressman analysis, [47]. The forecasted precipitation field is used as a first guess, further corrected for the difference between forecasted and observed precipitation amount at a given point. The weight function depends on the distance between the grid points and the observation point, as well as on the difference in their elevation.
EMEP-CTM results for yearly, monthly and daily values of various species on a grid covering Europe with horizontal resolution 0.1 • x 0.1 • were extracted from [34]. We used data for dry and wet depositions of S, RDN and OXN from model version rv4.33 with input for meteorology and emissions for the year of simulation (in our case 2016 and 2017). The meteorological driver for EMEP-CTM is the ECMWF-IFS model, the Integrated Forecast System model (IFS) of the European Centre for Medium-Range Weather Forecasts (ECMWF) (https://www.ecmwf.int/en/research/modeling-and-prediction, accessed on 10 February 2021). The updates of EMEP-CTM (rv.4.33) included the following significant changes with respect to previous model versions: new gas-phase chemical mechanism, carbon-bond scheme, schemes for secondary organic aerosols, new methods for emissions input [48]. For this work, the EMEP-CTM data for the territory of Bulgaria was extracted and interpolated to the BgCWFS grid with 9 km resolution.
Both models differ in their parameterizing schemes, chemistry mechanisms, initial and boundary conditions, and emissions input. As emissions play a key role in air quality modeling, we discuss in the next Section the differences between the two modeling systems relative to emissions.

Emissions
The emissions input for the two modeling systems was for different years, as follows: EMEP-CTM used emissions for 2016 and 2017 [49], while BgCWFS used emissions based on the national emission inventory for 2015, and the TNO inventory for 2011 for the rest of Europe. Table 1 shows officially reported total emissions for two domains-country (Bulgaria) and European-wide (27 countries: EU-27). The years in Table 1 correspond to the years used for emissions in the two models. Thus, the emissions with labels "BG 2015" and "EU-27 2011" correspond to BgCWFS, while the emissions with labels "BG 2016 and 2017" and "EU-27 2016 and 2017" correspond to EMEP-CTM. The data were extracted from the emissions database of EMEP [50].
For Bulgaria, the SO x emissions used in EMEP-CTM (for the years 2016-2017) were about 27% lower than the emissions in BgCWFS (for the year 2015). The significant reduction in national SO x emissions from 2015 to 2016 was due to sources from the public power sector, mainly coal-fired TPPs. The differences for the national total emissions of NOx and NH 3 in the years mentioned are much smaller. Thus, the differences between BgCWFS and EMEP-CTM for NOx and NH 3  The distribution of the emission sources over the territory of the country (spatial allocation) is also important in air quality modeling. There were some notable differences, mainly for SO x source allocations, between the two models. The locations of the main SO x large point sources in Bulgaria (TPP and industrial facilities) are shown in Figure 1. The bars in Figure 1 indicate SO x emissions (kilotons per year) as used by the two models and as reported in the European Pollutant Release and Transfer Register (EPRTR) [51] for 2017. While BgCWFS overestimated the emissions from TPPs with site numbers 1-3, the EMEP-CTM model underestimated these emissions and overestimated largely the emissions at TPP site number 4 ( Figure 1a). For the industrial facilities ( Figure 1b) the main overestimation by EMEP-CTM compared to reported SO x emissions was for the plant for copper production and processing of metals "Aurubis" (site number 7), and for the plant KCM producing non-ferrous metals (site number 9). The SO x emissions for the facilities with site numbers 1 to 10 are listed in Table 2. The last column in Table 2 includes satellite derived SO 2 emissions for the TPPs with site numbers 1 and 3 in 2017. They were extracted from the global catalogue of large SO 2 emission sources based on the onboard Ozone Monitoring Instrument (OMI) [52,53]. It can be noticed that the EMEP-CTM model underestimated by a factor of 6 the emissions from the TPP complex "Maritsa East" (site number 1, south-east Bulgaria), which is one of the biggest coal-fired power stations on the Balkan Peninsula. For Bulgaria, the SOx emissions used in EMEP-CTM (for the years 2016-2017) were about 27% lower than the emissions in BgCWFS (for the year 2015). The significant reduction in national SOx emissions from 2015 to 2016 was due to sources from the public power sector, mainly coal-fired TPPs. The differences for the national total emissions of NOx and NH3 in the years mentioned are much smaller. Thus, the differences between BgCWFS and EMEP-CTM for NOx and NH3 emissions in Bulgaria were small. For the EU-27 domain, the SOx emissions in BgCWFS (for the year 2011) were about 41% higher than the SOx emissions in EMEP-CTM (for the years 2016-2017). Thus, BgCWFS used higher SOx emissions than EMEP-CTM, both for the national and the European-wide domains. Another difference between EU-27 emissions and BG emissions was in the trend between sulfur and nitrogen emissions. While in EU-27 the SOx emissions in 2016-2017 were much lower than the emissions of NOx and NH3, in Bulgaria SOx emissions were comparable to NOx emissions.
The distribution of the emission sources over the territory of the country (spatial allocation) is also important in air quality modeling. There were some notable differences, mainly for SOx source allocations, between the two models. The locations of the main SOx large point sources in Bulgaria (TPP and industrial facilities) are shown in Figure 1. The bars in Figure 1 indicate SOx emissions (kilotons per year) as used by the two models and as reported in the European Pollutant Release and Transfer Register (EPRTR) [51] for 2017. While BgCWFS overestimated the emissions from TPPs with site numbers 1-3, the EMEP-CTM model underestimated these emissions and overestimated largely the emissions at TPP site number 4 ( Figure 1a). For the industrial facilities ( Figure 1b) the main overestimation by EMEP-CTM compared to reported SOx emissions was for the plant for copper production and processing of metals "Aurubis" (site number 7), and for the plant KCM producing non-ferrous metals (site number 9). The SOx emissions for the facilities with site numbers 1 to 10 are listed in Table 2. The last column in Table 2 includes satellite derived SO2 emissions for the TPPs with site numbers 1 and 3 in 2017. They were extracted from the global catalogue of large SO2 emission sources based on the onboard Ozone Monitoring Instrument (OMI) [52,53]. It can be noticed that the EMEP-CTM model underestimated by a factor of 6 the emissions from the TPP complex "Maritsa East" (site number 1, south-east Bulgaria), which is one of the biggest coal-fired power stations on the Balkan Peninsula.   Table 2. Table 2. SO x emissions (kt.y-1) from main point sources in Bulgaria from different sources: as used in the models (BgCWFS, EMEP), as provided by the EPRTR register [51], and as estimated based on retrievals of the satellite on-board OMI instrument [53]. TPP denotes coal-fired thermal power plants.    Table 2. SOx emissions (kt.y-1) from main point sources in Bulgaria from different sources: as used in the models (BgCWFS, EMEP), as provided by the EPRTR register [51], and as estimated based on retrievals of the satellite on-board OMI instrument [53]. TPP denotes coal-fired thermal power plants.   Wet and dry depositions in Sofia and Ahtopol were sampled with an automatic wet only device (WADOS, Kroneis GmBH). The sampler contains a humidity sensor which controls the lid of wet and dry collector compartments automatically. During the wet deposition events, the sensor moves the lid onto the dry collector and after the sensor surface becomes dry, the lid on the dry collector goes onto the wet collector. Wet deposition samples were sampled on daily base (used for monthly estimates), while dry deposition samples were collected on a monthly basis in Sofia and Ahtopol. Deionized water was added into the bucket for dry deposition before sampling. A passive bulk sampler was operated at Cherni Vrah. After each sampling, the funnel was washed with deionized water. Wet and dry depositions in Sofia and Ahtopol were sampled with an automatic wet only device (WADOS, Kroneis GmBH). The sampler contains a humidity sensor which controls the lid of wet and dry collector compartments automatically. During the wet deposition events, the sensor moves the lid onto the dry collector and after the sensor surface becomes dry, the lid on the dry collector goes onto the wet collector. Wet deposition samples were sampled on daily base (used for monthly estimates), while dry deposition samples were collected on a monthly basis in Sofia and Ahtopol. Deionized water was added into the bucket for dry deposition before sampling. A passive bulk sampler was operated at Cherni Vrah. After each sampling, the funnel was washed with deionized water.
In order to calculate the wet deposition of sulfur and nitrogen species, the measured concentration (mg·L −1 ) of each species (SO 4 2-, NO 3 − , NH 4 + ) was multiplied by the observed precipitation amount (L·m −2 ). The concentration of non-sea salt (non-marine) sulfate (nss_SO 4 2− ) was estimated by a correction based on the assumption that sodium is a sea salt tracer: [nss_SO 4 2− ] = [SO 4 2− ] − (0.25 × [Na]), following WMO recommendations [54]. In this work, data for the period from June to December 2017 were used.

Evaluation of Model Performance
The wet depositions as estimated by the two modeling systems were analyzed in two steps: model to model, and models to observations. The models were inter-compared on a quarterly and annual basis for the years 2016 and 2017, focusing on domain mean values (the average of the depositions for the grid nodes located in Bulgaria). The spatial distribution of S, RDN, and OXN wet and dry depositions was visualized on maps and compared qualitatively. The performance of the two models for precipitation was discussed in terms of annual domain mean and maximum values, quarterly domain mean values, and qualitative comparison of annual precipitation maps.
The comparison of modeled to observed depositions at the three sites was carried out for mean monthly values for the period from June to December 2017. Daily mean wet depositions at Cherni Vrah were analyzed only for three short-term periods, characterized by particular meteorological situations. The performance of the models for the selected periods was commented in view of mean values and normalized mean bias (NMB), where negative values indicate underestimation and positive values overestimation by the models.

Precipitation
Wet deposition estimates depend strongly on the precipitation amounts. Thus, we comment firstly differences between the modeled precipitations. Figure 3 shows the spatial distribution of the annual accumulated precipitation in 2017 as simulated by BgCWFS with application of the precipitation bias adjustment (PBA) and by EMEP-CTM, and as constructed using observations. Both models provided higher precipitation amounts for the mountain areas. To note also that the models indicated elevated precipitation in the most south-eastern part of the country where the sampling site Ahtopol is located.
In order to calculate the wet deposition of sulfur and nitrogen species, the measured concentration (mg.L −1 ) of each species (SO4 2-, NO3 − , NH4 + ) was multiplied by the observed precipitation amount (L.m −2 ). The concentration of non-sea salt (non-marine) sulfate (nss_SO4 2− ) was estimated by a correction based on the assumption that sodium is a sea salt tracer: [nss_SO4 2− ] = [SO4 2− ] − (0.25 × [Na]), following WMO recommendations [54]. In this work, data for the period from June to December 2017 were used.

Evaluation of Model Performance
The wet depositions as estimated by the two modeling systems were analyzed in two steps: model to model, and models to observations. The models were inter-compared on a quarterly and annual basis for the years 2016 and 2017, focusing on domain mean values (the average of the depositions for the grid nodes located in Bulgaria). The spatial distribution of S, RDN, and OXN wet and dry depositions was visualized on maps and compared qualitatively. The performance of the two models for precipitation was discussed in terms of annual domain mean and maximum values, quarterly domain mean values, and qualitative comparison of annual precipitation maps.
The comparison of modeled to observed depositions at the three sites was carried out for mean monthly values for the period from June to December 2017. Daily mean wet depositions at Cherni Vrah were analyzed only for three short-term periods, characterized by particular meteorological situations. The performance of the models for the selected periods was commented in view of mean values and normalized mean bias (NMB), where negative values indicate underestimation and positive values overestimation by the models.

Precipitation
Wet deposition estimates depend strongly on the precipitation amounts. Thus, we comment firstly differences between the modeled precipitations. Figure 3 shows the spatial distribution of the annual accumulated precipitation in 2017 as simulated by BgCWFS with application of the precipitation bias adjustment (PBA) and by EMEP-CTM, and as constructed using observations. Both models provided higher precipitation amounts for the mountain areas. To note also that the models indicated elevated precipitation in the most south-eastern part of the country where the sampling site Ahtopol is located. The domain mean annual precipitation, as estimated by observations, was higher in 2017 than in 2016 (Table 3). This was correctly simulated by BgCWFS, while EMEP-CTM provided almost identical domain mean precipitation amounts. BgCWFS overestimated The domain mean annual precipitation, as estimated by observations, was higher in 2017 than in 2016 (Table 3). This was correctly simulated by BgCWFS, while EMEP-CTM provided almost identical domain mean precipitation amounts. BgCWFS overestimated the observed domain mean precipitation in 2017 by about 8%, while EMEP-CTM underestimated it by about 5%. A significant difference between the models was noted for the domain-wide maximum annual precipitation, as follows: BgCWFS results were with NMB of about −4%, while EMEP-CTM results were with NMB of about −29%. The differences between the models on a quarterly basis (Table 4) indicate that BgCWFS simulated more precipitation than EMEP-CTM in the period from April to June (AMJ) and less precipitation in the period from July to September (JAS). Compared to observed precipitations, EMEP-CTM performed better with NMB from −2% to 9%, while the NMB for BgCWFS was from −23% to 25%.  observed  234  254  112  137  162  232  167  248  BgCWFS  204  296  86  150  203  270  150  254  EMEP  246  263  122  142  150  228  155  236 3.

Domain Mean Wet Depositions
Annual and quarterly S, OXN and RDN wet depositions (WD) simulated by the two models for 2016 and 2017 were compared as domain mean values, Table 5. Both models indicated prevalence of sulfur over nitrogen wet depositions. Annual values for S-WD by BgCWFS were higher by a factor of 1.5 than the ones by EMEP-CTM, most likely linked to the higher emissions used in BgCWFS. The annual values of RDN-WD and OXN-WD by BgCWFS were, however, lower than the EMEP-CTM values by about 12% and 28%, respectively. Both models simulated higher S-WD for the period from January to June, in correspondence to higher emissions and precipitations in this period. For the last three months of 2017 the S-WD by both models were relatively high, most likely due to the higher precipitation amounts by the models. For nitrogen depositions, the most significant difference was in the period from June to September, when EMEP-CTM values for RDN and OXN-WD were twice as high as the respective depositions by BgCWFS.

Spatial Distribution of Wet Depositions
The spatial distribution of annual wet depositions in 2017 showed significant differences between the two models ( Figure 4). In general, BgCWFS wet depositions were higher in regions with more precipitation (mountain areas), while EMEP-CTM wet depositions were higher near the emission sources used in this model. The maps for S-WD clearly indicated the difference in the locations of significant SO x emission sources: on the BgCWFS map this is TPP complex "Maritsa East", while on the EMEP-CTM map this is the industrial plant "Aurubis". As discussed in Section 2.1.2, EMEP-CTM had short-comings in the SO x emission quantities for these two facilities-the fingerprint of the first one is invisible on the EMEP-CTM maps, while the second one appears as a hot-spot. The maps for the differences between BgCWFS and EMEP-CTM wet depositions (Figure 4c) showed that S-WD and RDN-WD by BgCWFS were higher than the ones by EMEP-CTM over the whole country, while for OXN-WD the opposite was true. Bigger urban areas were noticed as hot spots of OXN-WD on the EMEP-CTM map, while such urban hot-spots were missing on the BgCWFS map. This suggested a deficit in NO 2 emissions input in BgCWFS. The spatial distribution of annual wet depositions in 2017 showed significant differences between the two models ( Figure 4). In general, BgCWFS wet depositions were higher in regions with more precipitation (mountain areas), while EMEP-CTM wet depositions were higher near the emission sources used in this model. The maps for S-WD clearly indicated the difference in the locations of significant SOx emission sources: on the BgCWFS map this is TPP complex "Maritsa East", while on the EMEP-CTM map this is the industrial plant "Aurubis". As discussed in Section 2.1.2, EMEP-CTM had short-comings in the SOx emission quantities for these two facilities-the fingerprint of the first one is invisible on the EMEP-CTM maps, while the second one appears as a hot-spot. The maps for the differences between BgCWFS and EMEP-CTM wet depositions (Figure 4c) showed that S-WD and RDN-WD by BgCWFS were higher than the ones by EMEP-CTM over the whole country, while for OXN-WD the opposite was true. Bigger urban areas were noticed as hot spots of OXN-WD on the EMEP-CTM map, while such urban hotspots were missing on the BgCWFS map. This suggested a deficit in NO2 emissions input in BgCWFS.
For the most south-eastern part of the country, characterized by numerous protected areas of nature, both models indicated higher values of S-WD and OXN-WD.
The maps for 2016 (not shown here) had similar patterns.

Spatial Distribution of Dry Depositions
The spatial distribution of modeled dry depositions depends mainly on the location of the sources. As both models exploited different emission inventories, significant differences were noted in the corresponding maps ( Figure 5). BgCWFS represented better the position of the national source with highest SO x emissions-the TPP complex "Maritsa East" in the southern part of the country (marked by circle in Figure 5). The map for the differences in S-DD indicated that BgCWFS provided much higher values than EMEP-CTM over almost the whole country. For OXN-DD, both models simulated higher values in the eastern part of the country, where coastal big cities were better evident on the EMEP-CTM map. The main differences in the spatial pattern for RDN-DD were in the northern part of the country, where BgCWFS values were higher by a factor of two.

Total Depositions of Sulfur, Reduced Nitrogen and Oxidized Nitrogen
The total (wet plus dry) depositions for domain mean annual values estimated by BgCWFS and EMEP-CTM models are presented in Figure 6. Sulfur total depositions were higher than the reduced nitrogen and the total oxidized nitrogen depositions for both models. The most significant difference between the models was for sulfur deposition, e.g., the domain mean total S deposition by BgCWFS was higher than the one by EMEP-CTM by a factor of about 3.5. For reduced nitrogen depositions, both models provided similar values, while for total oxidized nitrogen deposition the difference showed lower BgCWFS values by about 28%. A difference between the models was noted when comparing sulfur total depositions with total nitrogen depositions (reduced plus oxidized). BgCWFS indicated a notable prevalence of sulfur over nitrogen depositions, while EMEP-CTM suggested a prevalence of nitrogen deposition over sulfur deposition, by about 21% in 2016 and 14% in 2017. EMEP-CTM map. The main differences in the spatial pattern for RDN-DD were in the northern part of the country, where BgCWFS values were higher by a factor of two.

Total Depositions of Sulfur, Reduced Nitrogen and Oxidized Nitrogen
The total (wet plus dry) depositions for domain mean annual values estimated by BgCWFS and EMEP-CTM models are presented in Figure 6. Sulfur total depositions were higher than the reduced nitrogen and the total oxidized nitrogen depositions for both models. The most significant difference between the models was for sulfur deposition, e.g., the domain mean total S deposition by BgCWFS was higher than the one by EMEP-CTM by a factor of about 3.5. For reduced nitrogen depositions, both models provided similar values, while for total oxidized nitrogen deposition the difference showed lower BgCWFS values by about 28%. A difference between the models was noted when comparing sulfur total depositions with total nitrogen depositions (reduced plus oxidized). BgCWFS indicated a notable prevalence of sulfur over nitrogen depositions, while EMEP-CTM suggested a prevalence of nitrogen deposition over sulfur deposition, by about 21% in 2016 and 14% in 2017.

Chemical Analysis of Precipitation Samples
The percentage contribution of ions and elements to the total mass of all analyzed species in precipitation samples for the period from June to December 2017 are presented in Figure 7. The mean concentrations of all elements were generally higher for the mountain site. An exception was noticed for the coastal site Ahtopol where concentrations of Cl and Na were higher than for the other sites. Generally, nss_SO4 2-was found to be the dominant anion in precipitation samples from Sofia and Cherni Vrah (39.6% and 31.2%, respectively), followed by NO3 − (21.6% and 25.6%, respectively) and Cl − (4.9% and 10.4%, respectively). At the coastal site Ahtopol, the Cl anion showed a higher contribution (41%), as was expected for a site near the sea, followed by nss_SO4 2-(12.3%) and NO3 − (11.1%). The contribution of sea/marine aerosol to the total concentration of sulfates was also high for samples from Ahtopol. On average for Ahtopol, the sea salt sulfate contribution in wet deposition samples was 22%, and in dry deposition samples it was 18%. These values were much higher than for Sofia and Cherni Vrah, where in wet deposition samples the sea salt sulfate contribution was only about 1% of the total sulfate concentrations.

Chemical Analysis of Precipitation Samples
The percentage contribution of ions and elements to the total mass of all analyzed species in precipitation samples for the period from June to December 2017 are presented in Figure 7. The mean concentrations of all elements were generally higher for the mountain site. An exception was noticed for the coastal site Ahtopol where concentrations of Cl and Na were higher than for the other sites. Generally, nss_SO 4 2was found to be the dominant anion in precipitation samples from Sofia and Cherni Vrah (39.6% and 31.2%, respectively), followed by NO 3 − (21.6% and 25.6%, respectively) and Cl − (4.9% and 10.4%, respectively). At the coastal site Ahtopol, the Cl anion showed a higher contribution (41%), as was expected for a site near the sea, followed by nss_SO 4 2-(12.3%) and NO 3 − (11.1%). The contribution of sea/marine aerosol to the total concentration of sulfates was also high for samples from Ahtopol. On average for Ahtopol, the sea salt sulfate contribution in wet deposition samples was 22%, and in dry deposition samples it was 18%. These values were much higher than for Sofia and Cherni Vrah, where in wet deposition samples the sea salt sulfate contribution was only about 1% of the total sulfate concentrations. tain site. An exception was noticed for the coastal site Ahtopol where concentrations of Cl and Na were higher than for the other sites. Generally, nss_SO4 2-was found to be the dominant anion in precipitation samples from Sofia and Cherni Vrah (39.6% and 31.2%, respectively), followed by NO3 − (21.6% and 25.6%, respectively) and Cl − (4.9% and 10.4%, respectively). At the coastal site Ahtopol, the Cl anion showed a higher contribution (41%), as was expected for a site near the sea, followed by nss_SO4 2-(12.3%) and NO3 − (11.1%). The contribution of sea/marine aerosol to the total concentration of sulfates was also high for samples from Ahtopol. On average for Ahtopol, the sea salt sulfate contribution in wet deposition samples was 22%, and in dry deposition samples it was 18%. These values were much higher than for Sofia and Cherni Vrah, where in wet deposition samples the sea salt sulfate contribution was only about 1% of the total sulfate concentrations. The predominant cation for Sofia and Cherni Vrah was Ca representing, with a share of, respectively, 13.4% and 17.4%, while for Ahtopol the predominant cation was Na (19.9%). The share of ammonium ions was highest for Sofia (7.12%), for Cherni Vrah it was 4.75%, and for Ahtopol it was 2.9%.

Comparison to Wet Depositions at Three Sites
The comparison for S, RDN and OXN wet depositions was carried out on monthly basis for the period from June to December 2017 at all three stations. Hereafter, the observed S-deposition is corrected for sea salt contribution. Both models correctly simulated high S-WD at the three sites ( Figure 8). However, EMEP-CTM significantly underestimated S-WD at two of the sites. The absolute value of the mean NMB was about 15% for BgCWFS (with overestimation for Cherni Vrah), and about 43% for EMEP-CTM (with underestimation at all sites). In another study, based on comparison of EMEP-CTM to data from the EMEP monitoring network in 2017 [55], an underestimation of sulfur depositions was also reported, to the value of about 27% on average. The analysis in [55] further showed EMEP-CTM model results for OXN-WD without significant bias with respect to observations, while the RDN-WD were overestimated by 17%. For the period investigated in this work, the EMEP-CTM wet depositions at the three stations were characterized also by underestimation for S-WD (varying from 28% to 54%), an overestimation for RDN-WD (varying from 59% to 98%), and underestimation for OXN-WD, (varying from 51% to 69%). These higher biases with respect to the ones reported in [55] were partly due to the shorter period of analysis in this study.
observations, while the RDN-WD were overestimated by 17%. For the period investigated in this work, the EMEP-CTM wet depositions at the three stations were characterized also by underestimation for S-WD (varying from 28% to 54%), an overestimation for RDN-WD (varying from 59% to 98%), and underestimation for OXN-WD, (varying from 51% to 69%). These higher biases with respect to the ones reported in [55] were partly due to the shorter period of analysis in this study.  Figure 8 indicates that BgCWFS performs, in general, better than EMEP-CTM at the sampling sites, despite the usage of outdated emissions. The accumulated precipitation amount by EMEP-CTM is higher than the observed one for the studied period (Figure 8d). This suggests that the better agreement of BgCWFS results to observed wet depositions might be linked to better allocation of the main national emission sources. The significant underestimation of S-WD at stations Sofia and Cherni Vrah by EMEP-CTM is most probably due to underestimation of emissions from the large emission sources in the western part of the country, as shown in Figure 1.  Figure 8 indicates that BgCWFS performs, in general, better than EMEP-CTM at the sampling sites, despite the usage of outdated emissions. The accumulated precipitation amount by EMEP-CTM is higher than the observed one for the studied period (Figure 8d). This suggests that the better agreement of BgCWFS results to observed wet depositions might be linked to better allocation of the main national emission sources. The significant underestimation of S-WD at stations Sofia and Cherni Vrah by EMEP-CTM is most probably due to underestimation of emissions from the large emission sources in the western part of the country, as shown in Figure 1.

Comparison to Dry Depositions at Two Sites
The percentage contribution of ions and elements to the total mass of all analyzed species in dry deposition samples for the periods from June to September, and from October to December 2017 is presented in Figure 9.

Comparison to Dry Depositions at Two Sites
The percentage contribution of ions and elements to the total mass of all analyzed species in dry deposition samples for the periods from June to September, and from October to December 2017 is presented in Figure 9. For Sofia, the contributions of NO3 − and nss_SO4 2− were higher in the colder period (October to December) with a share of 30.7% and 27%, respectively, in comparison to 14.4% and 6.4%, respectively, in the warmer period (June to September). The contributions of Cl, Na, K, NH4 + were also higher in the colder period (10.6%, 5.4%, 30.5%, 4%, respectively) than the shares in the warmer period (3%, 1.6%, 6.7%, 2.6%, respectively,). For Ahtopol, the contributions of Cl and Na were higher in the warmer period (30.3% and 9%, For Sofia, the contributions of NO 3 − and nss_SO 4 2− were higher in the colder period (October to December) with a share of 30.7% and 27%, respectively, in comparison to 14.4% and 6.4%, respectively, in the warmer period (June to September). The contributions of Cl, Na, K, NH 4 + were also higher in the colder period (10.6%, 5.4%, 30.5%, 4%, respectively) than the shares in the warmer period (3%, 1.6%, 6.7%, 2.6%, respectively,). For Ahtopol, the contributions of Cl and Na were higher in the warmer period (30.3% and 9%, respectively) than in the colder period (18.4% and 7%, respectively). The nss_SO 4 2− contribution in the warmer period (11.5%) was lower than for the colder period (22.4%). The contribution of NO 3 − (3%) in Ahtopol for both periods was much smaller than in Sofia. For the whole period from June to December 2017, the dominant anion for Sofia was NO 3 − (18%), followed by nss_SO 4 2-(12.5%) and Cl − (9.8%). For Ahtopol, the dominant anion was Cl − (24.8%), followed by nss_SO 4 2-(16.8%) and NO 3 − (3%). There was a notable share of the crustal elements Ca and K in the samples from Sofia. For both sites, the predominant cation in dry deposition samples was Ca, with a share of 23.6% and 33.6%, respectively. The share of ammonium ions was higher for Sofia (3.4%) than for Ahtopol (1.1%).
The percentage contribution of different elements in the dry deposition samples for Sofia and Ahtopol for the period from June to December 2017 was in the following order: Sofia: Ca > K > NO 3 − > nss_SO 4 2-> Cl > Na > Mg > NH 4 + > Si > Fe > Zn > Cu, Ahtopol:Ca > Cl > nss_SO 4 2-> Na > Mg > K > Si > NO 3 − > NH 4 + > Zn > Fe > Cu. The comparison between model and observed dry depositions accumulated in the period June-December 2017 ( Figure 10) shows that BgCWFS largely overestimated S-DD in Sofia, and underestimated S-DD in Ahtopol, while EMEP-CTM significantly underestimated S-DD in Ahtopol. The underestimation of the S-DD in Ahtopol by both models was most probably linked to the resolution of the models, limiting the representation of local circulations in this complex coastal area. Another cause might be related to local emission sources-e.g., intensive road traffic in summer and wood burning from residential heating in autumn-not accounted for in the model for the studied period. The observed S-DD at the coastal site Ahtopol were higher than in Sofia. One of the reasons for this was the difference in the amount of accumulated precipitations for the studied period, with values at Ahtopol lower by more than 20% compared to values in Sofia. The number of days with precipitation in Ahtopol was also less (16) than in Sofia (25), suggesting that the removal of pollutants in Sofia was mainly through the wet deposition. Higher dry depositions of sulfur in Ahtopol in summer months (June to August) could be linked also to influence from the industrial plant "Lukoil Neftochim" Burgas, the largest refinery in the Balkans, located at about 60 km northwards. The SO2 concentration observed in Burgas at the background urban station BG0056 "Meden Rudnik" (not far away from the refinery plant "Lukoil Neftochim") had a mean value of 8.90 µg·m −3 for the The observed S-DD at the coastal site Ahtopol were higher than in Sofia. One of the reasons for this was the difference in the amount of accumulated precipitations for the studied period, with values at Ahtopol lower by more than 20% compared to values in Sofia. The number of days with precipitation in Ahtopol was also less (16) than in Sofia (25), suggesting that the removal of pollutants in Sofia was mainly through the wet deposition. Higher dry depositions of sulfur in Ahtopol in summer months (June to August) could be linked also to influence from the industrial plant "Lukoil Neftochim" Burgas, the largest refinery in the Balkans, located at about 60 km northwards. The SO 2 concentration observed in Burgas at the background urban station BG0056 "Meden Rudnik" (not far away from the refinery plant "Lukoil Neftochim") had a mean value of 8.90 µg·m −3 for the studied period.
At the same time, the mean SO 2 concentration in Sofia, observed at the background urban site BG0079 Sofia "Mladost" in the vicinity of the sampling site, was significantly lower (3.12 µg·m −3 ). The prevailing synoptic flow along the Bulgarian Black Sea coast might have contributed to higher dry depositions in Ahtopol. During the months of July and August, the prevailing winds at the southern Bulgarian Black Sea coast are from N to NE, the so called Etesian winds [56]. These rather intense winds contribute to higher turbulence and enhance the dry deposition of pollutants. The second half of the period (September to December) was characterized by more frequent passages of atmospheric perturbations from NW and SW, also reaching the coastal area. Thus, higher dry depositions in Ahtopol for the period June to December could be linked to influence from both national emission sources and sources outside of the country.

Effects of Long Range Transport on Depositions at Mountain Site Cherni Vrah
A combined analysis of back-trajectories and data from the chemical analysis of precipitation samples collected at the high mountain site Cherni Vrah was carried out aiming to identify the effects of long-range transport on precipitation chemical composition.
We briefly discuss three periods of 2017 characterized by elevated values of sulfate and nitrate concentrations in the samples: 16-17 June, 12-13 August and 26-28 September. Back-trajectories were calculated at three heights above ground level (500 m, 1000 m, and 2000 m) using the model HYSPLIT [57,58].
The wet depositions of sulfur, reduced nitrogen and oxidized nitrogen were estimated from the two models, BgCWFS and EMEP-CTM, and also based on the chemical analysis of the wet samples (observations). Then, the NMB was calculated for each model and each period. Figure 11 shows the pathway of the air masses and the NMB for BgCWFS and EMEP-CTM for the three cases. In the first period, the synoptic situation was characterized by the passage of a fast cold front from NW with light to moderate precipitations in the western part of the country. North-westerly outbreaks are typical for the western part of the country, where both Cherni Vrah and Sofia are located. In such types of outbreaks, the wet depositions of sulfur prevail over wet depositions of nitrogen [29]. For the two days 16-17 June, the non-sea salt sulfur deposition was 36.1 mg·m −2 , and the total nitrogen deposition was 31.2 mg·m −2 , as estimated based on observations. Both models simulated correctly the prevalence of S-WD over N-WD, but underestimated S-WD and RDN-WD. The NMB for depositions by BgCWFS was lower (Figure 11a), likely due to lower bias in precipitation (−35% compared to −70% by EMEP-CTM). The pathway suggests influence from TPPs located north-west of the country. The second period is typical for August in Bulgaria when anticyclonic weather conditions prevail, and the precipitations are often associated with local phenomena or the passage of low-pressure systems southward of the country. In this case, both models overestimated the precipitation by a factor of 2-3, but the NMB by the models for different types of deposition was variable (Figure 11b), suggesting higher impact by emission sources. The chemical composition showed relatively high amount of K (0.46 mg·L −1 ), likely due to influence from forest fires in the western part of the Balkans. Satellite images (https://firms.modaps.eosdis.nasa.gov/, accessed on 10 February 2021) indicated numerous forest fires south-west of the country in the period 09-13 August. The third period was characterized also by easterly winds (Figure 11c). The chemical analysis of the precipitation sample showed elevated concentrations of Cu (on average 0.014 mg·L −1 ), most likely due to influence from the Aurubis copper smelter plant located eastward of the sampling site (site number 7 in Figure 1). The precipitation was underestimated by both models; however, the NMB for the depositions were of opposite signs by the two models, with a better score for EMEP-CTM. 09-13 August. The third period was characterized also by easterly winds (Figure 11c). The chemical analysis of the precipitation sample showed elevated concentrations of Cu (on average 0.014 mg·L −1 ), most likely due to influence from the Aurubis copper smelter plant located eastward of the sampling site (site number 7 in Figure 1). The precipitation was underestimated by both models; however, the NMB for the depositions were of opposite signs by the two models, with a better score for EMEP-CTM.

Discussion
Results from two modeling systems (BgCWFS, based on WRF-CMAQ, and EMEP-MSC-W rv4.33) were compared for wet and dry depositions of sulfur, reduced nitrogen and oxidized nitrogen for the territory of Bulgaria in 2016 and 2017. The models have similar grid resolutions (about 10 km), but different inputs and parameterization schemes.

Discussion
Results from two modeling systems (BgCWFS, based on WRF-CMAQ, and EMEP-MSC-W rv4.33) were compared for wet and dry depositions of sulfur, reduced nitrogen and oxidized nitrogen for the territory of Bulgaria in 2016 and 2017. The models have similar grid resolutions (about 10 km), but different inputs and parameterization schemes. Significant differences were noted in the emission input, as follows: in BgCWFS, the emissions were for the year 2015 at national level, and for 2011 for the European domain, while EMEP-CTM exploited up-to-date emissions for the years 2016 and 2017. The analysis of SO x emissions used by the two models and their comparison to data from the European Pollutant Release and Transfer Register for the year 2017 showed short-comings in EMEP-CTM relevant to the locations and intensity of the main point sources in Bulgaria. For example, annual emissions from the coal-fired thermal power plant complex "Maritsa-East" in south Bulgaria were very low with respect to the officially reported data. This complex remains one of the major SO x polluters in the Balkans, despite the reduction in emissions in recent years, and it is important to be accurately represented in gridded emissions. EMEP-CTM undergoes annual validation based on observations from the EMEP monitoring sites for precipitation chemistry, which are not available in Bulgaria. This might have contributed to the undetected problems with emissions allocation in the country. The differences in the emission input led to different values of depositions, both in magnitude and in terms of spatial distribution. The outdated emissions in BgCWFS resulted in higher depositions especially for sulfur dry and wet depositions. Domain-averaged annual total depositions of sulfur estimated by BgCWFS were, by a factor of 3.5, higher than the ones simulated by EMEP-CTM. Another significant difference was noted for the total deposition of oxidized nitrogen, in that the domain-averaged value by BgCWFS was 28% lower than the one by EMEP-CTM. Only for reduced nitrogen wet depositions did both models provide similar domain-averaged values. The maps for the spatial distributions revealed significant differences in the footprint of the biggest emitters of pollutants at national scale (e.g., TPP complex "Maritsa East").
Despite the differences in the results from the two systems, the model methodologies revealed some common features. Firstly, sulfur deposition in Bulgaria in 2016 and 2017 was higher than the depositions of reduced nitrogen and of oxidized nitrogen. Secondly, the most south-eastern part of the country was subject to both elevated S-WD and wet and dry OXN depositions. This could be related to effects of transboundary pollution. According to the country report by EMEP for 2017, the transboundary contribution to oxidized sulfur deposition and oxidized nitrogen deposition is about 80% for the south-eastern part of Bulgaria [26]. As this region has numerous protected areas of nature, further studies are necessary to investigate the depositions' variability and the possible causes, taking into account effects both from national sources and from transboundary pollutant transport.
Prevailing sulfur depositions (wet plus dry) over the Balkans were also reported in the literature. In [24], based on model simulations for the period 2000-2007, the multi-year averaged S depositions in winter were more than 40% higher than the S deposition in summer. In our study, BgCWFS results averaged for 2016-2017 showed total S depositions for the period from January to March as being about 33% higher than for the period from June to September. This variability is linked to the intensity of the emissions sources in the region, where numerous coal-fired TPPs are operating. Additionally, in [15], the Balkans were identified as hot spots for S total depositions based on modeling results from 14 different CTMs for 2010. The trend of decreasing sulfur depositions in Europe in the period 1990-2010 was reported to be more pronounced before 2000, as a consequence of emission reduction strategies in Europe [14]. SO x emissions in Europe decreased by more than 80% during 2000-2019, but the observations in this period indicated an average reduction of 60% for sulfur wet deposition and that there were still many hot spots in South-East Europe [5]. The spatial distribution of depositions in Bulgaria for the period (2008-2014), discussed in [59], showed some similarities to the results in this work. In [59], the same CTMs were applied, but with different inputs. Similar to our results, the highest differences between the models in [59] was for sulfur deposition, both as magnitude and spatial distribution, indicating the importance of the SO x emission input to the models.
The chemical analysis of wet and dry atmospheric samples collected during June-December 2017 at three sites in Bulgaria showed prevalence of sulfate in the western part of the country, and Cl and Na at the coastal site Ahtopol. The sea salt contribution was also high at this site, with a share of 22% in wet precipitation samples and 18% in dry samples. The observed sulfur depositions (wet and dry) at all sites was higher than the oxidized and the reduced nitrogen depositions. Interestingly, the observed sulfur dry deposition at the coastal site was higher than the one in Sofia. Most likely this is due to the combined effect of local/regional circulations leading to more windy conditions, and influence from national and transboundary emission sources. A high share of sulfate and chlorine in wet samples at the south-eastern Bulgarian coast was also found in a previous study [30]. For the summer and autumn period of 2014, the predominant elements in wet samples from Ahtopol and Burgas were chlorine and sulphate, 39% and 19.5%, respectively. The analysis in [30] pointed out the effects of cyclone paths on the chemical composition of precipitation samples, showing that non-sea salt sulfate contribution increases for cyclones approaching the south-eastern part of Bulgaria from southern directions. The reported data for the observed dry deposition in Ahtopol in this work was in line with another study for the same site [60]. For July and August of 2018, the observed S dry deposition at Ahtopol was also higher than OXN-DD and RDN-DD, and the share of sea salt sulfate was 16% on average for the two months [60].
The collection of atmospheric samples at the three sites in Bulgaria (urban, mountain and coastal) does not have regular character, as they are carried out during experimental campaigns. The chemical analysis of these samples and the estimation of deposition, however, provide valuable background for the validation of modeling results in this part of South-East Europe, with a lack of precipitation chemistry data and dry deposition data.
The comparison of modeled to observed depositions showed that both BgCWFS and EMEP-CTM captured the prevalence of sulfur wet depositions in the period from June to December 2017 at all sites. S wet deposition by BgCWFS was with lower NMB (15%) compared to NMB by EMEP-CTM (58%) (NMB in absolute values). One possible reason for this better score of BgCWFS is the correct representation of the locations of national emission sources. The sulfur dry deposition at the coastal site Ahtopol was largely underestimated by both models (on average by more than 50%). In [60], results for S-DD at Ahtopol for the months from June to August of 2018 pointed to underestimation by EMEP-CTM of about 43%. Additional analysis with finer scale models suitable to represent the complex local circulations at these sites might be needed to have a more in depth insight for the depositions in this part of the country.
Variability in model results for depositions was reported by European wide modelintercomparison exercises [14,15]. Differences among the models are not surprising as they use different inputs and process parameterizations, and even a single model could provide different results due to different model versions. The differences between BgCWFS and EMEP-CTM depositions with reference to observations at three sites in Bulgaria were within the range reported in [14,15]. For S-WD, the NMB was between −51% and +38%, lower than the reported range in [14] (−70% to +82%). For RDN-WD, the NMB was in the range from −17% to +119%, compared to the one in [14] (−58%, +21%). For OXN-WD, the NMW was in the range from −67% to +32% within the range reported one in [14] (−86%, +72%). The results presented here for dry depositions showed a difference between model results by a factor of 2-3 for S-DD and OXN-DD and more agreement for the modeled RDN-DD, in line with the results discussed in [14].

Conclusions
We presented results for depositions in Bulgaria as estimated by two models: BgCWFS (the operational chemical weather forecasting system in Bulgaria) and EMEP-MSC-W rv3.33 (the model used for annual reporting on transboundary fluxes of particulate matter, photooxidants, acidifying and eutrophying components in Europe). The horizontal resolution of both models was similar (about 10 km), but emissions input, meteorological drivers and parameterization schemes were different. As a consequence, the spatial distribution of sulfur, reduced nitrogen and oxidized nitrogen had different patterns. However, some common features were noted as follows: both systems estimated sulfur depositions to be prevailing, with hot spots corresponding to big emission sources (coal-fired TPP or industrial facilities). Both models also indicated higher sulfur and oxidized nitrogen deposition in the most south-eastern part of the country-a region without significant emission sources and with numerous protected areas of nature. The variability of the modeled depositions was in agreement with findings from other studies for the Balkans and from well-known model-intercomparison exercises in Europe. The comparison of the model depositions to observed depositions at three sites in the country pointed out to some shortcomings in the models. For the region of the Balkans, where routine measurements of precipitation chemistry or dry deposition are missing, the data discussed here, though limited in time and space, provide the possibility to test model performance and thus contribute to the improvement of the models' applications.