Multi-source SO 2 emission retrievals and consistency of satellite and surface measurements with reported emissions

Reported sulfur dioxide (SO2) emissions from US and Canadian sources have declined dramatically since the 1990s as a result of emission control measures. Observations from the Ozone Monitoring Instrument (OMI) on NASA’s Aura satellite and ground-based in situ measurements are examined to verify whether the observed changes from SO2 abundance measurements are quantitatively consistent with the reported changes in emissions. To make this connection, a new method to link SO2 emissions and satellite SO2 measurements was developed. The method is based on fitting satellite SO2 vertical column densities (VCDs) to a set of functions of OMI pixel coordinates and wind speeds, where each function represents a statistical model of a plume from a single point source. The concept is first demonstrated using sources in North America and then applied to Europe. The correlation coefficient between OMI-measured VCDs (with a local bias removed) and SO2 VCDs derived here using reported emissions for 1 by 1 gridded data is 0.91 and the best-fit line has a slope near unity, confirming a very good agreement between observed SO2 VCDs and reported emissions. Having demonstrated their consistency, seasonal and annual mean SO2 VCD distributions are calculated, based on reported point-source emissions for the period 1980–2015, as would have been seen by OMI. This consistency is further substantiated as the emission-derived VCDs also show a high correlation with annual mean SO2 surface concentrations at 50 regional monitoring stations.


Introduction
Sulfur dioxide (SO2) is a designated criteria air pollutant that enters the atmosphere through anthropogenic (e.g., combustion of sulfur-containing fuels, oil refining processes, metal ore smelting operations) and natural processes (e.g., volcanic eruptions and degassing). Over the past three decades both the US and Canada have taken measures to reduce atmospheric emissions of SO2 in order to combat acidification of the ecosystem (e.g., acid rain) and fine particulate matter. As a result, 5 between 1990 and 2012, reported emissions of SO2 declined by 78 percent in the United States and 58 percent in Canada (IJC, 2014). In this study, we examined how well the changes in the reported emissions agree with the SO2 changes in North America observed by satellite and surface instruments.
Ground-based networks such as the US Clean Air Status and Trends Network (CASTNet) and Canadian Air and Precipitation Monitoring Network (CAPMoN) are specifically designed to monitor long-term trends of gaseous pollutants in 10 rural areas away from major pollution emission sources (Baumgardner et al., 1999;Park et al., 2004;Schwede et al., 2011).
Their measurements show that over the eastern US, reductions in regional SO2 emissions have led to significant reductions in monitored SO2 concentrations (Sickles II and Shadwick, 2015;Xing et al., 2013).
Satellites provide global measurements of SO2 vertical column densities (VCD): the total number of molecules or total mass per unit area (Krotkov et al., 2008;Li et al., 2013;Theys et al., 2015). They have been previously used to study 15 the evolution of SO2 VCDs over large regions such as Europe , China (Jiang et al., 2012;Koukouli et al., 2016;Li et al., 2010;Witte et al., 2009), India , and the U.S. (Fioletov et al., 2011). Satellite instruments can detect anthropogenic SO2 signals from large individual point sources such as copper and nickel smelters, power plants, oil and gas refineries , and other sources (Bauduin et al., 2014(Bauduin et al., , 2016Carn et al., 2004Carn et al., , 2007Fioletov et al., 2013;de Foy et al., 2009;Lee et al., 2009;McLinden et al., 2012McLinden et al., , 2014Nowlan et al., 2011;Thomas et al., 2005). An 11-year-long record 20 of satellite SO2 data over different regions of the globe, including the eastern US and southeastern Canada, was examined recently . The analysis shows a substantial (up to 80%) decline in the observed VCD values over that region.
These satellite measurements can also be used as an independent source to verify reported changes in emissions.
Methods for emission estimates from satellite measurements have been recently reviewed by . One such 25 method that does not require the use of atmospheric chemistry models has been commonly used in recent years. By first merging observations from the Ozone Monitoring Instrument (OMI) with wind information, the downwind decay of several pollutants can be analyzed, and in so doing estimates of the total SO2 (or NO2) mass (α) near the source and its lifetime or, more accurately, decay time (τ) can be derived (Fioletov et al., 2011(Fioletov et al., , 2015de Foy et al., 2015;Lu et al., 2013Wang et al., 2015). The emission strength (E) can be obtained using the expression E=α/τ if we assume a steady state for these 30 quantities. The mass can be derived directly from satellite measurements, while the lifetime can either be prescribed using known emissions (Fioletov et al., 2013) or estimated from the measurements based on the rate of decay of VCD with distance downwind (Beirle et al., 2014;Carn et al., 2013;de Foy et al., 2015). Model-based comparisons of different 3 methods to estimate E and τ demonstrate that such methods can produce accurate estimates of τ (de Foy et al., 2014). In our previous study (Fioletov et al., 2015), values of α and τ for anthropogenic point sources were derived from OMI measurements by fitting a 3-dimensional function of the geographic coordinates and wind speed.
These methods, however, are applicable to individual point sources. When this condition is not met, as is the case for multiple sources, the sources can either be combined together if they are close (Fioletov et al., 2015), or the fitting domain is 5 split and the sources are fit separately . Both approaches have their limitations. In this study, we derive a general relationship between emissions and VCDs that can be used for the estimation of emissions from multiple sources.
Moreover, the approach can be used in reverse: that is, VCDs can be estimated directly from reported emissions data, thus making it possible to study the link between VCDs and surface concentrations even for the period before satellite measurements became available. This study is focused on the eastern U.S. and southeastern Canada, where the majority of 10 large North American SO2 emissions sources (mainly coal-burning power plants) are located, where the changes in both reported emissions and measured VCDs are particularly large, and where emissions are measured directly at the stack for most sources. In this region, there is also a network with long-term records of uniform SO2 surface concentration measurements. All of this makes it possible to study consistency between the measurements of emissions, VCDs, and surface concentrations. Once the link between these measurements is verified, it is possible to estimate one measured quantity from 15 another. As an illustration, we demonstrate how European SO2 emissions can be estimated from OMI VCD data.

Satellite SO2 VCD data
OMI, a Dutch-Finnish UV-Visible wide field of view nadir-viewing spectrometer flying on NASA's Aura spacecraft (Schoeberl et al., 2006), provides daily global coverage at high spatial resolution (Levelt et al., 2006). OMI has the highest 20 spatial resolution and is the most sensitive to SO2 sources among the satellite instruments of its class (Fioletov et al., 2013).
OMI SO2 VCD data are retrieved for 60 cross-track positions (or rows). In order to use only data with the highest 25 spatial resolution, we excluded data from the first 10 and last 10 cross-track positions from the analysis to limit the acrosstrack pixel width from 24 km to about 40 km, while the along-track pixel length was about 15 km (de Graaf et al., 2016). In other words, a single OMI measurement represents an SO2 VCD value averaged over a 350-500 km 2 area.
Measurements with snow on the ground were excluded from the analysis as the OMI PCA algorithm presently does not account for the effects of snow albedo. Only clear-sky data, defined as having a cloud radiance fraction (across each 30 pixel) less than 20%, and only measurements taken at solar zenith angles less than 70° were used. Beginning in 2007, up to a half of all rows were affected by field-of-view blockage and stray light (the so-called "row anomaly") and those affected 4 pixels were also excluded. Additional information on the OMI PCA SO2 product can be found in other publications McLinden et al., 2015). SO2 VCD data from the Ozone Mapping Profiler Suite (OMPS) Nadir Mapper on board the Suomi National Polarorbiting Partnership (or Suomi NPP) satellite operated by NASA/NOAA and launched in October 2011 were also used in the study to verify a potential bias in some OMI data (see the Supplement, section S1). OMPS data were processed with the 5 same PCA algorithm as OMI data (Li et al., 2013;Zhang et al., 2017). OMPS has a lower spatial resolution than OMI, 50 km by 50 km, but better signal-to-noise characteristics.
To eliminate cases of transient volcanic SO2, periods when SO2 values observed over the eastern U.S. were affected by volcanic emissions; we determined and excluded such cases from the analysis. The range of analyzed SO2 VCD values was limited to a maximum of 3 DU. Since the average SO2 value over the largest SO2 source in the US is about 1 DU and the 10 standard deviation of individual measurements is 0.5 DU, such a limit corresponds to the 4 standard deviations level even over even the largest sources. Of the SO2 values over the eastern U.S. and southern Canada considered here, the years 2008 and 2009 are particularly problematic due to the eruptions of Kasatochi (Aleutian Islands, Alaska, August 2008, 52N) and Sarychev (Kuril Islands, Eastern Russia, June 2009, 48N). High volcanic SO2 values were also observed on several days in 2011. In addition to the filtering based on SO2 values, five time intervals were explicitly removed from the analysis to avoid 15 misinterpretation of volcanic SO2 as anthropogenic pollution. The intervals are: 07.07. 2008-23.07.2008, 08.08.2008-08.09.2008, 23.03.2009-10.04.2009, 16.06.2009-05.07.2009, and 22.05.2011-09.06.2011. To remove volcanic SO2 in the case of Europe, the analyzed data were divided into 5° by 5° cells, and for each cell, days with the 90 th percentile above a 5 DU limit were excluded from the analysis. Only about 1.5% of all data were removed by this screening.

Wind data 20
As in several previous studies (Fioletov et al., 2015;McLinden et al., 2016), wind speed and direction data for each satellite pixel were required for the analysis methods applied. European Centre for Medium-Range Weather Forecasts (ECMWF) reanalysis data (Dee et al., 2011) (http://apps.ecmwf.int/datasets/) were merged with OMI measurements. Wind profiles are available every 6 hours on a 0.75° horizontal grid and are interpolated in time and space to the location of each OMI pixel center. U-and V-(west-east and south-north, respectively) wind-speed components were first averaged in the vertical 25 between 0 and 1 km where the majority of the SO2 mass resides. The wind components were then interpolated spatially and temporally to the location and overpass time of each OMI pixel.
Note that to reconstruct annual mean VCD maps based on annual emissions (section 3.4), it is not necessary to have the actual year-specific meteorological information, as annual mean wind characteristics do not vary much from year to year (see the Supplement, section S6), and so for convenience, we simply used wind data from 2005 for all years prior to 2005. taken into account when emissions are allocated. Therefore, the MACC-III point source emission data should be regarded as estimates that may differ considerably from the actual emissions.

SO2 surface concentration data
In-situ SO2 ground-level measurements from the U.S. Clean Air Status and Trend Network (CASTNet) (Baumgardner et al., 1999;Park et al., 2004;Schwede et al., 2011), operated by the U.S. EPA (http://www.epa.gov/castnet), and the Canadian Air 30 and Precipitation Monitoring Network (CAPMoN: http://www.ec.gc.ca/rs-mn/default.asp?lang=En&n=752CE271-1) (Schwede et al., 2011), operated by Environment and Climate Change Canada (ECCC), were used in this study. Both 6 networks were established to assess regional trends in pollutant concentrations, atmospheric deposition, and ecological effects due to changes in air pollutant emissions. CASTNet started operations in 1987 and CAPMoN started in the late 1970s. Both networks employ filter packs to measure SO2, although CASTNet uses a one-week sampling period vs. a oneday sampling period for CAPMoN. It is important to note that the monitoring sites belonging to these networks are located in relatively remote areas, so that direct impacts of local pollution sources on the measurements are minimal. Annual mean 5 SO2 values in μg m -3 were used in this study.

Linking satellite SO2 VCDs and SO2 emissions
The method for linking OMI SO2 VCDs to SO2 emissions is based on a fit of OMI VCDs to an empirical plume model developed to describe the SO2 spatial distribution (as seen by OMI) near emission point sources (Fioletov et al., 2015), but unlike the previous studies, it is not limited to a single point source. The plume model assumes that the SO2 concentrations 10 emitted from a point source decline exponentially with time and that they are affected by turbulent diffusion that can be described by a 2D Gaussian function. The overall behaviour can be described as a combination of exponential and Gaussian random variables, also known as an exponentially modified Gaussian function (see the Appendix for details). Each satellite measurement (or pixel) is fit by a sum of plumes from all point sources. The distribution of SO2 emanating from each source is described by the plume model based on a known plume function Ω (θ, φ, ω, s, θi, φi) dependent on the satellite pixel 15 coordinates (θ, φ), pixel wind direction and speed (ω, s), and source coordinates (θi, φi) scaled by an unknown parameter (αi) representing the total SO2 mass from the source i. These unknown parameters are then estimated from the best fit of the OMI measurements. The emission rate for source i is E=αi/τ, where τ is a prescribed SO2 decay time. In other words, the method finds the emission rates that produce the best agreement with the observed OMI SO2 VCD values. The detailed formulas and prescribed seasonal decay times are given in the Appendix. 20 Thus, the fitting procedure allows for the isolation of the emissions-related "signal" in the data from known sources and can be used to check existing point-source emissions inventories. If all sources are included in the fit, it can be expected that the difference between the OMI data and the fit is within the noise level and the estimated emission rates E should agree with the reported emissions. We used OMI observations and emissions data for the eastern U.S. and southeastern Canada to confirm this expectation. Sources that are not included in the fit would appear as "hotspots" on the maps of the difference 25 between OMI VCDs that could be used for source detection. Furthermore, emissions from such sources could then be derived by adding their coordinates to the source list in the fitting procedure. The suggested method can thus be used as a source of independent emission estimates in regions where emissions values have large uncertainties.
The method requires information about the point-source locations. We used source location data available from the US and Canadian emission inventories mentioned in section 2.3. As discussed by (Fioletov et al., 2015), sources that emit 30 30 kt yr -1 or more can be detected by OMI. Since multiple smaller sources located in a close proximity can also be seen as a "hotspot" in OMI data, we lowered the minimum limit and included all SO2 point sources that reported emissions of 20 kt yr -7 detectability, it gives more accurate emission estimates for clusters of small sources where the point source algorithm is not really applicable.
Earlier versions of the OMI SO2 data product have some large-scale biases (Fioletov et al., 2011) that were largely removed in the present PCA version. However, we found that even the PCA version has some local biases that may interfere with the regression fit. The local bias can be accounted for by introducing functions that change slowly (compared to signal 5 from emission sources) with latitude and longitude. We used Legendre polynomials of latitude and longitude and their products that are orthogonal over the analysed domain as discussed in the Appendix.
The OMI data with and without the bias and the fitting results for four multi-year intervals are shown in the columns I and II of Figure 1. The additional plots of the bias itself and the residuals are available from the Supplement, Figures S1 -S4. Figure 1 is based on the annual estimates averaged over two-and three-year periods. Figure 1 confirms that 10 there was a large decline in SO2 VCD over the eastern U.S. and southeastern Canada in the period 2005-2015 . In contrast, the bias estimated from the fitting procedure appears to be fairly constant over time ( Figure S1), which suggests that it may be an artifact from the retrieval. The lack of this feature in OMPS observations further suggests it is a bias in OMI PCA data as discussed in the Supplement (Section S1 and Figure S2).
It should be mentioned that the use of an empirical plume model is appropriate when atmospheric 15 advection/diffusion can be considered to be the dominant process and meteorological conditions can be assumed to be quasisteady. This is a reasonable assumption for short time periods and transport distances and when chemical transformation and surface removal of SO2 can be well represented as simple first-order loss. The consistent mid-day overpass time for OMI means that the vast majority of the satellite measurements will be associated with a well-developed quasi-steady planetary boundary layer. A 3D atmospheric chemistry model, on the other hand, would be more appropriate for longer time periods 20 and transport distances and for emissions occurring at all times of day, but that is not the case for this analysis.

SO2 emission estimates from OMI data
The functions Ω (θ, φ, ω, s, θi, φi) decline very rapidly with distance from the source located at θi, φi. For an isolated point source (θi, φi) where other sources are located 100 km away or more, Ω (θ, φ, ω, s, θi, φi)is not correlated with any other Ω 25 (θ, φ, ω, s, θj, φj,), where i ≠ j, and the regression model (1) or (2) can be simply split into two parts: a model for point-source emission estimates for source i and a model for all other point sources. Then the estimate of αi is independent from estimates for all other sources. If, however, there is another source j located at (θj, φj) that is closer to source i than ~100 km, then the functions Ω (θ, φ, ω, s, θi, φi) and Ω (θ, φ, ω, s, θj, φj,) become correlated, as do their estimates of αi and αj. As the two Ω functions depend on the wind, the correlation coefficients also depend on the wind distribution and the locations of the 30 sources relative to the prevailing wind direction and to each other, but the separation distance is the dominant factor. Typical absolute values for the correlation coefficients are about 0.2, 0.6, and 0.8 for distances between sources of 100 km, 50 km, 8 and 25 km, respectively. A high correlation means that, in practice, emissions estimates for sources located in close proximity have large uncertainties as we may have difficulty separating signals from the individual sources. However, if sources i and j are located in close proximity to each other but far from all other sources, then their combined emissions can still be estimated accurately. Thus, such sources can be grouped into clusters, where the member sources are located in close proximity (20-40 km) but the clusters themselves are well-separated and total emissions from each cluster can be estimated 5 from satellite data.
Another way of grouping sources into clusters is to establish a grid over the analysis region and then sum up This grid-based approach can be potentially used for area sources or when the locations of sources are not well known. For illustration, we used VCD measurements over the same area, but assumed that that it is an area source with no 25 individual point sources. If we set a regular grid and assume that each grid point is a "source", we can estimate emissions from such "sources" as described above. VCD can then be calculated using these estimated emissions. Such reconstruction for a 0.5° by 0.5° grid is also shown in Figure 1 (column V) and demonstrates a good agreement with the measured VCD values. Note that the grid spacing should not be too large or else the areal emissions will be underestimated. Likewise, if it is too fine adjacent grid cells will be highly correlated and may result in artificial structure. As Figure 1 (column IV) shows, 30 the fitting results based on emissions are very close to the OMI fitted data (column II). We used the OMI data with local bias removed because with this approach, any instrumental local bias will be interpreted as an area source, resulting in overestimation of emissions.

9
Emissions estimated by this gridded method are also shown in Figure 2. Their uncertainties are higher than for the case of known source locations but are still reasonable. The standard deviation of the difference between the emission estimates for all 1° by 1° cells within the domain area shown in Figure 1 and reported SO2 emissions for the same cells are 54, 37, and 56 kt yr -1 for spring summer, and autumn, respectively. High measurement errors and data gaps prevent estimation of the emissions for winter. 5 The uncertainties of satellite-based emission estimates have been discussed in our previous studies (Fioletov et al., 2015(Fioletov et al., , 2016. They can be as high as 50%, but the two largest contributors to this uncertainty, the airmass factor (AMF; determined by the assumed vertical profile, surface reflectivity, and viewing geometry) and the prescribed lifetime, are related to site-specific conditions and can be considered primarily as systematic. They introduce a scaling factor in estimated emissions that affects absolute values but not relative year-to-year changes in emissions. Moreover, the constant, effective 10 AMF embedded in the OMI SO2 product is based on measurements taken in the eastern US, and the lifetime estimates used here are based on data from the US power plants as well, so these errors are minimal for this region. To further support this claim, AMF values were recalculated for all SO2 observations used in Fioletov et al. (2016) and its impact on these sources was found to minimal, typically less than 5%.

SO2 VCDs estimated from reported emissions 15
The equation that links emissions and VCDs (A1) can also be used for forward calculations: if coefficients αi are known, then SO2 VCDs can be calculated for any location for given wind conditions and these daily VCDs can be averaged to give annual or seasonal means for the analyzed area. Since αi = Ei·τ, and τ is prescribed in our calculations, the available emission inventories that contain Ei can be used to calculate αi. In this case, there is no need to do any fitting or to use any OMI measurements to calculate VCDs. In practice, we can simply use the reported emission data and available OMI pixel 20 locations merged with the wind information and calculate VCDs for each OMI pixel based on its center coordinates. OMI provides daily near-global coverage, and of course, no pixel screening is required for such forward calculations, so it would be essentially a reconstruction of daily VCD maps with spatial resolution of about 15 km by 35 km (approximately the average size of the OMI pixel used in this study) assuming a constant emission rate. Figure 1 (column IV) also shows the result of such annual reconstructions averaged over 2-to 3-year periods. 25 Annual point-source emissions from the EPA and NPRI inventories were used as inputs. The agreement of the reconstructed VCDs with OMI data (with the bias removed) is very good, and the agreement with the OMI data fitting results is truly remarkable. To characterize the overall agreement with the OMI data, fitting results, and reconstructed emissions-based VCDs, a 1° by 1° grid was established and various statistical characteristics were calculated for the gridded data. The standard deviation of the residuals ε for this grid is 0.025 DU, i.e., about 20 times less than the uncertainty of individual OMI 30 measurements. The standard deviation of the difference between the OMI-fitted and the reconstructed emissions-based VCDs is 0.016 DU. Figure 3 shows the scatter plots between the annual VCDs reconstructed from emissions and the three OMI-based data sets shown in Figure 1 for all years. The correlation between the VCDs reconstructed from emissions with the actual OMI data is 0.75, but it rises to 0.91 after the local bias is removed and to 0.97 after the emission-related signal is extracted from the OMI data by the fitting procedure (the first term of equation A4). Moreover, values of the latter correlation coefficient are above 0.88 for all seasonal averages (excluding winter) and they are substantially higher than the correlation 5 coefficients with the actual seasonal OMI data. This result could be used to extract an emissions-related SO2 signal from the OMI data when the signal is weak compared to the noise level but the source locations are known. Additional information is available from the Supplement, Section S2, including a figure of the difference between the fitted VCDs and the reconstructed VCDs as well as seasonal and annual statistics. Figure 1 shows the fitting results in geographical coordinates, i.e., the first term of equation (A4) from the 10 Appendix was calculated for each OMI pixel without any stratification by the wind speed and direction. However, the fitting itself is done in a four-dimensional space where the wind speed and direction are the other two coordinates. To illustrate the fitting results for different wind speeds, Figure 4 shows the original mean OMI SO2 values (with the bias removed) and the fitting results when the data are binned by the wind speed. Note that the fitting parameter estimate was done using data for all wind speeds and the binning applies only to the fitting outputs. In other words, the first term of equation (A4) from the 15 Appendix was calculated using only OMI pixels where the wind speed was within the selected range. The calculations were done for three wind-speed bins for the 2005-2007 period when the SO2 emissions were the highest and the measurements were not affected by the "row anomaly". The wind-speed modal value is about 10 km h -1 , and the first bin represent calm conditions, the second bin contains measurements taken within ±5 km h -1 from the modal value, and the last bin corresponds to relatively high wind speeds. As Figure 4 demonstrates the fitting results are able to capture the changes in SO2 distribution 20 at different wind-speed bins. When the wind speed is low, SO2 values are high over the sources, while the plume spreads out over a larger area when the wind speed is high. The figure also shows that SO2 VCD values measured over the sources, or integrated over a small area around the source, are not good proxies for the emissions because they depend on the wind speed.

Applications for other regions 25
Direct SO2 emissions measurements are not available for many regions of the globe. The described method can be used to verify or even estimate SO2 emissions for other regions. To test this method further, we applied it to the European region using European Pollutant Release and Transfer Register (E-PRTR) and TNO-MACC emission inventory data (see Section 2.3). Figure 5 is similar to Figure 1, but for a part of Europe where the majority of the SO2 sources are located. Sources that emitted more than 10 kt in any year between 2005 and 2014 are shown on the map as black dots. The limit was lowered to 30 Figure 5 also shows a good general agreement between the OMI data and VCDs estimated from emissions. Both show a substantial SO2 VCD decline over most regions, with SO2 values the highest at the beginning of the analysed period (Spain, Romania, Bulgaria, Greece), No major changes are observed by OMI for power plants in Serbia and in Bosnia and Herzegovina, and they are now producing the highest SO2 VCD values over the domain shown. As their emissions are not in the E-PRTR database, TNO-MACC emission inventory data were used instead. 5 The method produces estimates for individual sources that can be further grouped in different ways. Estimated and reported annual emissions for the period 2005-2014 were grouped by nation for nine countries with large SO2 emissions, as shown in Figure 6. There is good agreement qualitatively between the reported and estimated emissions. Some differences in absolute values are expected due to possible multiplicative biases in OMI-based estimates (from the airmass factor and potential errors in τ). In some cases, however, a possible deficiency in the reported emissions cannot be ruled out.

Reconstruction of the past VCD distribution
If detailed emission data are available, it is also possible to calculate emissions-based VCD maps using Equation (A3) for years before the launch of OMI. Figure 7 shows the annual mean VCD maps over the eastern U.S. and southeastern Canada 25 reconstructed from the emissions inventories available since 1980. All point sources (shown by black dots) that emitted more than 1 kt of SO2 in at least one year during the 1980-2015 period were included in the calculations for a total about 380 sources. Note that we slightly expanded the domain area in all directions to include sources that emitted large SO2 amounts prior to the OMI launch. There are two major periods of dramatic changes in SO2 VCD values: first, in the early 1990s,

SO2 surface concentrations and VCDs
Multi-year mean surface SO2 concentrations at stations belonging to the CASTNet and CAPMoN networks (see Section 2.4) were compared to the estimated VCD values. Maps of multi-year mean surface SO2 concentrations at stations belonging to 5 the CASTNet and CAPMoN networks are shown in Figure 8. The color scheme of Figure 8 was chosen to be comparable to that used in Figure 1  The surface-concentration-to-VCD ratio ultimately depends on the shape of the SO2 vertical profile. The shape could be affected by boundary-layer height, site elevation, and perhaps some local conditions. There are, however, some common features in the ratio distribution. As shown in Figure 9d, the ratio is low in areas of high emissions-based VCDs 30 and low in areas where emissions-based VCDs are low. Of course, it is not the mean VCD value itself that affects the ratio, but proximity to emission sources. Figure 9d is based on VCDs derived from emissions, but the same analysis for OMImeasured VCD demonstrates similar results (the Supplement, Section S5). It may be possible to reconstruct surface 13 concentration distribution from VCDs and additional information such as the planetary boundary layer height (Knepp et al., 2015), but such estimates are outside of the scope of this study.

Summary and discussion
Fitting OMI SO2 VCD data by a linear combination of functions, where each function represents the plume from an individual source, makes it possible to estimate emission from these sources or groups of sources. If the location of all 5 sources is known, it is expected that the fitting results and the actual OMI data will agree within the noise level as was found to be the case for the eastern U.S. and southeastern Canada. The same agreement is also observed for this region if the reported emissions are used to calculate VCDs. This suggests a simple way of interpreting satellite SO2 VCD data: they should agree with VCD estimates based on available emission inventories or the fitting results based on known source locations. 10 By applying a statistical plume model (developed from satellite SO2 measurements) to U.S. and Canadian annual SO2 point-source emissions inventories, we were able to reconstruct past annual mean VCDs for the period 1980-2015.
High correlation coefficients between the reconstructed VCDs and the OMI-based values (0.91 for OMI data with local bias removed) for the period 2005-2015 gives us confidence in both data sets. It also demonstrates that the reported changes in It can also provide emissions information for regions where there are no other information sources available. Unreported point and area sources can be detected and emissions from them can be estimates by subtracting VCDs calculated from available emission inventories from satellite VCD measurements, although emission inventories with good spatial resolution would be required for such an analysis. While this study is focused on SO2, the methods can be applied to other species with 25 relatively short lifetimes measured from space, particularly to NO2 and NH3.
We have also applied the method to Europe. The results strikingly illustrate the positive impact of EU legislation; the countries where no decreasing trends are observed are non-EU member states surrounded by EU countries with decreasing emissions. In general, the satellite-based results confirm the trends in reported SO2 emissions from EU member states over the period 2004-2014, but some discrepancies were found that deserve further attention. In one case, for 30 example, it seems that reported emissions already take into account certain planned or foreseen measures, but real-world (satellite-observation) estimates suggest that implementation of these measures was delayed by several years. Moreover, although the trend is clearly followed, the absolute emission levels suggested by the OMI SO2 VCD fitting method are 14 sometimes substantially above the reported emission levels for recent years (Figure 6). Whether these differences are due to underreporting or to methodological issues requires further study.
There are certain limitations to the suggested methods. Satellite SO2 VCD data may still contain local biases that will interfere with emissions estimates or will themselves be interpreted as a source. As the OMI and OMPS data show, these biases could be different from instrument to instrument. Moreover, data from the same OMI instrument could have different 5 biases if processed by different algorithms (Fioletov et al., 2016;their Figure 3). Although the biases could be partially removed using, for example, a constant (for a small fitting area) or polynomial (for larger areas) fit, further improvement of retrieval algorithms is required to eliminate the bias problem. The biases could be particularly large over regions of high SO2 VCD values such as the Persian Gulf and China, so the method should be applied there with caution. The method is also based on the assumption that all SO2 is located near the surface, which determines the wind data used for the fitting. This 10 may not always be the case for very large sources where SO2 can be lifted into the free troposphere. Finally, the plume model itself may not be optimal in some cases.

Data availability
OMI PCA SO2 data used in this study have been publicly released as part of the Aura OMI Sulphur Dioxide Data Product (OMSO2) and can be obtained free of charge from the Goddard Earth Sciences (GES) Data and Information Services Center. 15 This appendix contains a description of the fitting algorithm used to estimate emissions from point and multiple sources. The algorithm for point sources was previously published by Fioletov et al., (2015), but we briefly repeated it for reader's convenience.
Fitting algorithm, point source 5 The first step of the fitting algorithm involves a rotation of the location of each OMI pixel about the source such that, after rotation, all have a common wind direction. Then, the method assumes that concentrations of SO2 emitted from a point source decline exponentially (i.e., as exp(-λt)) with time (t) with a constant "lifetime" (or decay rate) τ=1/λ. In the absence of diffusion and with a constant wind direction and speed (s), SO2 is transported downwind (along the -y axis in the chosen coordinate system) with a concentration that declines exponentially with the distance from the source. Since t= -y/s, this 10 decay is simply exp(λy/s) or exp(λ1y) where λ1= λ/s. Likewise, if the wind speed is zero, the distribution of SO2 near the source is governed by diffusion or, more generally, random fluctuations, and can be described by a two-dimensional Gaussian function of the distance from the source that depends on one parameter σ. As both exponential decay of the concentration along the y coordinate and diffusion take place, the overall behaviour can be described as a combination of exponential and Gaussian random variables, also known as an exponentially modified Gaussian function. Therefore, the 15 statistical model of the SO2 plume employed near the point source has the form of a Gaussian function f(x, y) multiplied by an exponentially modified Gaussian function g(y, s): , where x and y (in km) are the coordinates of the OMI pixel center across and along the wind direction, respectively, and s (in km h -1 ) is the wind speed at the pixel center.
The model depends on two parameters, the decay time (τ), and the plume width (σ). It should be multiplied by a scaling factor α that is proportional to the emission strength. and .
The Gaussian function f (x, y) represents the distribution across the wind direction line. The function g (y, s) is essentially a convolution of Gaussian (determined by the plume width σ) and exponential functions (determined by λ1 related to the lifetime) and represents an exponential decay along the y axis smoothed by a Gaussian function: when σ is close to 0, then g(y, s) ≈ λ1 exp (λ1y) where y≤0. The wind speed s is included in (A1) only through λ1= λ/s. Note that σ1 was used in f (x, y) instead of σ. The value of σ1 increased with the distance from the source to reflect an additional spread of the plume downwind (i.e., when y<0). 5 Parameters σ, λ, and α, can be derived from the fit of the OMI observations by the function α Ω(x,y,s) , i.e., from a nonlinear regression model. However, if the values for σ and τ=1/λ are prescribed, then the remaining value, α, can be determined from a simple linear regression model. The function Ω depends on pixels coordinates in the Cartesian coordinate system related to the wind direction with the center at the analysed source. These coordinates can be derived from pixel latitude (θ) and longitude (φ), the wind 15 direction (ω), and the source latitude (θ0) and longitude (φ0), i.e., Ω(x,y,s) = Ω (θ, φ, ω, s, θ0, φ0). As OMI measurements were merged with the wind data, OMI SO2 VCD at each pixel can therefore be interpreted as a four-dimensional function OMI SO2 (θ, φ, ω, s). The dependence of Ω on the model parameters τ and σ is rather complex and we can simplify this approach by assuming that τ and σ are identical for all sources in the analysed region and only the parameter α differs from source to source (see sensitivity analysis in reference (Fioletov et al., 2016)). Values of τ and σ were selected based on 20 previous estimates for point sources in the eastern U.S. (Fioletov et al., 2015) with some seasonal adjustments: τ values were =5.6, 6.3, 7.7, and 6.3 hours for winter, spring, summer, and autumn respectively. The plume width σ=18 km is dependent on multiple factors, but mostly on the OMI pixel size.

Fitting algorithm, multiple sources
In case of multiple sources with prescribed τ and σ, OMI SO2 VCD can be expressed as a sum of contributions αi·Ωi from all 25 individual sources (i). If (xi, yi) and (x′i, y′i ) are the pixel's Cartesian coordinates (km) in the system with the origin at the source i before and after the wind rotation respectively, then they can be calculated from the pixel and source latitudes and longitudes as: xi= r·(φ-φi)·cos(θi); yi= r·(θ-θi); 30 x′i = xi · cos(-ω) + yi · sin(-ω); 2 SO OMI y′i = -xi · sin(-ω) + yi · cos(-ω); where r=111.3 km·π /180; φ and θ are the pixel longitude and latitude; ω is the pixel wind direction (0 for north); φi and θi are the source i longitude and latitude (all in radian).
Then, similarly to equation (S1), the contribution ai·Ωi from the source i can be expressed as αi·Ωi = αi·f(x′i, y′i) Thus, OMI SO2 VCD can be expressed as a sum of contributions from all individual sources (i) plus noise (ε): where only parameters ai are unknown. Equation (S3) represents a linear regression model where the unknown parameters αi can be estimated from the measured variable (OMI SO2) at many pixels and known regressors Ω (θ, φ, ω, s, θi, φi).
Calculations can be done on an annual or seasonal basis (i.e., using all data for a particular year or for a particular season of a year respectively). Emission estimates for shorter time intervals, e.g., monthly emissions, may be possible for large sources, but they appear to be too noisy for the eastern U.S. and southeastern Canada for practical applications. 15 Earlier versions of the OMI SO2 data product have some large-scale biases (Fioletov et al., 2011) that were largely removed in the present PCA version. However, we found that even the PCA version has some biases that may interfere with the regression fit if equation (S3) is used. If the fit is done for a relatively small area, the bias can be accounted for by adding a parameter α0 to the equation (S3) and estimating it from the fit: For a larger area, for example for the eastern U.S. and southeastern Canada, geographic variations in the bias can be accounted for by introducing functions that change slowly with latitude and longitude. We used Legendre Polynomials (Pn(x)) that are orthogonal on the interval from -1 to +1.

25
To make the polynomials orthogonal on the analyzed domain, the following transformation was applied: Lj(θ) = Pj(2·(θ -θmin)/(θmax-θmin)-1); where φmin, φmax, θmin, and θmax are latitudes and longitudes that define the domain area. Then Lj(θ) and Lk(φ), and their products were added to the fit: 5 where αi and βj,k are the estimated coefficients. The first sum represents the emission-related fitting and the second sum is the large-scale bias. Equation (S4) also represents a linear regression model and the unknown coefficients can be estimated from the available observations. Polynomials up to the 6th degree were used for each one-year or one-season fit for the selected domain (the eastern U.S. and southeastern Canada), although a higher (or lower) degree may be more suitable for a larger 10 (smaller) area (see also section S6). Note that the biases are related to retrieval effects such as imperfection of account for the ozone absorption and therefore are not related to SO2 abundances and not affected by the winds. For this reason, no dependence of the bias on s is considered. Figure A1 illustrates the method by using SO2 data from 2005-2007 near the Bowen power plant in Georgia, U.S.
There are 13 sources within the ±200 km square area around the Bowen facility. The fitting was done for every year, 15 estimated values ai·Ω (θ, φ, ω, s, θi, φi) were calculated for each satellite pixel, then summed up to obtain a SO2 VCD value for the fit for that pixel. For Figure S1, the actual OMI data and the fitting results were averaged over the 2005-2007 period and smoothed by the pixel averaging technique with a 30 km radius. The maps of estimated values for individual sources smoothed in the same way are also shown. The map of the residuals, or the difference between the OMI data-based map and the fitting results is also shown. 20  C. R., Chance, K., Liu, X., Lee, C. and Martin, R. V.: Application of OMI, SCIAMACHY, and GOME-2 satellite SO2 retrievals for detection of large emission sources, J. Geophys. Res. Atmos., 118(19), 11,399-11,418, doi:10.1002/jgrd.50826, 2013 included in the fit (they are shown as the black dots). Results of the fitting of OMI data by the set of functions that represent "sources" as 0.5° by 0.5° grid cells (shown as the black dots) using estimated emissions (see text) are shown in column V.
The maps are smoothed by the pixel averaging technique with a 30 km radius (Fioletov et al., 2011). Averages for four multi-year periods, 2005-2006, 2007-2009, 2010-2012, and 2013-2015 over the area 32.5°N to 43°N and 75°W to 89°W are 10 shown.   The area 30°N to 48°N and 70°W to 90°W is shown.