Long-Range Transport of Sulfur Dioxide Emissions From External Sources to Tehran

Large-scale emissions of sulfur dioxide (SO 2 ) from the combustion of heavy fuel oils are deteriorating the air quality in Tehran and regularly causing complex atmospheric pollution situations and human health concerns. Our analysis of the long-term SO 2 emission data in Tehran conrmed that the magnitude of local SO 2 emission sources is not adequate to reach SO 2 concentrations to their present levels. Tehran is predominantly affected by regional transport of SO 2 from exterior sources further away located in Iraq, Saudi Arabia, and adjacent provinces neighboring Tehran. Approximately 80% of total SO 2 emissions in Tehran were observed to have impacts from the external hotspots outside of Tehran. While local emission sources only contribute around 20% of the total SO 2 emissions. Bivariate polar plots, k-mean cluster, pairwise polar correlation, and PSCF analysis provided evidence for the impact of large-scale transport of SO 2 emissions from external locations from the west/northwest, north/northeast, and south/southwestern areas of the region. Further observations of these hotspot areas observed in our analysis with TROPOMI satellite data conrmed signicant SO 2 emissions resulting from the consumption of heavy fuel oils in thermal power plants and oil/gas reneries. Overall, the results suggested that the regulatory strategies for controlling local trac emissions of SO 2 in Tehran would not be benecial for reducing public health exposures to SO 2 in Tehran. Such improvements can be attained mainly by diminishing the emission sources located further away from Tehran.


Introduction
High sulfur dioxide (SO 2 ) emission from burning heavy fuels oils in power plants and industries in Tehran and surrounding regions are deteriorating the air quality in that middle eastern megacity in recent years (Amini et al., 2014;Hassanzadeh et al., 2009;Nazari et al., 2012). SO 2 is a highly reactive gas with the largest emission sources from fossil fuel combustion in power plants and industrial processes (Nazari et al., 2012). Smaller shares of SO 2 emissions would correspond to roadway tra c sources and high sulfur content fuels in industries and non-tra c transportations (Nazari et al., 2012). According to the latest emission inventories of Tehran, SO 2 emissions are estimated at around 13,000 tons per year, with 94% of its contributions being related to stationary sources, while mobile sources only contribute to 6% of total SO 2 emission (Shahbazi et al., 2019). The reports of fuel quality from Tehran Air Quality Control Company (AQCC) showed that the sulfur content of diesel fuel consumed by mobile and roadway sources decreased signi cantly by more than 98% since the end of 2016 (Sulfur content decreased from 4000 ppm to around 75 ppm) causing substantial improvements in the production and emissions of SO 2 from mobile sources (Ghadiri et al., 2017;Shahbazi et al., 2019). However, despite the reduction in SO 2 emissions of mobile sources, large amounts of SO 2 were released into the atmosphere, mainly from uncontrolled stationary sources in Tehran and its surrounding regions. Recent studies by NASA OMI satellites capturing more than 500 major point sources of SO 2 emissions across the globe con rmed Iran as one of the greatest anthropogenic SO 2 emission hotspots in the world (Dahiya and Myllyvirta, 2019; Zhong et al., 2020). The increased use of coal combustion and heavy fuel oil for increasing energy production demands in the area, partly due to a lack of an effective ue gas desulfurization (FGD) process and slow implementation of emission standards, raises concern about complex air quality situations in the region. (Nazari et al., 2012).
In some extreme cases, unpleasant odors were experienced in large areas of Tehran, which were later understood to be linked to high SO 2 production and emission transportation in Tehran. Our analysis during the pollution episodes in Tehran con rmed extremely high peaks of SO 2 concentrations lasting for a relatively long period were observed in large areas of the city, indicating the persistence of strong sources of SO 2 in the city and the surrounding areas.
To date, emission regulations and monitoring analysis were far too weak to su ciently evaluate and quantify SO 2 productions and their high emitting hotspot sources in the area. Source apportionment analysis in Tehran showed that sulfate could take up to 40% of total PM 2.5 contributions in the city, most of which were due to high primary SO 2 emission sources (Esmaeilirad et al., 2020). SO 2 has remained a critical air quality concern in Tehran; therefore, a comprehensive investigation of the extent of local and regional transport of SO 2, including their hotspot emission sources, is crucial for policymaking and developing future air pollution mitigation plans in Tehran.
To the best of our knowledge, comprehensive analyses on SO 2 emission and their local and regional hotspot sources are still poorly constrained to be su ciently comprehensive and quanti able. Chemical transport models cannot provide a reliable estimate of SO 2 concentration due to the lack of a proper emission inventory in the study area. In addition, accessing such a long-distance SO 2 emission database was impossible for the regions outside Tehran. Therefore, remote sensing observations and the spatial distributions of columnar SO 2 over Iran and its surrounding area were used for the whole study period.
Furthermore, a series of source apportionment analyses were conducted, and SO 2 pro les were closely evaluated at various regions in Tehran to quantify long-term regional-scale contributions using air quality and long-term SO 2 monitoring networks. Our analysis will clarify how long-range regional transport of SO 2 and the different local emission sources and spatial-scale contributions will govern SO 2 pollution in the region.

Study area
Located in an urban region including 15 million inhabitants in its larger municipal areas, Tehran is considered Iran largest and one of the most populous cities in the world (Bayat et al., 2019;Yeganeh et al., 2021). Tehran's climate can feature a cold and semi-arid characteristic with a yearly mean temperature of 17°C and annual precipitation of about 230 mm. The average annual highs will be 25°C, and the annual lows averages at around 10°C. Extreme temperatures reaching up to 40°C could also be seen during July, while minimum temperatures of up to -10°C could also happen during January.
Comprehensive information on the time series of meteorological data is provided in Fig. S1 in the supporting information.

SO 2 monitoring
Tehran has a network of 35 air quality monitoring stations mainly administered by the Department of Environment (DOE) and the Tehran AQCC. The stations belonging to AQCC were equipped with regular gas and PM analyzers (e.g., CO, NO, NO 2, NO x , SO 2 , O 3 , PM 2.5, and PM 10 analyzers) operating continuously over the city (Yeganeh et al., 2021). SO 2 measurements data conducted at several governmental centers and community districts within Tehran were designated to represent diverse environments in Tehran and were categorized as urban-residential (TARB area), tra c hotspots (SHAD and R21 areas), industrial (R20 area), and residential (AGD and PIR areas). UV uorescence analyzers (Model Serinus 50, Ecotech, Australia and Model AF22, Environment S.A, France) were used to analyse SO 2 in all of the monitoring networks in Tehran. The SO 2 data at six representative areas at various administrative divisions in Tehran were selected. The available database includes SO 2 emission data for the period 2020-2021. The location of the monitoring sites and detailed descriptions and codes can be found in Fig. S2 and Table S1 in the supporting information. The overall data quality for the study period in all of our monitoring areas was reported good. We obtained an average of approximately 96% of data coverage with less than 4% of data missing values. More information on the data quality of this study is provided in Table S2 in the supporting information. Maintenance and calibration of the analyzers are performed regularly according to BS EN 14212: 2012 standard.

Bivariate Polar Plots (BPPs) and k-mean cluster analysis
Bivariate polar plots (BPPs) could be created by apportioning wind speed, wind direction, and measurement data in different wind speed-direction intervals (Khuzestani et al., 2017). The u and v wind component variables will be obtained by deconvoluting the wind components according to Where u is the hourly average wind speed, and θ is the wind direction in degrees. A Generalized Additive Model (GAM) framework has nally been applied to u, v, and concentration (C) data to generate a continuous surface to de ne the measurement levels as a function of the wind constituents to obtain valid source descriptions. GAMs can be expressed as shown in Eq. 2.
Eq. 2 C i = β 0 + ∑ n j = 1 S j X ij + e i Where C i is the pollutant i concentration, β0 is the average of response sj (xij) is the smooth function of value i of covariate j, n is the entire amount of covariates, and ei is the residual i. For better model ( ) ( ) √ ( ) diagnostics and improved residual distributions, the concentration data Ci will be square root transformed. Analysis of a wide range of scenarios suggests that appropriate details of surface concentration distributions could be obtained by binning 10-degree intervals of wind direction and 30 wind speed intervals. This strategy could guarantee an active data reduction approach without sacri cing the reliability of the polar plot. Due to the signi cant variability in wind direction data, it is usually recommended to use this data from various months or years to create a more nely resolved and robust surface interpolation of BPPs (Carslaw et al., 2006;Westmoreland et al., 2007).
To further analyze the BPP's identi ed features regarding the potential source characteristics, we performed a k-mean cluster analysis. This strategy can obtain similar cluster features in BPPs based on the similarity or dissimilarity between the data points. Polar clusters can select interesting features in BPPs as clusters based on different wind components. The fundamentals of polar clusters are k-mean clustering the BPPs as follows: Initially, random k points representing the initial group centroids will be selected from the surface generated by BPPs. All the objects and points will then be allocated to the groups that have the closest centroid. The positions of the k-points will then be repeatedly recalculated until the centroids no longer move. k-mean cluster analysis is similar to the idea of clustering trajectories described in detail in Carslaw and Beevers (2013). However, wind components and concentration data rather than the air mass origins will be used to identify similar groups.

TROPOspheric Monitoring Instrument (TROPOMI) data analysis
TROPOMI was launched on October 13, 2017, by the Sentinel-5 Precursor (S5P) from northern Russia. The satellite orbits the Earth at 817 km and has a repetition cycle of 17 days (KNMI, 2017). The SP5 has four spectrometers, one covering the short-term infrared and the other three covering the ultraviolet-near infrared ranges. The satellite data processing could obtain high-resolution (3.5 x 7 km 2 ) daily coverage data of O 3 , NO 2 , SO 2 , CO, CH 4 , HCHO, cloud and aerosol properties. Detailed information on the methodology of SO 2 data processing has been described elsewhere in KNMI, 2017. We analyzed the spatial distributions of columnar SO 2 over Iran and its surrounding region for the whole year of 2020.
In summary, The OFFL/L3-SO 2 offers high-resolution o ine images of SO 2 concentration levels.
With 1-day revisit time and the high spatial resolution (3 x 5.5 km), S5P/TROPOMI is able to provide su cient details, including the recognition of minor SO 2 plumes. For the L3-SO 2 product, the absorbing aerosol index is computed with measurements at 340 and 380 nm wavelengths. The OFFL/L3 products for the ares within the product's bounding box were constructed using the harp convert tool. The areaaverage data for the pixels with different values for different times were then merged together into a greater mosaic containing sets of different tiles with orthorecti ed raster data.

Potential Source Contribution Function (PSCF)
The hybrid single-particle Lagrangian integrated trajectory model (HYSPLIT version 5) was implemented to calculate the hourly-interval 48-h backward trajectories arriving at heights of 500 m above ground level at the monitoring areas in Tehran using a vertical velocity model. The Global Data Assimilation System (GDAS) one-degree archive data combined measurements into a gridded space with observed data (Khuzestani et al., 2017). Fig. S3 in the supporting information illustrates the frequencies of the trajectories arriving in Tehran, highlighting that most of the back trajectories originated from the west of the study area. The Potential Source Contribution Function (PSCF) analysis could determine probable geographical locations and hotspot emission sources by calculating the probability elds resulting in high concentrations at the monitoring area. Detailed information about the PSCF analysis has been described elsewhere in Khuzestani et al. (2017). PSCF values representative of the potential hotspot areas could be calculated as follows: Where n ij is the total number of endpoints that pass through the grid cell (i, j), and m ij is the number of endpoints related to the values that exceed the threshold criterion values in the same grid cell.

SO 2 variation pro les
The general daily variation pro le of SO 2 in Tehran is displayed in Fig. 1a. Almost all monitoring areas showed two peaks for SO 2 concentration when the rst peak coincides in the morning rush hours and the during the morning rush hours in addition to the HDDVs in parallel with reducing the PBL height, respectively. Considering that most SO 2 in Tehran was generated from stationary sources, we can expect that the operation of thermal power plants and the available industries in the region would also largely contribute to higher SO 2 concentrations in the area. In addition to the stationary sources, mobile sources would also contribute to the SO 2 pro les in Tehran.
Diurnal variation patterns of SO 2 were considerably different across the sampling sites in Tehran. Fig. 1d illustrates the time series of the diurnal variations of the concentrations of SO 2 across all the sampling sites in Tehran. SO 2 concentrations at R21 and SHAD areas increased at around 6:00 local time and moved to their highest levels at around 9:00 local time. During the afternoon, concentrations decreased again and reached their lowest values at around 17:00 local time. The concentration of SO 2 during the nighttime period started to increase and reached its maximum at around 20:00 local time. R21 and SHAD areas are close to roadway tra c, and therefore their diurnal variation is greatly in uenced by the tra c activities. SO 2 concentration at these stations could also be affected by the northwestern emission sources in Tehran, and therefore, higher concentrations were experienced in these stations.
In the R20 area, we observed maximum SO 2  The weekly variation pro les of SO 2 were also illustrated in Fig. 1c

Bivariate Polar Plot (BPP) analysis
The role of wind components on SO 2 concentration variations for analysis of local and regional potential sources in uences in Tehran are evaluated via the bivariate polar plots (BPPs) and illustrated in Fig. 2a-f. As displayed in Fig. 2, relatively similar pro les were observed in BPPs of all the monitoring areas in Tehran, demonstrating two primary potential hotspot sources for the SO 2 emissions. The rst identi ed hotspot area in the BPPs in all the monitoring areas showed high potential SO 2 sources under low wind speed (<5 m/s). Little wind direction dependence slightly expanded toward the southern directions in Tehran. We expect that these areas mainly dominate local emission sources originating from roadway tra c emissions and domestic heating under stagnant atmospheric conditions that help the potential accumulation and build-up of SO 2 emissions from these local hotspot areas, resulting in higher SO 2 concentrations (Hama et al., 2020).
In addition to the local hotspot areas illustrated in BPPs in all our monitoring areas, high potential SO 2 emission sources were observed in northwestern areas indicating high potential SO 2  We have performed a non-parametric wind regression (NWR) approach of modeling the surfaces in BPPs using kernel smoothers that weigh SO 2 on a surface according to the proximity to de ned wind speed and direction intervals described by Henry et al. (2009). NWR approach provides the means to evaluate the stability of our analysis and con rm the consistency of the long-range transported hotspot areas illustrated in BPPs in all the monitoring areas. This analysis is similar to the BPP analysis described earlier, but it has signi cant advantages when su cient data are not available to generate surfaces in some areas of BPPs. Fig. S4 in the supporting information shows the results of the NWR analysis of all the monitoring areas in Tehran. As illustrated in Fig. S4, our NWR approach shows relatively similar results to our BPP analysis. In addition to the local (central areas in the plot) and regional hotspots (located in the northwestern areas), high potential regional sources under high wind speeds (≈10 m/s) originated from the southern areas were also observed in NWR analysis, illustrating probabilities of existence of high emitting SO 2 sources located further away in southern parts of Tehran. The NWR approach generally showed similar results for our monitoring areas east of Tehran (e.g., AGD, PIR, and R20), pointing to northwestern areas as hotspot sources. However, other mentoring areas located in the western parts of Tehran mainly showed hotspot sources originating from the south/southwestern areas (Fig. S4). This again con rmed the availability of 2 main regionally transported emission sources of SO 2 located further away in the west and southwest of the study area.

k-mean cluster analysis
To further analyze the potential hotspot areas identi ed by BPP and NWR analysis, we performed k-mean clustering of the polar plots allowing various sources of SO 2 in Tehran to be revealed in terms of the clusters. This method can characterize SO 2 features into various homogenous clusters by minimizing the sum of the squared euclidean distances between the observed data and the cluster center. We performed k-mean cluster analysis for 2 -10 cluster numbers, and the results are provided in Fig. S5 in the supporting information. Our analysis and comparison of the clustering database with the results obtained from BPP and NWR analysis showed that the 8-cluster solution could su ciently de ne all the characteristics of SO 2 variation features. The results from our 8-cluster analysis as well as the corresponding NWR and BPP plots for the general SO 2 patterns in Tehran are illustrated in Fig. 3a-c.
According to Fig. 3c, cluster 3 of our 8-cluster solution mostly represents the local sources; however, cluster 2 and cluster 7 could probably correspond to the external distant hotspot locations transported in Tehran from the southern and western areas, respectively, as observed in BPP and NWR plots ( Fig. 3a and  b). All other sources, rather than from clusters 3 and 6, may also correspond to long-range external locations emitting SO 2 in Tehran (will be discussed in the following section). However, our postprocessing of SO 2 pro les in all of the monitoring stations in Tehran (Fig. 4) con rmed that clusters 2, 3, and 7 had the highest contributions compared with other clusters observed in Fig. 3c. As described earlier, cluster 3 mainly corresponded to local SO 2 emissions generated from the roadway tra c emissions; however, clusters 2 and 7 mainly corresponded to distant external hotspot emission sources located further away from the region.
The constructed diurnal variation pro les of cluster 3 (Fig. S6) agreed reasonably well with the variations of general SO 2 diurnal pro les in Tehran observed in Fig. 1, implying that cluster 3 could remarkably illustrate the characteristics of local emission sources in all the monitoring stations in Tehran. The diurnal variation pro les in clusters 2 and 7 were not so profound and did not agree well with that observed for cluster 3 and the general SO 2 pro les implying that other emission sources rather than local tra c emissions could also impact SO 2 pro les in other clusters in Tehran.
Our observations showed that large SO 2

Pairwise polar correlation analysis
Pairwise regression of SO 2 concentrations for all the individual sites was compared to evaluate the impacts from local/regional sources of all SO 2 emissions across all the monitoring areas in Tehran. A linear t was integrated with each regression analysis to de ne the slope, slope uncertainty, intercept, intercept uncertainty, and R 2 value (Table S3). To statistically evaluate the effects of the local/regional sources between the monitoring areas in Tehran, the slopes and intercepts of the regression lines generated from each pair were analyzed to see if they are statistically signi cant by deciding if 1 minus the slope is larger than two times the uncertainty, or the absolute intercept value is larger than two times of intercept uncertainty (Khuzestani et al., 2017;McGinnis et al., 2014). Detailed information on the pairwise regression analysis is provided in Table S2 in the supporting information. Most of our location pairs showed the in uence of local emissions by indicating statistically signi cant slopes and intercepts. However, some of our pairs illustrated the regional transport of SO 2 emissions in Tehran.
To visually illustrate the effects of regional hotspot sources in the region, we consider illustrating the Pearson correlations coe cients (r) derived from the pairwise regression analysis between each pair to vary against wind components. The combination of pairwise correlation analysis with polar plots could provide signi cant insight into the source apportionment information (Grange et al., 2016). This analysis involves weighting r values obtained from pairwise regression analysis by wind speed/direction bins. A Gaussian kernel model used for polar plots was implemented to guarantee that maximum data is used in every wind speed-direction interval (Grange et al., 2016). The results of our pairwise polar correlation analysis are provided in Fig. 5. Very high correlation coe cients between all the pairwise cases were experienced under higher wind speeds and different wind directions. Signi cant sources of SO 2 in Tehran were observed from external emission sources, again illustrating the effects of long-range transport of SO 2 in Tehran. Comparison analysis of the cases depending on the location, characteristics, and the distance between our monitoring areas could reveal different hotspot areas for SO 2 emission sources.
However, regardless of the pairwise correlation analysis, all the cases were able to detect relatively similar source origins corresponding to large long-range transport of SO 2 emissions from external locations in the west/northwestern areas (cluster 4 and 7), north/northeast areas (cluster 5 and 8) and southern/southwest areas (cluster 1 and 2). This agreed perfectly with the results from our k-mean clustering analysis, BPP, NWR, and nally PSCF analysis, con rming that all the clusters (except clusters 3 and 6) could correspond to signi cant regional long-range transport of SO 2 emission from external locations in Tehran. None of our pairwise correlation polar plots detected low wind speed central areas (clusters 3 and 6), implying signi cant impacts of long-range transport of SO 2 from several external high emitting sources located further away from Tehran.

Effects of regional transport
The TROPOMI data analysis illustrated an extensive line of SO 2 emission sources covering western parts of the country with main hotspots related to Iraq, Kuwait and Saudi Arabia. Some other high emitting sources located in Turkey and Uzbekistan were also observed in TROPOMI data analysis (Fig. 6). PSCF analysis results for all the monitoring sites in Tehran are also shown in Fig. 7. Our PSCF analysis showed a high probability of emission sources in the Iraq and Kurdistan area in the west, Saudi Arabia in the southwest, and adjacent parts in the north, east and southeastern areas. Interestingly, our further observations of these hotspot locations showed that the observed source regions agreed perfectly with the results of the TROPOMI SO 2 properties of the spatial distributions of columnar SO 2 over Iran and its surrounding region for the whole year of 2020 (Fig. 6) as well as the NASA global SO 2 emission source catalog as identi ed power plants and oil and gas re neries with high SO 2 emissions in the region. This was also consistent with the results of our BPP, k-mean clustering, and polar correlation analysis, further demonstrating the evidence of large-scale transport of SO 2 from distant hotspots located away from Tehran. This suggests that high SO 2 emissions from the industries in Iran's neighboring countries (i.e., Iraq and Saudi Arabia as identi ed in our PSCF and TROPOMI SO 2 data analysis) in addition to the high emitting SO 2 sources located in surrounding areas of Tehran would signi cantly affect SO 2 emissions in our monitoring sites.

Conclusions
SO 2 emission is considered a critical air quality concern in Tehran, often causing complex atmospheric pollution situations in that middle eastern megacity. Herein, comprehensive and quantitative studies on SO 2 emission production and their local and regional hotspot sources in Tehran were evaluated using various analysis methods, including bivariate polar plots, k-mean clustering, pairwise correlation analysis, potential source contribution function as well as the remote sensing TROPOMI SO 2 spatial distribution.
Our analysis of different monitoring areas in Tehran con rmed that local sources originated from roadway tra c emissions, and domestic heating only contributes a minor fraction (20%) to the total SO 2 productions in Tehran. Approximately 80% of the total SO 2 emissions were transported in Tehran throughout long-range transport from different external hotspot areas located further away from Tehran. Our analysis with BPP, NWR, pairwise polar correlation, and PSCF analysis con rmed that similar source origins corresponding to large long-range transport of SO 2 emissions from external locations related to Iraq and Kurdistan area in the west, Saudi Arabia in the southwest, and adjacent areas in north, east and southeastern areas. Our further observations of those hotspot areas with remote sensing TROPOMI SO 2 data processing illustrated large SO 2 emissions resulting from the consumption of heavy fuel oils in thermal power plants located in the west and southwestern areas of Tehran. Lack of proper control measures, including ue gas desulfurization process to control SO 2 emissions from thermal power plants and oil/gas re neries, has resulted in signi cant SO 2 emissions in Tehran from external hotspot locations transported in the city. The variability in contributions of local/regional emission sources has implications for potential control strategies that could reduce public exposures. Our analysis showed that limiting the high emitting vehicles and other regulatory strategies for controlling tra c-related emissions could be bene cial in reducing the exposure to many pollutants in Tehran. However, such strategies would be less effective for reducing the exposure to SO 2 because of the signi cant contribution of longrange SO 2 transport from external locations affecting the city.
Declarations Figure 1 Hourly, monthly and weekly concentration pro les of SO 2 in Tehran. a, b, and c panels represent concentrations pro les for the average SO 2 concentrations, while d, e, and f panels represent the variation patterns for each speci c monitoring site in Tehran

Supplementary Files
This is a list of supplementary les associated with this preprint. Click to download. SIv8.docx