SO2 Emission Estimates Using OMI SO2 Retrievals for 2005-2017.

SO2 column densities from Ozone Monitoring Instrument provide important information on emission trends and missing sources, but there are discrepancies between different retrieval products. We employ three Ozone Monitoring Instrument SO2 retrieval products (National Aeronautics and Space Administration (NASA) standard (SP), NASA prototype, and BIRA) to study the magnitude and trend of SO2 emissions. SO2 column densities from these retrievals are most consistent when viewing angles and solar zenith angles are small, suggesting more robust emission estimates in summer and at low latitudes. We then apply a hybrid 4D-Var/mass balance emission inversion to derive monthly SO2 emissions from the NASA SP and BIRA products. Compared to HTAPv2 emissions in 2010, both posterior emission estimates are lower in United States, India, and Southeast China, but show different changes of emissions in North China Plain. The discrepancies between monthly NASA and BIRA posterior emissions in 2010 are less than or equal to 17% in China and 34% in India. SO2 emissions increase from 2005 to 2016 by 35% (NASA)-48% (BIRA) in India, but decrease in China by 23% (NASA)-33% (BIRA) since 2008. Compared to in situ measurements, the posterior GEOS-Chem surface SO2 concentrations have reduced NMB in China, the United States, and India but not in South Korea in 2010. BIRA posteriors have better consistency with the annual growth rate of surface SO2 measurement in China and spatial variability of SO2 concentration in China, South Korea, and India, whereas NASA SP posteriors have better seasonality. These evaluations demonstrate the capability to recover SO2 emissions using Ozone Monitoring Instrument observations. ©2019. The Authors.

Introductions

Sulfur dioxide (SO2) is a reactive gas‐phase, air pollutant. Exposure to SO2 affects human health by increasing respiratory (Chen et al., 2007; Chen et al., 2012; Clark et al., 2010; Ranzi et al., 2014; Rodriguez‐Villamizar et al., 2015) and cardiovascular morbidity (Chen et al., 2012). In the atmosphere, SO2 is predominantly oxidized by reactions with the OH radical in the gas phase, or by hydrogen peroxide (H2O2) and ozone (O3) in the aqueous phase to form sulfuric acid, the latter which can lead to acid rain formation and particulate sulfate (Seinfeld & Pandis, 2006). Sulfate is a key component of particular matter with aerodynamic diameter less than 2.5 μm (PM2.5) and comprises 10%–67% of aerosol mass in the midlatitude Northern Hemisphere (Zhang et al., 2007). Exposure to ambient PM2.5 is estimated to cause 4.2 million premature deaths in 2015 (Cohen et al., 2017). SO2 also leads to the formation of haze that decreases visibility (Haagen‐Smit, 1952).

The primary anthropogenic sources of SO2 are combustion of fossil fuels and smelters (Kato & Akimoto, 1992; Arndt et al., 1997; Garg et al., 2001; Serbula et al., 2015). Major natural SO2 sources are oxidation of biogenic dimethyl sulfide (Davis et al., 1999) and volcanic eruptions (Carn et al., 2015; Ge et al., 2016; Graf et al., 1997). Bottom‐up anthropogenic SO2 emission inventories from MACC/CityZEN (MACCity, http://ether.ipsl.jussieu.fr/eccad), Emissions Database for Global Atmospheric Research (EDGAR; http://edgar.jrc.ec.europa.eu/), the Greenhouse Gas and Air Pollution Interactions and Synergies (GAINS, https://gains.iiasa.ac.at/gains/),REanalysis of the TROpospheric chemical composition over the past 40 years (RETRO, http://accent.aero.jussieu.fr/RETRO_metadata.php), and Hemispheric Transport of Air Pollution (HTAP, Janssens‐Maenhout et al., 2015) have tracked the magnitude and trend of global SO2 emissions from human activities over the past decades. Bottom‐up inventories show a decrease of SO2 due to emission regulation in the United States (de Gouw et al., 2014), China (Liu et al., 2015), and Europe (Crippa et al., 2016), which also contributes to the decline of sulfate and organic aerosol concentrations (Geng et al., 2017; Marais et al., 2017; Paulot et al., 2016). However, SO2 emissions in India are still on the rise, and will likely continue to increase during the 2020–2050 period if emissions remain unregulated, driven by the growth of electricity demand and industrial production (Venkataraman et al., 2018). Although the bottom‐up inventories provide important information on the magnitude and trend of SO2 emissions, they may have systematic bias and their availability lags in time. For instance, the National Emission Inventories for the United States is released every three years; the most comprehensive anthropogenic emission inventory for China, Multi‐resolution Emission Inventory (MEIC), has only been released for 2008, 2010, and 2012. In addition, uncertainties in bottom‐up SO2 emissions can be significant—in Asia up to 35% (Kurokawa et al., 2013; Lu et al., 2011; Zhang et al., 2009; Zhao et al., 2011). Total SO2 emission estimates in Asia differ by 29% between MIX and EDGAR inventories; the discrepancies are even larger on a regional scale and when considering individual sectors (Li, Zhang, et al., 2017; Saikawa et al., 2017).

Over the past two decades, retrievals of SO2 columns from satellite remote sensing instruments have provided important information about SO2 sources (Burrows et al., 1999; Eisinger & Burrows, 1998; Fioletov et al., 2015; Krueger, 1983) and more timely updates on SO2 magnitudes and trends compared to bottom‐up inventories (Krotkov et al., 2016; Li, McLinden, et al., 2017; Wang et al., 2016; Zhang et al., 2017). These instruments include the Total Ozone Mapping Spectrometer, Global Ozone Monitoring Experiment (GOME), GOME‐2, SCanning Imaging Absorption SpectroMeter for Atmospheric CHartographY (SCIAMACHY), Ozone Monitoring Instrument (OMI), and Ozone Mapping and Profiler Suite (OMPS).

Among these instruments, OMI distinguishes itself by having the finest spatial resolution (footprint of 13 km × 24 km at nadir) and daily global coverage (Levelt et al., 2006). Its SO2 retrievals have been demonstrated to have better performance in detecting smaller sources than SCIAMACHY and GOME‐2 (Fioletov et al., 2013). There are two SO2 column retrieval products available for OMI. The National Aeronautics and Space Administration (NASA) product employs a Principal Component Analysis (PCA) approach, whereas the Belgian Institute for Space Aeronomy (BIRA‐IASB) product uses the Differential Optical Absorption Spectroscopy (DOAS) algorithm (Platt & Stutz, 2008). The two products have been compared and evaluated in Theys et al. (2015). Although both algorithms show comparable SO2 column over industrialized and populated areas, the BIRA product has generally larger values than the NASA product, with a mean difference of 0.2 Dobson unit (DU; 1 DU = 2.69 × 1016 molec/cm2; Theys et al., 2015).

SO2 satellite observations from different instruments have also been shown to have different magnitudes, especially over China. For instance, Fioletov et al. (2013) found similar mean SO2 column concentrations from OMI (using a previous version of the NASA product (Krotkov et al., 2006)), SCIAMACHY, and GOME‐2 retrievals at locations close to SO2 sources, but the magnitude of SO2 columns from the latter two instruments are ~1.5 times larger than the OMI NASA retrieval over parts of China. Yang et al. (2013) found that OMI SO2 columns, calculated using a linear fit algorithm (Yang et al., 2007), are mostly smaller than OMPS data in China, with the largest OMPS SO2 column being 1.65 DU in China, whereas the OMI linear fit value was only 0.8 DU at the same location. Another comparison of OMPS (Yang et al., 2013) with BIRA products over East China shows that OMPS SO2 is a factor of 2 to 4 lower in January of 2013 (Theys et al., 2015).

Differences in SO2 retrievals can propagate into posterior emissions derived using these observations. For instance, using the OMI NASA SO2 product (PCA algorithm), Wang et al. (2016) estimated posterior anthropogenic SO2 emissions in China to be 843 Gg/S in April 2008, whereas Koukouli et al. (2018) obtained a posterior SO2 emission of ~1,542 Gg/S for the same month using the OMI BIRA product. The magnitude of posterior SO2 emissions in Wang et al. (2016) is more consistent with another posterior estimates of 815 Gg/S constrained by MODIS AOD observations for April 2008 (Xu et al., 2013). However, the magnitude of posterior SO2 emissions constrained by the BIRA product is closer to prior emissions from INTEX‐B (Wang et al., 2016; Xu et al., 2013). Differences in the emission estimate method (4D‐Var in Wang et al. (2016), mass balance in Koukouli et al. (2018)), the atmospheric models employed, or other factors such as meteorological inputs may also contribute to these ranges of estimates for SO2 emissions.

In light of these discrepancies in SO2 column retrievals, the factor of ~2 difference in posterior SO2 emissions from different OMI products, and the lack of detailed comparisons and evaluations of the magnitude, seasonality, and interannual variation of the newest OMI SO2 retrievals, in this study we derive posterior SO2 emissions using NASA and BIRA products and evaluate their robustness by analyzing the two retrievals and posterior emissions. We then compare posterior simulations with in situ measurements. We focus on measurements from OMI owing to its high spatial resolution, daily global coverage, the availability of data over 14 years, and access to a new prototype OMI SO2 product from NASA with more detailed vertical profile information than the current NASA standard product version 3.

The goal of this work is to estimate long‐term SO2 emissions using constraints from OMI SO2 products. Given the differences between several OMI SO2 retrievals and the propagation of these discrepancies to emission estimates, we start by comparing SO2 columns from NASA and BIRA products and evaluate causes of differences between them in section 3. In section 4, posterior SO2 emissions are derived at three domains: globally at the 2° × 2.5° resolution, nested U.S. domain at 0.5° × 0.667°, and nested East Asia domain at 0.5° × 0.667°. The magnitude, seasonality, and interannual variation of posterior SO2 emissions derived from the two OMI retrievals are compared in section 4.2. In section 5, the posterior SO2 emissions are used in the GEOS‐Chem chemical transport model to simulate surface SO2 concentrations, which are further compared with measurements.

Observations and Model

OMI SO2 Observations

OMI is a nadir‐viewing ultraviolet/visible (UV/VIS) spectrometer measuring backscattered sunlight in two UV regions (264–311, 307–383 nm) and one visible region (349–504 nm). It is aboard the Earth Observing System (EOS) Aura satellite, with a local equator crossing time of approximately 13:45 in the ascending node and has daily global coverage. OMI was launched in August 2004, and its operational data have been available since October 2004. However, since January 2009, parts of the OMI cross‐track positions are affected by row anomalies, which became more severe in July 2011 (http://projects.knmi.nl/omi/research/product/rowanomaly-background.php). This effectively changes the OMI coverage from daily global coverage to two‐day daily coverage.

We use three different level 2 OMI SO2 retrieval products in this study: the NASA standard product (operational product) OMSO2 (version 3; Li et al., 2013), the NASA prototype product, and the BIRA product (Theys et al., 2015). The OMI SO2 retrieval in the NASA standard product and prototype product employ a principal component analysis of radiances with a fitting window of 310.5–340 nm (Li et al., 2013). Based on different assumptions of the vertical profile and center of mass altitude (CMA) of SO2, four products, that is, planetary boundary layer (PBL) SO2 column with CMA = 0.9 km, lower tropospheric SO2 column with CMA = 2.5 km, middle tropospheric SO2 column with CMA = 7.5 km, and upper tropospheric and stratospheric SO2 column with CMA = 17 km, are provided. Given our focus on constraining anthropogenic SO2 emissions, which is emitted within the PBL, we use the PBL product. We also use two months of NASA prototype data (January and July of 2010—the only two months available at the time of our analysis). This new algorithm uses geometry, O3, and cloudiness information to calculate the air mass factor (AMF) for each individual retrieval and updates snow/ice pixels, changes treatment of small cloud fractions (treats cloud fraction less than 0.05 as cloud free), and improves the cloud radiance fraction calculation compared to NASA standard product. BIRA SO2 retrievals use a radiative transfer calculation based on the DOAS algorithm to fit a window of 312–326 nm (Theys et al., 2015). Its retrievals have good agreement (difference within 12%) with ground‐based MAX‐DOAS SO2 measurements in the county of Xianghe, China when aerosol optical depth is low, but are significantly different from other instruments (as much as 10–15 DU) under severe haze. This is explained by the compensation of reduction and enhancement of SO2 due to aerosols (Theys et al., 2015). In addition to different retrieval algorithms, the NASA and BIRA products also use different cloud and climatology data, listed in Table 1. However, NASA and BIRA retrievals have consistent detection limits of 0.5–0.6 DU (Theys et al., 2015). Treatment of AMFs and scattering weight from these three products are discussed in Appendix A.

Table 1

Summary of Major Differences in NASA and BIRA OMI SO2 Retrievals

	NASA standard product	NASA prototype	BIRA
Retrieval method	PCA	PCA	DOAS
Fitting window	310.5–340 nm	310.5–340 nm	312–326 nm
Cloud product	Assumes cloud free	OMCLDRR	OMCLDO2.003
Surface reflectivity	Fixed value of 0.05	TOMS climatology	OMI climatology
Radiative transfer model	VLIDORT	VLIDORT	LIDORT
Cloud radiative fraction (terms defined in section 3)	0	CRF=fcI_cloudI_measured	CRF=fcfcI_cloudfcI_cloud+1−fcI_clear
Uncertainties	0.5–0.9 DU (https://SO2.gsfc.nasa.gov/Documentation/OMSO2Readme_V120_20140926.pdf)		0.5–1.2 DU (Theys et al., 2015)

We use consistent filtering criteria for all three products for comparisons in section 3. We only include data with values between −10 and 5 DU, not affected by the row anomaly, on descending nodes (due to NASA retrieval algorithm performance requirement), with solar zenith angles (SZA) less than 65°, have a cloud fraction less than 0.2, O3 slant columns less than 1,500 DU, and only use the 10th–50th pixel on the cross‐track dimension. These 10 pixels at each edge have larger footprint (~30 km × 160 km) and uncertainties, and are therefore excluded from our assimilation. The use of different cloud products inevitably filters out different pixels from the two retrievals. For the posterior emission estimates in section 4, the same filtering criteria are applied to the NASA standard product; the O3 slant column and descending node criteria do not affect the BIRA retrieval quality and are not included in our filtering. Data at the north and south edges of the swath are found to be noisy and are excluded from inversion in section 4. Data in regions affected by the South Atlantic Anomaly and dominated by nonanthropogenic sources are also excluded.

GEOS‐Chem and Adjoint Model

The GEOS‐Chem adjoint model (Henze et al., 2007, 2009) was developed for inverse modeling of aerosol and gas emissions using the 4D‐Var approach. It is employed in this study to derive posterior SO2 emissions at three domains. The simulation has a global horizontal resolution of 2° × 2.5° and 47 vertical layers. It uses the Modern‐Era Retrospective analysis for Research and Application, version 2 (MERRA‐2) meteorology fields. MERRA‐2 is produced by the NASA Global Modeling and Assimilation Office GEOS‐5 data assimilation system version 5.12.4. The data are available from 1979 to present. Both nested East Asia (70°E–150°E, 0°N–50°N) and nested U.S. (126°W–66°W, 13°N–57°N) simulations have a horizontal resolution of 0.5° in latitude and 0.667° in longitude, and 47 vertical layers. Meteorological inputs for both nested domain simulations are from Global Modeling and Assimilation Office GEOS‐5 reanalysis fields, which have a temporal coverage from 2004 to mid‐2013, and native horizontal resolution of 0.5° × 0.667° and 72 vertical layers. The sulfur cycle simulation (emission, chemistry, advection, convection, diffusion, dry deposition, and wet deposition) is based on the Goddard Chemistry Aerosol Radiation and Transport model (Chin et al., 2000) and was implemented in GEOS‐Chem by Park et al. (2004). Thermodynamic equilibrium of sulfate, nitrate, and ammonium is calculated through RPMARES (Park et al., 2004). A resistance‐in‐series model is applied to compute dry deposition (Wang et al., 1998; Wesely, 1989). Wet deposition is described in Liu et al. (2001). Previous studies have found slightly smaller SO2 lifetimes in GEOS‐Chem simulations than lifetimes estimated from in situ measurement due to faster SO2 oxidation in the model (Lee et al., 2011). Although there is not yet scientific consensus on the precise chemical mechanism of sulfate formation during extreme winter haze episodes, including aqueous phase oxidation of S (IV) by NO2 only leads to less than 3% changes in mean wintertime surface SO2 concentrations in the United States (Sarwar et al., 2013) and less than 20% changes over East Asia in January 2010 in our tests. These differences and the impact of uncertain chemical mechanisms on SO2 concentrations are much smaller than the impact of uncertainties in bottom‐up and satellite observations in terms of monthly and seasonal averages.

Anthropogenic emissions of NOx, SO2, NH3, CO, non‐methane volatile organic compounds (NMVOCs), and primary aerosols are taken from the HTAP version 2 ( Janssens‐Maenhout et al., 2015). This inventory is compiled using different regional inventories, including Environmental Protection Agency (EPA)'s for the United States, U.S. EPA (for data in the south of 45°N and diurnal variation) and Environmental Canada's for Canada, the Model Inter Comparison Study in Asia's for Asian countries, and the European Monitoring and Evaluation Programme and Netherlands Organisation for Applied Scientific Research (TNO)’s for Europe; the EDGAR v4.3 is used for the rest of the world. The HTAP inventory has a horizontal resolution of 0.1° × 0.1°. Three hourly biomass burning emissions from Global Fire Emissions Database, version 4 (Giglio et al., 2013) are used. These are preprocessed for the GEOS‐Chem adjoint using the HEMCO stand‐alone model (Keller et al., 2014). We apply 2010 anthropogenic emissions for all simulated years. Time‐varying natural emissions of SO2 include oxidation of dimethyl sulfide from oceanic plankton to form SO2 follows Park et al. (2004), where oceanic dimethyl sulfide emissions are calculated as the product of local seawater dimethyl sulfide concentration and sea‐to‐air transfer velocity. We implement volcanic SO2 emissions from Ge et al. (2016) (using OMI SO2 linear fit retrievals from Yang et al. (2007)), although this still results in an underestimation of modeled SO2 columns compared to OMI SO2 BIRA product and NASA standard product (Figure S1). We attribute this remaining discrepancy to the different magnitude of retrieved SO2 columns in Yang et al. (2007) and other products used in this study. Volcanic regions show the largest model‐observation discrepancy. This leads the optimization algorithm to focus on emission adjustments in these regions first. Since we aim to constrain anthropogenic SO2 sources, we mask regions (shown in Table S1, by excluding observations from the assimilation) affected by volcanic emissions. We also replace monthly OMI SO2 column concentrations in months (August 2005 and June 2011 in East Asia and August 2008 and June 2009 in United States) when volcanic eruption lead to underestimates of simulated SO2 over the whole of East Asia and U.S. domain by columns linearly interpolated between OMI column concentrations in the same month in adjacent years.

In Situ Measurements

We use surface SO2 measurements from China National Environmental Monitoring Center, U.S. EPA Air Quality System, India Central Pollution Control Board, and South Korea National Institute of Environmental Research for posterior SO2 emission and simulation evaluations. Pulse UV fluorescent analyzers are used for SO2 measurement in United States and South Korea (Han et al., 2017; NIER, 2016). The measurement methods in China and India may be different at different sites during the studied period. The uncertainties of SO2 measurements are estimated to be 9.1% in South Korea (KCL, 2015), and 6.0%–9.3% based on a monitor in Cuyaghoga County of Ohio State in the United States (EPA, 2017). Uncertainties of SO2 measurements in China are not well documented. Locations of surface measurements and their comparisons with model simulations using bottom‐up emissions are shown in Figure S2.

Aircraft measurements of SO2 during KORea‐United States Air Quality (KORUS‐AQ) campaign (May–June 2016) and Deriving Information on Surface Conditions from COlumn and VERtically Resolved Observations Relevant to Air Quality (DISCOVER‐AQ) campaign (June–August 2011) are employed to evaluate GEOS‐Chem SO2 vertical profile. These measurements are made by chemical ionization mass spectrometer and modified trace‐level pulsed fluorescence analyzer (He et al., 2016), respectively. Uncertainties of SO2 measurements are about 15% at 16.5 pptv for the chemical ionization mass spectrometer (Huey et al., 2004) and range from 92% at 30 pptv to 16% above 300 pptv for the pulsed fluorescence analyzer (Luke, 1997). We compared model simulations with measurements from 23 flights from KORUS‐AQ and 24 flights from DISCOVER‐AQ. Flight tracks of these measurements are shown in Figure 1. One issue for airborne SO2 measurements during the summer when SO2 lifetime is short (less than one day in the eastern United States) is that the measurement is significantly influenced by the stochastic probability of hitting (with a probability of ~10%) or missing plumes from power plants (He et al., 2016; Taubman et al., 2006). These tropospheric SO2 measurements provide important evaluations on our posterior simulations using emissions constrained by satellite observations, which are more sensitive to SO2 aloft (Krotkov et al., 2006).

Figure 1

Flight tracks of KORUS‐AQ measurements overlapped on surface SO2 concentration from the GEOS‐Chem simulation (a) at 2° × 2.5° in May 2016, (b) at 0.5° × 0.667° in May 2010, and (c) flight tracks of DISCOVER‐AQ overlapped on surface SO2 concentration from the GEOS‐Chem simulation 0.5° × 0.667° in July 2011.

Evaluation of SO2 Columns in GEOS‐Chem and OMI Retrievals

Given the importance of satellite retrieval accuracy for constraining SO2 emissions, we compare NASA and BIRA SO2 retrieval products to understand how their differences may indicate the robustness of estimated SO2 emissions in different times and regions. All products consistently show the highest SO2 loadings over China and India. These values are compared to simulated SO2 column concentrations, using different scattering weights and air mass factor treatment as described in Appendix A. We further derive posterior SO2 emissions using equations (A8) and (A9) in section 4.1 (Figure 6 and Table 2) and compare simulated surface SO2 concentrations using these posterior SO2 emissions with surface measurements in section 5 (Table 3). The comparison of the NASA SP and the NASA prototype product in this section provides a means of evaluating the model‐observation differences calculated using constant AMF with the NASA SP.

Figure 6

Monthly SO2 emissions in China, South Korea, India, and United States (normalized to annual mean) in 2010 from nested simulations (0.5° × 0.667°). Posterior emissions are from 4D‐Var estimates.

Table 2

Annual Budget of SO2 Emissions in 2010 (Posteriors Are From 4D‐Var; Tg/S)

	Prior	BIRA posterior	NASA SP posterior	NASA SP posterior (equation (A8))
China	12.39	11.83	11.21	9.48
India	4.36	3.24	3.56	2.99
United States	0.42	0.38	0.38	N/A

Table 3

Comparisons of Annual Mean SO2 Surface Measurements in China With GEOS‐Chem Surface Layer SO2 Concentrations in 2010

		In situ	Prior (HTAP)	NASA SP posterior	BIRA posterior	NASA posterior (equation (A8))
China	SO2 (μg/m3)	43.57	48.34	30.07	39.66	27.61
	R		0.19	0.22	0.28	0.17
	NMB		0.11	−0.31	−0.09	−0.37
	NMSE		0.77	0.82	0.59	0.94
China provincial capitals	SO2 (μg/m3)	46.80	76.35	34.93	48.03	30.85
	R		0.16	0.44	0.30	0.25
	NMB		0.63	−0.25	0.03	−0.34
	NMSE		1.12	0.33	0.38	0.55
East China provincial capital	SO2 (μg/m3)	44.82	72.96	35.03	59.51	40.26
	R		0.71	0.73	0.76	0.79
	NMB		0.63	−0.22	0.33	−0.10
	NMSE		0.50	0.22	0.27	0.13
South Korea	SO2 (ppbv)	5.04	3.25	2.47	3.05	2.95
	R		0.29	0.32	0.37	0.30
	NMB		−0.35	−0.51	−0.40	−0.41
	NMSE		0.36	0.74	0.42	0.48
South Korea large cities	SO2 (ppbv)	5.12	3.79	3.17	3.28	3.51
	R		0.79	0.76	0.73	0.78
	NMB		−0.26	−0.38	−0.36	−0.31
	NMSE		0.13	0.26	0.24	0.18
India	SO2 (ppbv)	10.95	23.40	10.99	9.01	11.77
	R		0.34	0.14	0.36	0.35
	NMB		1.14	0.004	−0.18	0.08
	NMSE		2.21	0.89	0.58	0.62
United States	SO2 (ppbv)	1.93	3.59	3.44	3.58
	R		0.18	0.19	0.18
	NMB		0.86	0.78	0.85
	NMSE		2.62	2.54	2.58
U.S. large cities	SO2 (ppbv)	2.04	2.03	1.67	1.74
	R		−0.15	−0.17	−0.17
	NMB		−0.003	−0.18	−0.15
	NMSE		2.44	2.35	2.35

Any sites that have measurements in 2010 are included. Statistics are calculated over all monitoring sites, that is, by taking the annual mean of surface measurements and GEOS‐Chem simulation sampled at each monitoring location, and calculating statistics across the sites. The 10 large cities in South Korea are Seoul, Busan, Incheon, Daegu, Daejeon, Gwangju, Suwon, Goyang‐si, Seongnam‐si, and Ulsan (http://www.geonames.org/KR/largest-cities-in-south-korea.html). The 10 large cities in the United States are New York City, Los Angeles, Chicago, Houston, Philadelphia, Phoenix, San Antonio, San Diego, Dallas, and San Jose (https://public.opendatasoft.com/explore/dataset/1000-largest-us-cities-by-population-with-geographic-coordinates/table/?sort=-rank).

Figure 2 shows SO2 SCDs from three OMI retrievals (second row) and GEOS‐Chem (first row). GEOS‐Chem SCDs calculated with the NASA AMF (equation (A8); first column) and scattering weight (equation (A4); second column) are larger than those using the BIRA scattering weights (equation (A4); third column). The use of a constant AMF of 0.36 in the NASA standard product neglects different viewing geometry and cloud conditions and lead to uncertainties in the posterior SO2 emissions. The mean percent differences between OMI AMFs from the NASA prototype product and 0.36 are 162% for China and 124% for India in January, 2010. Fioletov et al. (2016) estimated new OMI AMFs for large SO2 sources at 491 sites using site‐specific surface reflectivity, solar zenith angle, viewing geometry, surface pressure, cloud fraction, pressure, and SO2 profile shape. The mean percent differences between these OMI AMFs and 0.36 are 33% over all sites, 40% over sites in China, and 20% over sites in India. The mean percent differences between AMFs calculated following Wang et al. (2016) and 0.36 depends on the month, and amount to 341% over China and 173% over India in 2010. However, the accuracy of the AMFs calculated from the look‐up table is subject to further evaluation.

Figure 2

East Asia SO2 SCDs from (first row) GEOS‐Chem, (second row) OMI, and (third row) the difference between the two for January 2010. GEOS‐Chem SCDs are sampled within half an hour of the OMI overpass time and at grid cells that contain OMI footprints. This comparison in the base year reflects differences in SO2 columns only caused by emissions.

In addition to scattering weight and AMFs, GEOS‐Chem SO2 columns are also affected by filtering using different cloud products. For instance, differences between GEOS‐Chem simulated SO2 SCDs (calculated using scattering weights from BIRA retrievals) filtered using NASA and BIRA cloud products vary between −0.70 and 0.71 DU, with a mean of 0.007 DU over East Asia domain in January 2010. BIRA OMI retrievals generally have larger SO2 SCDs than NASA retrievals, which is consistent with the difference between these two of +0.2 DU over southeastern China reported in Theys et al. (2015). The combination of these two effects leads to larger model‐observation differences in the NASA retrieval than the BIRA retrieval in China. The magnitudes of OMI SO2 column concentrations in the prototype retrieval are closer to BIRA than the NASA SP retrieval. The retrievals and model estimates are more consistent in India. Similar analyses for July, as well as the U.S. and global domains for January and July, are provided in section S1.

From these global and regional domain comparisons, the model‐observation difference is more consistent in all three retrievals in low‐latitude regions (smaller view zenith angle) and in July in North Hemisphere (smaller SZA), suggesting that more robust posterior emission estimates may be obtained in these regions and seasons. This difference in China (highly polluted with high aerosol loadings) has the most obvious variability across the three products in both January and July, consistent with previous comparisons of SO2 columns from different instruments (Fioletov et al., 2013; Theys et al., 2015; Yang et al., 2013) and posterior SO2 emission estimates using different OMI products (Koukouli et al., 2018; Wang et al., 2016).

There are several factors that contribute to the differences in NASA and BIRA SO2 retrievals. To understand under which conditions these retrievals are the most consistent and may lead to more robust SO2 emission estimates, we descriptively analyze discrepancies in cloud products and cloud radiative fraction calculations as well as their impacts on the retrievals. A critical factor in the differences shown across the retrieved and simulated SO2 SCDs are the AMFs, which are evaluated in section S2. In brief, NASA AMFs are higher in both January and July than BIRA AMFs owing to different scattering weights, for several reasons.

Both the cloud product and the cloud radiative fraction calculation lead to differences in the vertical profile of scattering weights . Cloud products used in the different retrieval algorithms have different cloud top pressures and heights (Figures S13–S15). Scattering weights below clouds are treated as zero, leading to zero GEOS‐Chem SCDs below clouds. Also, both NASA and BIRA algorithms treat aerosol optical effects implicitly as “effective” clouds, leading to larger uncertainty in retrievals in regions with high aerosol loading, that is, China (Liu et al., 2018). The cloud radiance fraction (CRF) is used to weight scattering weights and AMFs between cloudy and clear parts of a pixel, and is calculated differently in two products. The NASA product uses where f is effective cloud fraction, I is the top‐of‐atmosphere radiance caused by cloud, and I is the measured top‐of‐atmosphere radiance. The BIRA product uses where I is the top‐of‐atmosphere radiance for a cloud free scene. As shown in Figure 3, SO2 sensitivities from NASA are larger than BIRA at most altitudes for all cloud fractions, leading to larger GEOS‐Chem SO2 SCDs in NASA retrievals.

Figure 3

Examples of vertical distribution of SO2 sensitivity (m) in the BIRA and NASA retrievals under (a) clear‐sky conditions (cloud fraction = 0, VZA = 56.61°), (b) cloud fraction between 0 and 0.1 (VZA = 56.62°), and (c) cloud fraction between 0.1 and 0.2 (VZA = 62.08°).

Another factor correlated with differences between the OMI SCDs from the NASA and BIRA products is the view zenith angle. Locations where BIRA OMI SO2 column values are large usually correspond to view zenith angles greater than 30° (Figure S17). This correlation is possibly caused by different radiative transfer calculations. Differences in climatology and cloud fraction input to these calculations (Table 1) can also lead to different scattering weights.

We next compare GEOS‐Chem SO2 SCDs calculated using NASA (horizontal axis) and BIRA (vertical axis) scattering weights of the first retrieval falling in each model grid cell in East Asia domain (Figure 4). The use of first retrieval ensures that only one GEOS‐Chem SCD per grid cell is used in the comparison instead of the average of several SCDs. The smallest NMSE and largest correlations between the two simulated SO2 SCDs were found for cases without clouds. Since the BIRA retrieval sets the cloud fraction to zero when the cloud fraction is between 0 and 0.1, the differences in the second column reflects the impact of clouds on the NASA scattering weights, which are zero below clouds, leading to lower SO2 columns. The third column shows the impact of different cloud top height and cloud radiative fraction from both products. In contrast, the correlation of the actual OMI SO2 SCDs (Figure S18) is smaller and the bias is larger than GEOS‐Chem SCDs. This is caused by the use of a different retrieval method and fitting window (310.5–340 nm for the NASA product and 312–326 nm for the BIRA product), and the use of a background correction scheme in the DOAS algorithm which is not used in the principal component analysis method (Theys et al., 2015).

Figure 4

Comparison of GEOS‐Chem SO2 SCDs (DU) for different cloud fractions (CF) for January 2010.

In summary, we find that different cloud top heights and CRF calculations in the NASA and BIRA retrievals lead to different scattering weights and thus GEOS‐Chem SO2 SCDs, which contribute to model‐observation discrepancies. OMI SCD differences stem from use of different radiative transfer models and several other elements in the algorithm as listed in Table 1. Based on our evaluation, the differences between model and observations are more consistent at lower latitudes and during the summer, times and locations where we thus expect to obtain more robust emission estimates.

Posterior Anthropogenic SO2 Emissions

In this section, we derive global posterior SO2 emissions (2° × 2.5°, 2005–2017) using both the NASA and BIRA products. We also derive higher resolution (0.5° × 0.667°) SO2 emissions over East Asia and the United States from 2005 to 2012 (limited by the availability of GEOS‐5 meteorology fields) and use these for evaluation with surface measurements in section 5. Reasons for deriving these high‐resolution inversions are the large SO2 burden in East Asia and the United States, the interest in the spatial distribution of SO2 emissions inside these regions, and the impact of finer model resolution on comparisons with surface measurements. As only two months (January and July 2010) of the NASA prototype data were available, and the sign and spatial distribution of the model‐observation differences are quite similar when using the NASA prototype and SP (Figure 3 in section 3), we use the SP retrieval to estimate SO2 emissions.

For all three domains, we use the hybrid inversion approach developed in Qu et al. (2017), wherein we perform a 4D‐Var inversion for a base year (2010) and then use posterior SO2 emissions in 2010 as the prior emissions for all years in which we preform mass balance inversions. This facilitates decadal‐scale emission estimates to be more efficient than the 4D‐Var approach yet more accurate than the traditional mass balance approach.

4D‐Var Inversion in the Base Year

In our 4D‐Var inversions, we calculate emission scaling factors as the ratio of optimized emissions to prior emissions. The cost function, J, measures the sum of the error‐weighted squared departure of emission scaling factors from their prior estimates and the error‐weighted squared difference between model estimated and observed SCDs as where H is an observation operator that maps the species concentration vector () to the observation space, where only model concentrations within half an hour before and after the satellite overpassing time are sampled to reduce temporal representativeness error. SCD obs are OMI SO2 slant column densities. S is the observation error covariance matrix. Theoretically, it is the sum of the retrieval errors of the OMI SO2 columns and spatial representational errors. However, given the average of ~30% uncertainty in OMI SO2 retrievals and the smaller size of the model grid cells compared to the smearing length scale of SO2 (e.g., 260 km in the summer and 960 km in the winter using the definition of smearing length from Palmer (2003) and SO2 lifetimes from Lee et al. (2011)), we expect the measurement error to dominate and therefore only include this error in S . γ is a regularization parameter that is tuned to balance model error and prior constraints. are prior emission scaling factors, equal to 1. Ω is the domain (in time and space) where observations are available. S is the prior emission scaling factor error covariance matrix. We consider the monthly mean SO2 SCDs because the SO2 satellite retrievals have large random errors that cause a large amount of negative SO2 SCDs. As shown in Figure S19, the percentage of negative SO2 SCDs reduces when considering monthly means compared to individual values in all months using the BIRA retrieval and in most months using the NASA retrieval in China. In India, this percentage reduces in all months using the NASA retrieval and most months using the BIRA retrieval. Taking the seasonal mean of SO2 SCDs in India further reduces the percentage of SCDs that are negative from 21% (January 2010, BIRA) to 14% (winter mean of 2009 and 2010). While this percentage varies from month to month and region to region, the BIRA retrievals generally have fewer negative values than the NASA SP retrievals.

The observation error covariance matrix S is diagonal, because correlated data due to row anomalies have been filtered out, and we assume that the remaining observations are independent. Following Qu et al. (2017), we calculated the error for each observation as the sum of an absolute error and a relative error: where k is the index for observations in this grid cell and is the relative error, calculated as 110% for the NASA product (Krotkov et al., 2016) and as the sum (in quadrature) of slant column system error and slant column random error for the BIRA product. is the absolute error. For all retrieved observations, has a value of 0.252 DU for the NASA product (https://SO2.gsfc.nasa.gov/Documentation/OMSO2Readme_V130_20160616.pdf) and is calculated as the product of 0.35 DU and AMF for the BIRA product (Theys et al., 2015). At global and nested domains, all gridded BIRA AMFs are greater than 0.3, and 75% (at global domain, 67% at East Asia nested domain, 83% at North America nested domain) of them are greater than 0.4 in 2010, therefore ensuring that the absolute error is larger than 0.105 DU. Each diagonal element of the monthly mean observation error covariance matrix ( ) is calculated as where N is the total number of SO2 observations in the grid cell in a month.

For S , we assume that emissions errors are independent across grid cells, since major SO2 anthropogenic emissions come from power plants, which are estimated using unit‐specific parameters for each generation unit (Li, Zhang, et al., 2017). Therefore, emission errors in one grid cell are independent of emission errors in adjacent grid cells at ~50‐ and 250‐km resolution. We use a constant of 0.4 for all diagonal elements based on estimates from Li, Zhang, et al. (2017). Further adjustment of the weighting from errors in observations and prior estimates is fulfilled through regularization discussed below.

To optimize SO2 emissions, the quasi‐Newton Limited‐memory Broyden‐Fletcher‐Goldfarb‐Shanno Bounded optimization routine is applied to find the minimum of the cost function with respect to emission scaling factor. We calculate regularization parameter values for January 2010 in the East Asia, United States, and global domains using an L‐curve test and error minimization method (Figures S20–S22). The γ values in other months are estimated by scaling γ from January using the ratio of the number of SO2 observations in each month to that number in January.

In the global domain (Figure 5), adjustments to SO2 emissions in the 4D‐Var posteriors mainly occur in China, India, the Middle East, and the United States, where the magnitudes of SO2 column densities from anthropogenic sources are largest. China and India account for the largest anthropogenic sources of SO2, with annual emissions of 12.4 Tg/S for China and 4.36 Tg/S for India in 2010 in the prior inventory. SO2 emissions in other countries are much smaller (e.g., 0.42 Tg/S in the United States and 0.17 Tg/S in South Korea). We focus our following discussion on China, India, United States, and South Korea.

Figure 5

The 4D‐Var updates to SO2 emissions (posterior‐prior) when constrained using the BIRA retrieval products in 2010.

Compared to the prior emissions, SO2 emissions constrained by BIRA have smaller values in the Western United States and India in all months in 2010. Bottom‐up emissions in China are found to be generally overestimated in the southeast regions and underestimated in the North China Plain compared to BIRA posterior emissions. However, NASA posterior emissions are smaller in broader areas of East China (Figure S23). The increase of SO2 emissions in the North China Plain in the BIRA posterior emissions but not in NASA posterior emissions is consistent with the prior model SO2 column concentrations being smaller than the BIRA product yet larger than the NASA product over this region (see section 3). Although the model‐observation discrepancies are different in north and west China between the different retrieval products, posterior SO2 emissions are not (Figures 2 and S4) because OMI retrievals in these regions are below the detection limit. The signs of posterior SO2 emission changes agree well between the two posteriors in southern China and the United States due to small differences in model‐observation differences in these products at lower latitudes (as concluded in section 3). The normalized mean square differences between the global NASA and BIRA posterior emissions are smallest in August and September, corresponding to months when differences between SO2 retrievals are small. Comparisons with the nested East Asia simulations are shown in Figures S24 and S25, where emission changes at higher resolution show similar but more detailed spatial pattern in corresponding months than coarse resolution posteriors.

We next consider the seasonality of the 4D‐Var posterior emissions in China from the nested inversions, shown in Figure 6. For nested East Asia results, we also include NASA SP posteriors from equation (A8). The BIRA and NASA SP posterior emissions (AMF applied following equation (A9)) have seasonalities similar to the prior emissions, peaking in the winter due to residential and commercial heating. NASA posterior emissions using the AMF following equation (A8) are 12%–36% smaller than HTAP emissions and peak in summer, which is inconsistent with the current understanding of China energy usage. The annual budget of SO2 emissions in China (Table 2) reduces by 5% in the BIRA posterior and by 10% in the NASA posterior (AMF applied following equation (A9)). The difference between the posterior SO2 budget in China using the AMF calculation following equations (A8) and (A9) amounts to 15% in 2010, reflecting the large impact of different AMF treatment on emission estimates. The reduction in posterior emissions is consistent with Wang et al. (2016) and OMI‐HTAP estimates (Liu et al., 2018). In comparison, posterior SO2 emissions constrained by BIRA SO2 in Koukouli et al. (2018) increase from prior values of 13.2 Tg/S (the MEIC v1.2) to 16.3 Tg/S in China (15–55°N, 102–132°E). The different sign of the changes could be explained by the use of different chemical transport models (CHIMERE versus GEOS‐Chem) and inversion methods (domain‐wide mass balance versus 4D‐Var). It is beyond the scope of this paper to further investigate this aspect.

Posterior SO2 emissions in India have larger seasonal variation than the prior emissions, with a peak in summer, opposite to the peak in the prior inventory (bottom left panel of Figure 6). The differences in magnitude (Table 2) and seasonality between simulations using prior and posterior emissions are related to different seasonalities in the simulated versus measured SCDs in India (as shown in Figure S28) and a less accurate prior inventory due to limited and/or inaccurate statistics, as described in Liu et al. (2018). The 2010 SO2 annual budget in India reduces by 18% in the NASA posterior and by 26% in the BIRA posterior (Table 2).

Posterior SO2 emissions in the United States (bottom right panel of Figure 6) also peak in the summer due to cooling requirements (Abel et al., 2017; Nsanzineza et al., 2017), but the magnitude of the posterior emissions are 11%–27% (NASA) and 10%–28% (BIRA) smaller than the prior emissions in summer months. Posterior emissions mainly decrease in the Eastern United States and the annual budget decreases by 11% in both posteriors (Table 2), consistent with the difference between posterior and HTAP SO2 emissions in Liu et al. (2018). The larger HTAP emissions are likely due to omission of updates on fuel quality and technologies, discontinuous monitoring of SO2 sources, and misplacement of point sources over larger areas (Liu et al., 2018).

Posterior SO2 emissions in countries other than the above three barely change, and their seasonalities are less accurate, since the SO2 concentration magnitudes are close to the OMI noise level. For instance, changes of annual budget of posterior SO2 emissions in 2010 in Japan and South Korea are less than 2% compared to the prior. SO2 emissions in South Korea also have smaller seasonal variations compared to the other three (Figure 6).

Long‐Term Hybrid Posterior Emissions

The hybrid 4D‐Var/mass balance approach (Qu et al., 2017) is applied to derive posterior SO2 emissions for the nested East Asia and U.S. domains (2005–2012) and the global domain (2005–2017 for NASA, 2005–2016 for BIRA, given data availability). Consistent with the 4D‐Var inversion, a uniform error of 40% (Li, Zhang, et al., 2017) is ascribed to the prior emissions in China, and a uniform error (e.g.,Koukouli et al., 2018 ; Martin et al., 2003) of 30% in monthly average SO2 is applied to top‐down emissions in the mass balance calculation.

Emission Trends in Polluted Regions

China and India have the largest anthropogenic sources of SO2. SO2 emissions in other countries are much smaller (e.g., 0.17 Tg/S in South Korea and 0.42 Tg/S in United States); hence, we separate our trend discussion into polluted and clean regions, and focus in this section on posterior emissions in China and India. At 2° × 2.5°, posterior SO2 emissions in India increase by 35% (NASA) and 48% (BIRA) from 2005 to 2016 (Figure 7). These increasing trends are consistent with those from EDGAR v4.3.2 when also accounting for biomass burning emissions (magenta line). However, the growth rates in posterior emissions (13% for NASA and 22% for BIRA) are much smaller than the 76% increase in the EDGAR inventory from 2005 to 2012. The absolute magnitude of NASA (BIRA) posteriors are 14%–42% (16%–47%) lower than EDGAR (plus biomass burning) emissions.

Figure 7

Annual budget of SO2 emissions in (left) China and (right) India from 2005 to 2017. Coarse resolution refers to 2° × 2.5°; fine resolution refers to 0.5° × 0.667°.

Posterior SO2 emissions in China start to decrease from 2008 onward, which is in agreement with the MEIC inventory but not in the EDGAR v4.3.2 inventory. At 2° × 2.5°, the BIRA posterior emissions decreases by 33% from 2007 to 2016, which is larger than the decrease in the NASA posterior emissions of 24%. In comparison to the 27% decrease of OMI SO2 column densities from 2005 to 2014 (Liu et al., 2018), NASA and BIRA posterior SO2 emissions decrease by 16% and 23% and bottom‐up emissions (green line) decrease by 49%. The smaller magnitude in posterior emission reduction could be explained by prior constraints in 2010 and the omission of negative values in the hybrid approach, the latter of which has a larger impact on the trend of SO2 emissions in clean areas, as discussed in detail in section 4.2.2.

Posterior SO2 emissions from the global versus nested resolutions show consistent trends and similar growth rates, with negligible impact of resolution and meteorology on country‐level SO2 emissions. The annual budget of finer resolution posterior SO2 emissions is 5%–21% smaller than coarse resolution posterior emissions in India, and are −3%–15% smaller in China. By performing inversions at both resolutions, we benefit from the availability of the coarse‐resolution meteorology over longer time periods to exam long‐term trends, while also being able to perform a more careful comparison to in situ measurements using the fine‐resolution model.

The estimated changes of posterior SO2 emissions in this study are based on national total emissions. We note that these are smaller than estimates based on trends in 47 and 82 large SO2 point sources in India and China, respectively, in Li, McLinden, et al. (2017) and trends in 10 provinces with the highest SO2 concentrations of van der A et al. (2017). As shown in Tables S2 and S3, SO2 emission trends are different for sources of different magnitude. In China, larger sources generally have larger reductions from 2007 to 2016. In India, we find the largest increases in midsized SO2 sources. If we consider only grid cells that have emissions larger than 5 × 107 kg/S per grid (at 2° × 2.5°) in the prior inventory, our posterior SO2 growth rate would be −52% (NASA SP posterior emissions) and −65% (BIRA posterior emissions) for China from 2007 to 2016, which is close to the 75% decrease over the same period calculated in Li, McLinden, et al. (2017). Koukouli et al. (2018) also report SO2 emission trends at national scale; they find a decrease of 33% from 2005 to 2015 in China (constrained by BIRA observations), which is consistent with our decrease of 30% (BIRA posterior) over the same period. The annual budget of our posterior SO2 emissions in China is only 6% different when constrained by BIRA versus NASA products in 2010. Therefore, the nearly factor of two differences in SO2 emissions estimated in Wang et al. (2016) and Koukouli et al. (2018) are mostly related to the different treatment of SO2 retrievals, inversion methods, and chemical transport models.

Emission Trends in Cleaner Regions

Posterior SO2 emissions in less polluted areas are much harder to estimate, given the comparable magnitude of SO2 column concentrations to the OMI detection limit and the large noise in OMI SO2 retrievals, which results in a large number of negative SO2 retrievals. Negative SO2 retrievals are not used to scale emissions in the mass balance approach, since emissions cannot be negative. Here we use the United States as an example to analyze issues that arise when applying our hybrid inversion in cleaner areas, and we evaluate the threshold of SO2 concentrations above which this method can be applied for reasonable trend estimates. We benefit from the availability of surface SO2 measurements and relatively accurate bottom‐up SO2 emissions in the United States for evaluating the trends and magnitudes of our inversion results.

A significant challenge for inverse modeling in cleaner conditions is that a large number of OMI SO2 retrievals are negative, which can lead to negative country‐level monthly mean SO2 column densities. An example of this over the United States is shown in Figure S26. To further evaluate the capability of the hybrid method to recover trends, we separate our U.S. domain of study into East, Central, and the West United States. In Figure 8, left panel, the SO2 column density in the East United States shows a decreasing trend from 2005 to 2012, consistent with EPA Air Quality System surface measurements, National Emission Inventories emissions, and previous SO2 column (Krotkov et al., 2016) and top‐down studies (Fioletov et al., 2017). However, even in the largest SO2 concentration region (East United States), the trend of SO2 column is still impacted by the large amount of negative SO2 columns—as shown in the right panel of Figure 8, the decreasing trend in black lines when using both positive and negative retrievals to calculate the annual mean is not reflected in the magenta lines when only positive SO2 retrievals are used for the calculation. Therefore, the decrease of SO2 could not be reflected in our hybrid posterior emissions from both NASA and BIRA products, given that only positive OMI SO2 columns are used to recover the trend of emission.

Figure 8

Trends of SO2 column concentrations (left) including both positive and negative retrievals in U.S. regions and (right) over the East United States using only positive values or both positive and negative values.

We also find that the trend of SO2 emissions in cleaner area (e.g., United States, South Korea, and Japan; shown in Figure S27) are less consistent between NASA SP and BIRA posteriors (R: 0.72, 0.18, and 0.46, respectively) than those in China and India (R: 0.86 and 0.75). Evaluations with surface measurement in sections 5 and S3 further demonstrate the greater robustness of posterior emissions in China and India than in less polluted countries, that is, South Korea. This is mainly hindered by large amounts of negative SO2 retrievals over small sources and the comparable magnitude of SO2 column densities over clean areas with the noise level of the retrievals. The percentage of grid cells that are partly over the ocean is also larger in South Korea and Japan due to their smaller areas, which dilutes our calculation of national mean SO2 concentrations for these regions. We therefore only focus on applying the hybrid method to estimate trend of nation‐wide SO2 emissions in China and India in this work. Other emission trend estimates focusing on large point sources instead of nation‐wide emissions using a plume model approach (e.g., Fioletov et al., 2016; Fioletov et al., 2017; Liu et al., 2018; ; McLinden et al., 2016) take long‐term averages of SO2 observations and are less hindered by negative SO2 retrievals. Still, SO2 emissions below 15 kt S/year (determined by the detection limit of OMI SO2 column densities, same for all methods) could hardly be detected (Fioletov et al., 2016).

Evaluation With In Situ Measurements

To evaluate the magnitude and seasonality of posterior SO2 emissions, we use them in GEOS‐Chem to simulate SO2 concentrations to compare with in situ measurement. All posterior simulations in this section are from the hybrid (rather than 4D‐Var) approach.

Surface SO2 Concentration

Evaluation of spatial and temporal variability of SO2 concentrations in one year is representative of the model performance, given the similar model‐observation discrepancies at the same month of each year and the similar seasonality in simulated and retrieved SO2 column densities (Figure S28). Therefore, we simply choose the year where bottom‐up emissions (HTAP) are available (2010), so that we can make a direct evaluation of the model performance when using bottom‐up and top‐down emissions. Furthermore, inconsistencies in the availability of the surface measurements from different locations in different years preclude analysis of interannual variability in India. Table 3 shows statistics when comparing the spatial variability of GEOS‐Chem surface layer SO2 concentration with surface measurements in 2010. In China, the magnitude of the NMB and NMSE has decreased compared to the prior using the BIRA posterior but has increased using the NASA SP posterior; posterior simulations constrained using AMFs as described by equation (A8) have a larger error and bias and worse correlation than simulations using the prior emissions and simulations using the NASA SP and BIRA posterior emissions when compared with surface measurements. Spatially, posterior simulations show reduced SO2 concentrations in North China Plain, Yangtze River Delta, Sichuan, Chongqing, Hubei, and Hunan, which better match surface measurements, but the reduced concentrations in the Pearl River Delta in the posterior deviate more from the measurements. Simulated SO2 concentrations using bottom‐up emissions at monitoring sites closer to emission sources, that is, in provincial capitals (Table 3) have a larger bias than the national mean. This suggests that the HTAP inventory is likely overestimating SO2 emissions in cities but underestimating emissions in rural areas.

In South Korea, both NASA SP and BIRA posterior simulations have a better spatial correlation but a larger negative bias than prior simulations when compared to surface SO2 measurements. The posterior annual mean surface SO2 concentrations also deviate more from the measurements compared to the prior simulations, which are initially biased low. Posterior SO2 simulations in the 10 large cities in South Korea also deviate more from the measurements. Further consideration of causes due to vertical distributions is discussed in section 5.2.

In the United States, simulated surface concentrations are larger than measured concentrations. Posterior simulations have a reduced error and bias. However, in the 10 large cities, prior simulations have a smaller bias but larger error than posterior simulations.

We also evaluate the interannual variation of SO2 concentrations in China. In Figure 9, SO2 measurements decrease by 28% from 2005 to 2012, whereas simulations using fixed anthropogenic emissions (HTAP 2010) decrease by 0.2%, reflecting small impacts from meteorology. Biomass burning SO2 emissions are less than 0.2% of total SO2 emissions in China, and are not affecting the trends here. Simulated SO2 concentrations at 0.5° × 0.667° resolution decrease by 10% (NASA SP posterior) and 14% (BIRA posterior) from 2005 to 2012, and by 5% (NASA SP posterior) and 14% (BIRA posterior) from 2005 to 2010. Measured SO2 concentrations decrease by 28% from 2005 to 2012 and 17% from 2005 to 2010. The greater inconsistency between posterior simulations and measurements from 2010 to 2012 may be caused by the reduced number of OMI observations due to the presence of the row anomaly in later years. The generally smaller trends in the simulated concentrations using posterior emissions compared to the measurements are related to the underestimation of SO2 emission trends as discussed in section 4.2.1. The surface measurements decrease from 2005 onward rather than from 2008 onward as with the posterior emissions, possibly due to inconsistent measurement methods throughout the period, measurement error, or different trends of SO2 concentrations at the surface versus the whole column in the real atmosphere. The correlation of surface SO2 growth rate between measurements and simulations is highest for BIRA posterior (R = 0.83), and lowest for prior (R = 0.34), with the correlation for the SP posterior in between (R = 0.68).

Figure 9

Percent changes relative to 2010 in SO2 concentrations in China from surface measurements (from the 272 sites with data in every year from 2005 to 2012) and co‐located estimates from prior and posterior GEOS‐Chem simulations.

The number of available monitoring sites in our studied period is rather limited in India, so we focus only on the magnitude and seasonality of SO2 concentrations in 2010. Unlike in China and South Korea, surface SO2 simulations in India using bottom‐up emissions are 94% larger than the measurements (Table 3). The posterior monthly simulations are smaller by 8% (NASA SP posterior) and 25% (BIRA posterior) than the surface measurements. As shown in Table 4, the monthly variability of domain‐average posterior simulations using NASA posterior emissions have a reduced error and a smaller bias compared to the BIRA posterior emissions. The better monthly correlation between the prior SO2 simulation and surface measurement suggests possible errors in the seasonality of posterior SO2 emissions related to uncertainties in OMI retrievals and GEOS‐Chem model errors.

Table 4

Statistics of Simulated SO2 Monthly Variability When Compared With Surface Measurements in China, South Korea, United States, and India

	Emissions	Prior (HTAP)	NASA SP posterior	BIRA posterior
China	R	0.94	0.92	0.92
	NMB	−0.36	−0.14	−0.31
	NMSE	0.27	0.04	0.17
South Korea	R	0.68	0.65	0.71
	NMB	−0.51	−0.45	−0.50
	NMSE	0.59	0.41	0.56
India	R	0.61	0.50	0.21
	NMB	0.94	−0.08	−0.25
	NMSE	0.33	0.19	0.48
United States	R	0.80	0.85	0.83
	NMB	0.62	0.15	0.28
	NMSE	0.26	0.07	0.12

Statistics are calculated from the monthly domain averages of surface measurements and GEOS‐Chem simulations throughout 2010.

In the United States, surface concentrations in the posterior reduce the NMSE relative to Air Quality System measurements by 54% (BIRA posterior) and 73% (NASA SP posterior) when compared to the prior, as shown in Table 4. The NMB also reduces by 55% (BIRA) and 76% (NASA) compared to the prior. Interannual variations of SO2 surface concentrations are not evaluated here due to a lack of robustness of trends of the hybrid posterior simulations in clean regions, as discussed in section 4.2.2.

Evaluation of SO2 Vertical Profiles

In South Korea, the prior simulations of surface SO2 concentrations are smaller than the measurement, but simulations using OMI posteriors are even smaller than the prior. The surface SO2 amounts and distributions can exhibit highly nonlinear relationships with respect to the emissions and a stronger dependence on the meteorological conditions (Calkins et al., 2016). To investigate SO2 vertical profiles in different regions and their impact on the comparison of simulated surface SO2 concentrations with measurements, we compare the SO2 vertical profile in GEOS‐Chem simulated using the prior emissions with KORUS‐AQ DC‐8 aircraft measurements from April to June 2016 (left panel of Figure 10). Due to the limitation of this version of the GEOS‐Chem adjoint model to use nested‐resolution GEOS‐5 meteorology that is only available up to 2013, GEOS‐Chem simulations in 2016 are performed at 2° × 2.5° resolution using MERRA‐2 and HTAP 2010 anthropogenic emissions. However, we tested the impact of model resolution on the simulated SO2 vertical profile by comparing 0.5° × 0.667° (GEOS‐5) and 2° × 2.5° (MERRA‐2) simulations in 2010 sampled at the corresponding time and location of KORUS‐AQ DC‐8 measurements in 2016. As shown in the green (0.5° × 0.667°) and red (2° × 2.5°) lines in Figure 10, the horizontal resolution has little impact on the SO2 vertical profile. All simulated SO2 concentrations are smaller than aircraft measurements at the surface layer, whereas model SO2 concentration at ~1.5–3 km are larger than in situ measurements possibly due to higher wind speed in the model during KORUS‐AQ, lower PBL in the real atmosphere at night, or potential errors in SOx partitioning in the model since differences in the vertical profile are opposite to that of aerosol extinction noted in van Donkelaar et al. (2013).

Figure 10

Comparison of GEOS‐Chem SO2 vertical profiles with (left) KORUS‐AQ DC‐8 aircraft measurements from 26 April to 18 June of 2016 and (right) DISCOVER‐AQ aircraft measurements over Virginia and Maryland in the United States from 8 June to 25 August 2011. The horizontal bars show the 25% and 75% quartiles of the measurements averaged at each height, and the year and horizontal resolution of the GEOS‐Chem simulations are indicated in the legends.

The larger model concentrations around 1.5–3 km potentially contribute to larger total SO2 column density in GEOS‐Chem than in OMI retrievals. This larger simulated column leads to a decrease of SO2 emissions, even though simulated surface SO2 concentrations are still smaller than surface measurements. This difference in SO2 vertical profile may explain why the bias in the posterior SO2 simulations has increased with respect to surface measurements in South Korea.

A similar effect is not seen for China, India, and United States partly due to the larger prior surface SO2 concentrations compared to the measurements. We also compare GEOS‐Chem SO2 vertical profile with DICOVER‐AQ measurements over Virginia and Maryland in 2011 (right panel of Figure 10). Prior GEOS‐Chem simulations are larger than measurements near the surface but smaller at higher altitude, opposite to the case in South Korea. Therefore, the decrease of SO2 emissions in United States reduces the bias and error compared to the measurements (Tables 3 and 4).

Discussion and Conclusions

We applied a hybrid 4D‐Var/mass balance approach to estimate SO2 emissions constrained by OMI observations. Given the differences in three recent OMI SO2 products, we also compared similarities and differences in their SO2 column densities and evaluated how they affect emission estimates. Compared to the prior in 2010, posterior SO2 emissions show consistent reductions in East United States, India, and southeast China when using either NASA SP or BIRA products. However, changes relative to the prior are opposite in sign in North China Plain when using these different retrievals. This discrepancy stems from differences between OMI SO2 retrieval, which are more consistent across products at low latitude and in the summer.

Even though the spatial distribution of OMI SO2 column densities are different across retrievals, satellite‐based national‐scale posterior SO2 emission trends are similar when constrained using either the NASA and BIRA products, or when constrained using a model at different resolutions. In China, we find decreasing SO2 emissions since 2008 in the posteriors (by 24% in NASA SP and 48% in BIRA by 2016), which is consistent with the decrease in the MEIC anthropogenic emissions and the timing of implementation of SO2 emission control measures, but opposite the trend in the EDGAR v4.3.2 inventory. In India, both SO2 posteriors and anthropogenic emission from EDGAR v4.3.2 increase from 2005 to 2017 (by 35% in NASA SP and 48% in BIRA till 2016), but the growth rate in the posterior emissions is significantly smaller. Indian emissions are expected to continue to increase due to electricity and industrial demand (Li, McLinden, et al., 2017; Venkataraman et al., 2018). Still, as of 2017, the annual budget of posterior SO2 emissions in China is approximately 2.2 times that of India. Annual budgets of posterior SO2 emissions are mostly lower than bottom‐up estimates in China and India, which is consistent with the factor of 6 high bias for power plant SO2 emissions in bottom‐up inventories compared to top‐down estimates from Liu et al. (2018). These are related to the distribution of SO2 sources over too few point sources in bottom‐up inventories (Liu et al., 2018), lower SO2 column densities in OMI retrievals compared to simulations, and small sources that are not detectable in OMI retrievals when the signal is comparable to the noise. The constraints from the prior emissions in the mass balance method (more details in Qu et al. (2017) and Qu et al. (2019) further limit the posterior emissions from deviating too much beyond the prior emissions, which were for 2010.

The performance of the hybrid method is degraded by the presence of negative SO2 retrievals from OMI, which are especially frequent over areas with smaller SO2 concentrations. Given this, posterior seasonality and trends are likely to be most informative in China and India, where SO2 loadings are high. The trends of posterior SO2 emissions in the other countries are less consistent between NASA and BIRA posteriors. The hybrid inversion performance could potentially be improved in future studies by assimilating observations with higher spatial resolution, larger signal‐to‐noise ratios, and a lower detection limit from TROPOMI (Theys et al., 2017) and upcoming geostationary satellites (e.g., Geo‐stationary Environmental Monitoring Spectrometer (GEMS), Tropospheric Emissions: Monitoring of Pollution (TEMPO), and Sentinel‐4).

In the emissions we reported, we excluded contributions from volcanoes but do not separate biomass burning and anthropogenic emissions. However, in the countries where the posterior emissions are robust (China and India), biomass burning emissions amount to less than 0.3% of anthropogenic emissions and their contribution to our estimate is negligible. We also calculated the percentage of monthly posterior SO2 emissions at grid cells where emissions change more than 5% compared to the total monthly posterior SO2 emissions in China and India. In 2010, this percentage ranges from 65% (18%) to 82% (60%) in the BIRA (NASA) posterior emissions in China, and from 51% (14%) to 75% (66%) in the BIRA (NASA) posterior emissions in India. Since a lack of change in emissions can also be indicative of correct prior emissions, we expect the percentage of SO2 emissions that OMI can detect to be higher than these numbers.

Compared with surface measurements, the NASA SP posterior has a smaller NMB, a smaller NMSE, and a larger correlation than the BIRA posterior in United States in terms of both spatial (Table 3) and seasonal (Table 4) variability, suggesting more accurate SO2 column and emissions in the United States from the NASA SP. In China, South Korea, and India, the NASA SP posterior has smaller NMB and NMSE than the BIRA posterior in terms of seasonality, but statistics of the spatial variability have worsened. The annual growth rate of simulated surface SO2 concentration from the BIRA posterior also has a better correlation and consistency in growth rate with surface measurements.

The national mean of posterior surface SO2 concentrations do not improve compared with in situ measurements in China and South Korea. However, at provincial capitals in East China the annual mean of posterior simulations improve compared with measurements. A possible reason for the lack of improvement in simulated surface SO2 concentrations comes from the SO2 vertical profile in GEOS‐Chem, which has lower concentrations than aircraft measurement at the surface but higher at 1–4 km, leading to larger SO2 columns but smaller surface SO2 concentrations. Further investigation of sulfate partitioning, vertical mixing and wind speed in the model, and assimilation of SO2 vertical profiles in the real atmosphere could help correct GEOS‐Chem SO2 vertical distribution for more accurate emission estimates. Additionally, errors in emissions of other species and simulated O3 and OH concentrations could lead to biased adjustments to SO2 emissions. The global mean OH concentration in GEOS‐Chem ranges from 10 × 1015 to 13 × 1015 molec/cm3 from version 8 to version 12 (http://wiki.seas.harvard.edu/geos-chem/index.php/Mean_OH_concentration). Our subsequent study incorporates multiple species observations to synergistically change O3 and OH concentration, and thus provide more accurate constraints on SO2 emission estimates (Qu et al., 2019).

Besides uncertainties in model simulation, several factors in satellite retrievals also contribute to the uncertainties of SO2 emission estimates. First, different SO2 column densities from NASA and BIRA products lead to different signs of SO2 emission changes at locations where SO2 loadings are high. Second, posterior emissions are sensitive to the different treatment of the AMF, so it is important that the next version of NASA operational product includes additional scattering weight and AMF information. Still, it is worthwhile to estimate top‐down emissions using a constant AMF of 0.36 for comparisons with existing trend and emission studies using this approach (e.g., Krotkov et al., 2016; Li et al., 2017; Li et al., 2018). The consistent estimates and evaluations of three SO2 products and two sets of AMF applications (equations (A8) and (A9)) in the NASA SP in this study provides a bridge between existing studies using the current NASA SP to the ones using products with retrieval‐specific scattering weights (BIRA and NASA prototype). The assumed center of mass altitude in the NASA product also contributes to uncertainties in the Jacobian and the AMF calculation and thus SO2 column densities, as the sensitivity of SO2 absorption strongly depends on this altitude.

We demonstrate in this study the benefit of using multiple satellite retrievals and multiscale in situ measurements to evaluate the robustness of posterior emissions and to identify potential issues in model processes and satellite retrievals. These analyses are especially important for posterior evaluations for an inversion system that does not account for biases in the observations and for which quantification of posterior uncertainties is nontrivial (Bousserez et al., 2015). Further application to emission estimates of other O3 and PM2.5 precursor gases could help improve understanding of pollutant sources and the performance of chemical transport model in predicting pollutant concentration and spatial distribution.

Supporting Information S1

Click here for additional data file.