Hybrid Mass Balance/4D-Var Joint Inversion of NO x and SO2 Emissions in East Asia.

Accurate estimates of NO x and SO2 emissions are important for air quality modeling and management. To incorporate chemical interactions of the two species in emission estimates, we develop a joint hybrid inversion framework to estimate their emissions in China and India (2005-2012). Pseudo observation tests and posterior evaluation with surface measurements demonstrate that joint assimilation of SO2 and NO2 can provide more accurate constraints on emissions than single-species inversions. This occurs through synergistic change of O3 and OH concentrations, particularly in conditions where satellite retrievals of the species being optimized have large uncertainties. The percentage changes of joint posterior emissions from the single-species posterior emissions go up to 242% at grid scales, although the national average of monthly emissions, seasonality, and interannual variations are similar. In China and India, the annual budget of joint posterior SO2 emissions is lower, but joint NO x posterior emissions are higher, because NO x emissions increase to increase SO2 concentration and better match Ozone Monitoring Instrument SO2 observations in high-NO x regions. Joint SO2 posterior emissions decrease by 16.5% from 2008 to 2012, while NO x posterior emissions increase by 24.9% from 2005 to 2011 in China-trends which are consistent with the MEIC inventory. Joint NO x and SO2 posterior emissions in India increase by 15.9% and 19.2% from 2005 to 2012, smaller than the 59.9% and 76.2% growth rate using anthropogenic emissions from EDGARv4.3.2. This work shows the benefit and limitation of joint assimilation in emission estimates and provides an efficient framework to perform the inversion. ©2019. The Authors.

Introduction

Nitrogen oxides (NO = NO + NO2) and sulfur dioxide (SO2) are precursor gases for PM2.5 (Seinfeld & Pandis, 2006). NO also leads to the formation of ozone (Seinfeld & Pandis, 2006). These pollutants decrease visibility (Haagen‐Smit, 1952), harm respiratory system (Burnett et al., 1995; Schwartz et al., 2012), and damage ecosystems. Quantification of long‐term NO and SO2 emission are important for evaluation of air quality regulations (e.g., de Foy et al., 2016). Their emissions are also used in chemical transport models to study the formation of PM2.5 (e.g., Tao et al., 2017; Zhang et al., 2015) and causes of haze (e.g., Gao et al., 2016; Zhang et al., 2016).

Bottom‐up NO and SO2 emission inventories from Streets et al. (2003), Zhang et al. (2009), Regional Emission inventory in ASia (REAS) (Kurokawa et al., 2013; Ohara et al., 2007), and MIX (Li, Zhang, et al., 2017) have been extensively used for air quality and human health studies in Asia (e.g., Gao et al., 2018; Kanaya et al., 2017; Uno et al., 2003; Wang et al., 2016). However, uncertainties in anthropogenic NO and SO2 emissions from these inventories are up to 49% and 35% (Kurokawa et al., 2013; Lu et al., 2011; Zhang et al., 2009; Zhao et al., 2011). The unified compilation framework used in Streets et al. (2003) and REAS (Kurokawa et al., 2013; Ohara et al., 2007) introduces uncertainties due to lack of knowledge of emission factors, chemical profiles, spatial proxies, temporal profiles, etc., whereas the mosaic approach, which harmonizes multiple emission inventories at different regions into one product, brings in errors from inconsistencies among different data sets (Janssens‐Maenhout et al., 2015; Li, Zhang, et al., 2017; Zhang et al., 2009). Bottom‐up emission inventories may also take a long time to compile, often with more than 1‐year lag in time until they become available for use.

Alternatively, satellite observations provide an additional means of assessing emissions that is often spatially comprehensive and potentially in near real time; as such, they have recently been used to study the trend of NO2 and SO2 column densities and emissions. Decreases of NO2 (de Foy et al., 2016; Duncan et al., 2016; Krotkov et al., 2016; Liu et al., 2016) and SO2 (Calkins et al., 2016; Krotkov et al., 2016; van der A et al., 2017) column concentration have been detected from Ozone Monitoring Instrument (OMI) retrievals over China since 2011 and 2008. However, trends of species column densities are not the same as their emissions trends due to the impact of chemistry, meteorology and transport (Qu et al., 2017). To further estimate emissions in East Asia, physical models have been combined with satellite observations through inverse modeling techniques. For instance, the four‐dimensional variational (4D‐Var) method (e.g., Elbern et al., 1997, 2000; Henze et al., 2007, 2009) has been compared with the mass balance approach (e.g., Boersma et al., 2008; Martin, Jacob, Chance, et al., 2003; Toenges‐Schuller et al., 2006) in Qu et al. (2017), and a hybrid 4D‐Var/mass balance method has been developed to derive top‐down NO emissions for China from 2005 to 2012. Extended Kalman filter (e.g., Ding, van der A, et al., 2017; Mijling et al., 2013) and ensemble Kalman filter (e.g., Miyazaki et al., 2017) are employed to estimate multiyear NO emissions in China. 4D‐Var (e.g., Qu et al., 2019; Wang et al., 2016; Xu et al., 2013), mass balance (e.g., Koukouli et al., 2018; Lee et al., 2011), and plume methods (e.g., Fioletov et al., 2013; McLinden et al., 2016) have also been applied to derive SO2 emissions from OMI SO2 and Moderate Resolution Imaging Spectroradiometer aerosol optical depth (AOD) observations.

Though posterior emissions estimated based on single‐species observations provide insights into the trend of air pollutants, chemical interactions among atmospheric species are often overlooked in the assimilation system, which can lead to errors in the derived emissions. To improve NO model simulations, Chai et al. (2006) assimilated O3, NO, NO2, HNO3, PAN, and RNO3 aircraft observations and found significantly better agreement of simulated NO with measurements than when only assimilating NO observations. Hamer et al. (2015) evaluated the ability of a photochemical box model to predict ozone and found assimilating multiple ozone precursors using 4D‐Var decreases the uncertainty of ozone forecasting than the prior. Miyazaki et al. (2017) applied multispecies assimilation of NO2, O3, CO, and HNO3 to estimate NO emissions using an ensemble Kalman filter and over India obtained posterior emissions 10% smaller and a posterior emissions trend 15% smaller than the NO2‐only assimilation. The improvement of NO emissions from the multispecies inversion is confirmed by the better agreement of simulated and observed O3 concentration (Miyazaki et al., 2017; Miyazaki & Eskes, 2013).

In this work, we extended a recently developed hybrid mass balance/4D‐Var method (Qu et al., 2017) to estimate long‐term (2005–2012) NO and SO2 emissions in East Asia simultaneously using joint NO2 and SO2 satellite constraints. Simultaneous assimilation of these two species helps reduce biases caused by interactions of NO and SO2 through ozone chemistry and aerosol thermodynamics. Since observations of both species are from OMI, the measurement techniques, resolution, observation operator, available time period, and locations are more consistent. In theory, better emission estimates of these two aerosol precursor gases can improve model simulation of PM2.5 concentration and AOD.

We first update the 4D‐Var data assimilation framework and develop a joint mass balance method to estimate emissions from multiple species observations. The performance of this joint 4D‐Var and joint mass balance approach are evaluated through pseudo observation tests over China in section 3. These joint inversion methods are then combined (hybrid method) to generate eight years of NO and SO2 emissions for China and India in section 4. In section 5, we evaluate posterior simulations with surface measurements.

Model and Observations

GEOS‐Chem and Its Adjoint Model

The GEOS‐Chem adjoint model (Henze et al., 2007) v35f is used to perform NO and SO2 emission joint inversions. We use Goddard Earth Observing System (GEOS‐5) reanalysis meteorology field from National Aeronautics and Space Administration (NASA) Global Modeling and Assimilation Office from 2005 to 2012 (Bey et al., 2001). The GEOS‐5 meteorological data have a native horizontal resolution of 0.5° × 0.667°, 72 vertical layers, and a temporal resolution of 3 or 6 hr. GEOS‐Chem nested‐grid (70–150°E, 0–50°N) simulations are performed at 0.5° × 0.667° horizontal resolution and 47 vertical layers, with 3‐hourly boundary conditions generated for the first three grid cells at each of the four sides from global 4° × 5° simulations.

A detailed O‐NO‐hydrocarbon chemical mechanism (Bey et al., 2001) is included in the GEOS‐Chem simulation. The sulfur cycle simulation (emission, chemistry, advection, convection, diffusion, dry deposition, and wet deposition) is implemented by Park et al. (2004) based on the Goddard Global Ozone Chemistry Aerosol Radiation and Transport model (Chin et al., 2000). The gas and particle phase partitioning of nitric acid (HNO3) and ammonia (NH3) is calculated through the RPMARES aerosol thermodynamics scheme from Park et al. (2004). Heterogeneous reaction of NO2, NO3 (Martin et al., 2003), and N2O5 (Evans & Jacob, 2005) to form HNO3 is included. Dry deposition in GEOS‐Chem is computed using a resistance‐in‐series model (Wang et al., 1998; Wesely, 1989). Wet deposition is described in Liu et al. (2001).

Anthropogenic emissions of NO, SO2, NH3, CO, non‐methane volatile organic compounds (NMVOCs), and primary aerosol from HTAP 2010 inventory version 2 (Janssens‐Maenhout et al., 2015) are used to drive all of our prior simulations from 2005 to 2012. We use nonanthropogenic emissions corresponding to each simulation year. Three hourly biomass burning emissions are from fourth‐generation global fire emissions database (Giglio et al., 2013). Lightning NO emissions are calculated using the cloud top height parameterization (Price & Rind, 1992), vertical distribution profile (Pickering et al., 1998), local redistribution (Sauvage et al., 2007), and satellite constraints (Murray et al., 2012). Yienger and Levy's algorithm (Yienger & Levy, 1995) and soil canopy reduction factors (Wang et al., 1998) are used for deriving soil NO emissions. SO2 emissions from oxidation of dimethyl sulfide follow Park et al. (2004). Treatment of volcanic SO2 emissions follows Qu et al. (2019).

OMI NO2 and SO2 Observations

OMI is a nadir spectrometer measuring visible and ultraviolet (264–504 nm) solar backscatter radiation onboard the NASA Aura satellite. It has a Sun‐synchronous ascending orbit overpassing the equator at approximately 13:45 local time. OMI observations have a daily global coverage with a footprint of 13 km along track and 24 km across track. Aerosol and gaseous column retrievals are available since October 2004.

There are three recent operational OMI NO2 retrievals available over East Asia domain. The NASA standard product (Krotkov et al., 2017), Quality Assurance for Essential Climate Variables product (Boersma et al., 2018), and Peking University OMI NO2 (POMINO) product (cover China and part of India; Liu et al., 2019). In this study, we assimilate the NASA standard product OMNO2 (Level 2, Version 3) tropospheric NO2 slant column density (Krotkov et al., 2017). It is retrieved based on measurement spectrum of 405–465 nm (Boersma et al., 2011; Bucsela et al., 2013). Errors of the retrieved tropospheric NO2 column come from the total slant column density (SCD), separation of the stratosphere and troposphere column concentrations, and calculation of the tropospheric air mass factor (AMF). We screen the data by the criteria of positive tropospheric column, cloud fraction <0.2, solar zenith angle <75°, and viewing zenith angle <65°. The bias introduced by only using positive NO2 column density is small, as described in Qu et al. (2017). We further filter the observations using retrieval quality flags and exclude data with row anomalies (http://projects.knmi.nl/omi/research/product/rowanomaly-background.php).

We use the Royal Belgian Institute for Space Aeronomy (BIRA) SO2 Level 2 product in this study (Theys et al., 2015). The choice of this product comes from the availability of scattering weight and better consistency in SO2 trend evaluated with surface measurement over East Asia (Qu et al., 2019). The BIRA product is retrieved using the Differential Optical Absorption Spectroscopy technique (Platt & Stutz, 2008) for the 312–326 nm wavelengths. Given the focus of this study on anthropogenic SO2 emissions, we only employ tropospheric SO2 VCDs smaller than 5 Dobson unit (DU) to exclude cases of transient volcanic SO2, following Lu et al. (2013). Only footprints with cloud fraction less than 0.2 are used for the inversion. Data affected by row anomaly and with SZA larger than 65° are excluded through screening of quality flags.

In Situ Measurements

To evaluate simulated surface layer pollutant concentration over China, we use daily SO2 and NO2 surface measurements from the China National Environmental Monitoring Center from 2005 to 2012. There are 1,118 sites, all of which are managed by local governments in China, with different vendors and calibrations. However, not all of these sites have NO2 or SO2 measurement at every month during the studied period. For the NO2 measurement, NO compounds are converted to NO and measured by chemiluminescence analyzer. Therefore, the reported NO2 concentrations are generally biased high, yet the portion of converted NOz is variable at different locations and seasons.

Surface measurements of NO and SO2 in India are from Central Pollution Control Board (CPCB) performed by India's Ministry of Environment, Forest and Climate Change. The measurement methods at different sites can be different depending on the laboratory conducting the measurements, but interlaboratory comparison are conducted using Ring Test Facility (http://www.cpcb.nic.in/newitems/7.pdf). We use hourly averaged measurements in 2010 from ~20 sites. Since most of these sites do not have long‐term measurements covering our entire period of study, we choose not to evaluate the trend of surface SO2 concentration with measurements because that may only reflect spatial differences in SO2 concentration.

Joint Inversion Methods

4D‐Var Joint Inversion

Updates of 4D‐Var Framework for Joint Inversion

We assimilate NO2 and SO2 SCDs from OMI retrievals following Qu et al. (2017) and Qu et al. (2019), where scattering weights are applied to convert both OMI and simulated vertical column densities to slant column densities. We use an assimilation window of one month for both NO2 and SO2 observations. The 4D‐Var framework is updated to simultaneously assimilate NO2 and SO2 observations and constrain both surface NO emissions, E (i,j), and surface SO2 emissions, E (i,j) at each grid cell with longitude index i and latitude index j. Any errors in lightning NO are compensated by adjusting surface NO emissions, which could introduce an additional source of uncertainty in interpreting posterior emissions trends as indicative of anthropogenic activity or in comparison to other bottom‐up anthropogenic inventories. The linear emission scaling factors for NO and SO2 are defined as and where E (i,j) and E (i,j) are prior emissions for NO and SO2. The cost function

measures the error weighted departure of NO and SO2 emission scaling factors ( and respectively) from their prior estimates ( and , all elements equal to 1) known as “parameter error”, and the sum of uncertainty weighted squared error between model and observations (“prediction error”) over time and the simulation domain. is an observation operator that maps species concentrations of NO2 ( ) and SO2 ( ) to observation space to be comparable with satellite NO2 (SCD ) and SO2 (SCD ) slant columns. To help mitigate the presence of unrealistic negative SO2 column densities and avoid introducing bias into the system by simply filtering these out, we averaged OMI () and GEOS‐Chem ( ) SO2 SCDs overpassing each grid cell over each month following Qu et al. (2019) and use these monthly averages ( and ) in the cost function. Calculations of monthly mean observation error covariance matrices also follow Qu et al. (2019). Monthly averaging is not needed for NO2 because less than 5% of the NO2 retrievals are negative, and these are simply excluded from the assimilation. S and S are the error covariance matrices of NO and SO2 emission scaling factors, respectively. A constant fractional error of 0.4 is used for all diagonal elements in these two matrices based on Li, Zhang, et al. (2017) and is further adjusted as described in the next paragraph. Use of the same uncertainty in all grid cells is an approximation owing to lack of readily available detailed emission uncertainty information. Imperfect error estimates could cause uncertainties in the magnitude of the emission changes relative to the prior. Major sources of anthropogenic SO2 emissions are power plants, which are estimated for each generation unit using unit specific parameters (Li, Zhang, et al., 2017). Therefore, emission errors from these point sources are not spatially correlated at the 0.5° × 0.667° resolution. NO emissions mainly come from transportation, inventories of which use county‐level emission factors and provincial‐level technology distribution, whose errors are spatially correlated. We therefore assume independent emission errors for SO2 and exponentially decaying error correlation for NO with a decay distance of 150 km. S and S are diagonal error covariance matrices of NO2 and SO2 observations, respectively. Each of their diagonal elements is the sum of a relative error of the corresponding footprint provided in the retrievals and an absolute error of the product of 0.35 DU and AMF for SO2 and 1 × 1015 molecules/cm2 for NO2. More details of the error covariance matrices are discussed in Qu et al. (2017) for NO2 and Qu et al. (2019) for SO2. Ω is the domain (in time and space) where observations and model simulations are available.

The linear addition of two species parameter and prediction error here assumes the emissions and observations errors are uncorrelated. In order to weight NO2 and SO2 prediction error equally in the cost function, we scale the SO2 prediction error by α, which is the ratio of the number of NO2 and effective SO2 observations (i.e., number of grid cell that have SO2 observations). The impact of α is further discussed in section 4.1. A regularization parameter, γ , is used to adjust the magnitude of the parameter errors, and its value is determined using an L‐curve test (Hansen, 1999) and minimization of the total error (Henze et al., 2009). Based on test results conducted for January 2010 (Figure S1), we choose a γ value of 50. For other months in 2010, γ values are adjusted based on the effective number of observations (α × number of SO2 observations + number of NO2 observations). These values are shown in Table S1.

Evaluation of Joint 4D‐Var Inversion Using Pseudo Observation Test

We design an inverse problem with known emission solution to evaluate the performance of the joint 4D‐Var inversion system. Bottom‐up emission inventories described in section 2.1 are used to generate hourly NO2 and SO2 column densities within 30 min of OMI overpass time using GEOS‐Chem. Random noises with mean of 1 and standard deviation of 0.2 (calculated based on monthly mean OMI NO2 error for January 2010) and mean of 1 and standard deviation of 0.3 (based on monthly mean OMI SO2 error for January 2010) are applied to each NO2 and SO2 pseudo observations.

We design two cases to evaluate the inversion performance when bottom‐up emissions are biased high or low: (a) both SO2 and NO emissions across the entire model domain are scaled to 0.5 times the true emissions with normally distributed random noise N(1,0.1) (only using values between 0.6 and 1.4) applied to SO2 and spatially correlated emission noise N(1,0.1) (only using values between 0.6 and 1.4) to NO emissions (referred to as “σ = 0.5” case); (b) both SO2 and NO emissions are scaled to 1.5 times the true emissions and the same noise as in the σ = 0.5 case are applied (“σ = 1.5” case). Spatially correlated noise is only applied to NO emissions but not SO2 emissions in order to mimic emission error correlations in the real case, as described in section 3.1.1. We use a standard deviation of 0.1, instead of 0.4 as specified in the prior emission error covariance matrix, because the equivalent error of bottom‐up emissions is 0.057 after adjusting by γ , and we find that when the random noise in the emissions are too large, the inversion hardly improves the emission estimates. More details of the pseudo observation setup are provided in section S1 in the supporting information. The pseudo observations are then used to optimize the prior emissions using the 4D‐Var method. We evaluate the performance of the inversion by comparing the posterior emissions with the true emissions. In the σ = 0.5 case, the normalized mean square error (NMSE) of NO (SO2) has reduced by 63.6% (82.9%) and normalized mean bias (NMB) of NO (SO2) has reduced by 24.0% (40.7%) after the joint inversion; in the σ = 1.5 case, the NMSE of NO (SO2) has reduced by 73.0% (34.8%) and the NMB of NO (SO2) has reduced by 39.0% (70.4%).

The magnitude of observation noise also affects the joint 4D‐Var performance. In Table 1, the ratio of joint to single species posterior SO2 (NO) NMSE and NMB decrease when the observation noise of SO2 (NO2) increases, demonstrating the accuracy of the joint posterior increases as observation noise of the optimized species increase (although the NO NMSE in the joint inversion is persistently higher than that of the single species inversion, possibly owing to the higher uncertainties in SO2 observations compared to NO2 observations). This improved performance of joint inversion benefits from incorporation of the other species observations. More comparisons of joint 4D‐Var with other posteriors are in section 3.4.

Table 1

Ratio of NMSE and NMB of Posteriors From Joint to Single Species Inversion Constrained by Pseudo Observations With Different Standard Deviations

Variable	Value
NO2 pseudo observation standard deviation	0.2	0.2	0.3	0.5
SO2 pseudo observation standard deviation	0.3	0.5	0.2	0.2
SO2 NMSE of joint/SO2 NMSE of single	0.80	0.44
NOx NMSE of joint/NOx NMSE of single			1.54	1.21
SO2 NMB of joint/SO2 NMB of single	0.91	0.80
NOx NMB of joint/NOx NMB of single			1.00	0.99

Mass Balance Joint Inversion

Joint Mass Balance Framework

The mass balance approach was originally proposed to estimate single species top‐down emissions by scaling the prior emissions using the ratio of observed and simulated column density based on a zeroth order Taylor expansion (Martin, Jacob, Chance et al., 2003). A finite difference mass balance was then put forward by Lamsal et al. (2011) and Cooper et al. (2017) to estimate changes in NO emissions using a first order Taylor expansion. We extend this finite difference mass balance to multiple species in order to simultaneously estimate NO and SO2 emissions using NO2 and SO2 observations, as derived below.

For any grid cell (i,j), the NO2 slant column density (SCDN(i,j)) is affected by NO emission (E N(i,j)) and SO2 (E S(i,j)) emission, and similar for SO2 slant column (SCDS(i,j)). We define slant column density matrix as and prior slant column density where SCDaN(i,j) and SCDaS(i,j) are prior slant column densities for NO2 and SO2. Top‐down and prior emissions are expressed as where E (i,j) and E (i,j) are prior emissions for NO and SO2, respectively. Sensitivities of NO2 and SO2 SCDs to emissions are expressed as

The Jacobian matrix is defined as

A first order Taylor expansion of SCD around prior model simulated NO2 and SO2 slant columns is

An expression for the top‐down emissions is thus

We use OMI NO2 and SO2 observations for SCD(i,j) and then weight top‐down (E(i,j)) and prior (E (i,j)) emissions for each species using maximum likelihood estimation as derived in Qu et al. (2017). Following Qu et al. (2017), uncertainties of 0.2 and 0.3 for all grid cells are calculated by averaging NO2 and SO2 observation errors over a month and are used in the weighting. The partial derivatives B are calculated using the central difference finite difference method with 10% perturbation to NO and SO2 emissions. When calculating the mass balance emission scaling factors for our simulation, we ignore the buffer region (the first three grid columns along each horizontal side) within which species concentrations are overwritten by those from the 4° × 5° simulations.

Evaluation Using Pseudo Observation Test

The same prior emissions and pseudo observations in 3.1.2 are used to test the performance of the joint mass balance inversion. In Table 2, we compare the NMSE and NMB of the true NO and SO2 emissions with the prior, mass balance (MB) posterior, and joint mass balance posteriors for the σ = 0.5 case in January 2010. Posterior NO and SO2 emissions have reduced bias (by 20.5–40.0%) and error (by 52.8–66.7%) using both methods. Joint mass balance posteriors have smaller error (by 8.1% for NO and 14.4% for SO2) and bias (by 6.9% for NO and 8.5% for SO2) than single species mass balance. The comparison for the σ = 1.5 case is shown in Table S2, for which joint mass balance have similar performance (difference in NMSE and NMB less than 2%) with mass balance.

Table 2

Error and Bias in NO and SO2 Emissions When σ = 0.5

Emission types	NMSE			NMB
	NOx	SO2	Total	NOx	SO2	Total
Prior	3.42	4.28	5.48	−50.47%	−49.93%	−100.40%
MB	1.24	2.02	2.37	−33.82%	−39.70%	−73.52%
Joint MB	1.14	1.73	2.07	−31.47%	−36.33%	−67.80%

Note. NMSE = normalized mean square error; NMB= normalized mean bias; MB =mass balance.

The limited improvement of NMB comes from the constraints from prior estimates in the mass balance. Without this constraint, the total NMB would have reduced to −33.36% for MB and −21.36% for Joint MB. Cooper et al. (2017) found NME reduce by 35.6% and 47.1%, respectively, without prior constraints in January when performing pseudo NO inversion test (hourly pseudo observations are generated using a different fossil fuel inventory) at a coarser resolution (2° × 2.5° and 4° × 5°). For comparison, our NME for posterior NO emissions from the single species mass balance have decreased by 52.3% without prior constraints.

The Role of Chemistry in the Joint Inversion

To study the reasons for the improved performance in the joint inversion, we first look at changes in NO, SO2, O3, and OH concentrations when uniformly decreasing NO and SO2 emissions by 30% at all grid cells in the domain. In Figure 1, when only decreasing NO emissions, O3 concentrations increase in populated areas (NO saturation regime). OH concentrations increase in the North China Plain (NCP), where populations and thus NO emissions are especially high. HO2 is therefore converted to OH through NO in these regions. OH decreases in the rest of the domain, where NO emissions are lower and OH is possibly converted to HO2 by carbon monoxide (CO) and volatile organic compounds, leading to increases of SO2 concentrations. In urban areas, only small increases or even decreases of SO2 occur, due to increases of OH. Changes in O3 and OH concentrations after decreasing SO2 emissions are about five times smaller than those after decreasing NO emissions. Decrease of SO2 leads to increase of OH through gas phase oxidation in most regions of China other than Yunnan, Beijing, Shandong, Hebei, and Henan. Increase of OH leads to increase of HO2 and O3 concentration in Southeast China and India, but a decrease of O3 in the NCP.

Figure 1

Percent changes of O3, OH and SO2 (NO) surface concentration after decreasing NO (SO2) emissions by 30% in January 2010.

In 4D‐Var, these chemical interactions are accounted for within the gradients calculated by the adjoint model. In mass balance, β values approximate the first order sensitivity of NO2 and SO2 columns to emissions. Negative β values mostly appear in remote regions (e.g., Tibet Plateau), where increase of NO leads to increase of OH and decrease of SO2 concentration. 81% of β over China and India are positive, reflecting decrease of OH after increasing SO2 emissions, which reduces oxidation of NO2. Sensitivities of NO2 and SO2 columns to emissions are different from sensitivities of surface layer concentrations due to different concentrations at each layer and impact of vertical mixing and meteorology. Calkins et al. (2016) also showed that in winters over China, columnar SO2 change is a better indicator of SO2 emission change than surface SO2 because of meteorological impacts.

Method Intercomparisons

We compare the performance of 4D‐Var and mass balance inversions by evaluating the total (NO + SO2) NMSE of the posterior emissions when compared to the true emissions. In Figure 2, blue bars show the decreases of NMSE when assimilating both NO2 and SO2 observations and optimizing NO and SO2 emissions over China in the σ = 0.5 case. Joint 4D‐Var posteriors have larger decreases of NMSE (by 15.1% on left panel and 5.4% on right panel) than joint mass balance posteriors. Orange bars sum up errors of two single species inversion—one assimilates NO2 and optimizes NO emission and the other assimilates and optimizes SO2. When using the same inversion setup as in section 3.1.2 (left panel), the joint mass balance has a larger decrease of error by 5.64% than summing up decreases from the separate single species mass balance inversions. However, the joint 4D‐Var has 4.18% less error reductions than the sum of two single species 4D‐Var inversions. This degraded performance of the joint 4D‐Var compared to single species 4D‐Var is caused by the large noise in observations and prior emissions. In a separate test with SO2 observation error reduced to 20% and zero noise in the prior emissions (first two bars in the right panel), the joint 4D‐Var decreases NMSE by 12.3% more than summing up two single species 4D‐Var inversions. This suggests that including observations with large uncertainties would possibly degrade the overall performance of the joint inversion.

Figure 2

Decrease of normalized mean square error (NMSE) for the σ = 0.5 case over China for the following scenarios: joint 4D‐Var inversion (VarJ), sum of error reductions of two individual 4D‐Var inversions (VarI), joint mass balance inversion (MBJ), sum of error reductions of two individual mass balance inversion (MBI), assimilate both species but only optimize SO2 emissions using 4D‐Var (OptS), assimilate both species but only optimize NO emissions using 4D‐Var (OptN), assimilate only SO2 but optimize both species emissions using 4D‐Var (ObsS), assimilate only NO2 but optimize both species emissions using 4D‐Var (ObsN). Sum of errors are added in quadrature. The left panel uses random noise with standard deviation of 0.2 for NO2 and 0.3 for SO2 when generating pseudo observations, whereas the right panel uses 0.2 for both NO2 and SO2. The same prior emission random noise as in section 3.1.2 is also applied to inversions in the left panel but not in the right one.

The decrease of NMSE when using both observations and optimizing both species (first four bars) are larger than observing both species but only optimizing single species (5th and 6th bars), and larger than observing single species but optimizing both species emissions (7th and 8th bars). Since joint 4D‐Var and joint mass balance inversions (the two blue bars) have better performance than other scenarios when observation errors are small, we combine these two in our hybrid joint 4D‐Var/mass balance framework to derive long‐term posterior NO and SO2 emissions.

Hybrid Joint Inversion

The hybrid 4D‐Var/mass balance approach was first derived to facilitate long‐term NO inversions (Qu et al., 2017). In a base year, 4D‐Var is applied to rigorously trace prediction error back to sources and identify emissions corrections. This 4D‐Var posterior emission inventory is then used as the prior inventory for mass balance inversions in other years. This approach is adopted as the 4D‐Var approach is too computationally expensive to be applied in all years, compared to the mass balance inversion, which is about 20 times faster per year. We perform pseudo observation tests in January of each year from 2005 to 2012 to evaluate the performance of this hybrid inversion to recover emission trends. We first perform joint 4D‐Var inversions for January 2010. The posterior NO and SO2 emissions in the base year are then used to simulate NO2 and SO2 columns in January of other years (2005–2012). We then apply the joint mass balance derived in section 3.2 to recover the emission trend for both NO and SO2.

For these pseudo observation tests of the hybrid joint inversion, the true NO and SO2 emissions have a 10% annual growth rates compared to 2010. All anthropogenic prior emissions have a 50% low bias compared to 2010 levels and also include random noise (same standard deviations as in section 3.1.2). Therefore, the NMSE in the prior emissions is increasing from 2005 to 2012 (Figure 3). In the first three years, the NMSE in both the prior and posteriors are small and differences among posteriors are trivial. From 2008 and after, the hybrid posterior has smaller NMSE (by 41–64%) than applying mass balance alone; NMSE in the joint NO and SO2 posteriors are 3–12% and 0–18% smaller than posteriors from single species inversion. Although the percent decrease of NMSE in the posterior compared to that in prior are similar (80–86% for single and joint hybrid NO, 57–76% for single and joint mass balance NO, 74–81% for single and joint hybrid SO2, and 27–60% for single and joint mass balance SO2), the absolute NMSE values in the posteriors increase each year. This can be attributed to restriction of posterior emission to their prior values in this Bayesian analysis, the impact of meteorology leading to opposite trends in column concentrations from the true emissions, and differences in observation density and quality throughout the period (Qu et al., 2017). Similar comparisons for NMB over the 8 years are shown in Figure S2.

Figure 3

Normalized mean square error (NMSE) of SO2 (left) and NO (right) emissions when compared to their true emissions over China for January, 2005 to 2012 (σ = 0.5). The true emissions have a 10% annual growth rate of anthropogenic emissions compared to 2010. The 2010 anthropogenic emissions with 50% low bias are used as prior emissions for all 8 years. Nonanthropogenic emissions from bottom‐up inventories are used for corresponding years in both prior and true emissions.

We generate pseudo observations using the same chemical transport model as the one we use to perform the inversion, which implicitly neglects uncertainties in chemistry and transport and thus reflects “best case” scenarios. When using real OMI observations in the inversion, inaccuracies in modeled transport and chemistry can increase uncertainties in the posterior emissions. Therefore, reductions in NMSE and NMB in emissions using real satellite observations are expected to be smaller than those reported here.

Joint Inversion With OMI Observations

In this section, we use OMI NO2 and SO2 retrievals and the hybrid joint inversion method to derive posterior emissions for East Asia. Given the large uncertainties in OMI SO2 retrievals in countries other than China and India (Qu et al., 2019), we only focus on emissions in these two countries. We start with evaluating 4D‐Var performance in 2010 and then apply the hybrid method to study the trends of these two pollutants from 2005 to 2012.

4D‐Var Inversion in 2010

One challenge of the joint 4D‐Var inversion is to balance the contribution of observation terms for different species in the cost function. A hypothesis being tested in this work is that jointly assimilating these species has synergistic impacts and thus giving equal weight to their observations is by design. While a break from a strict Bayesian derivation of the cost function, we subjectively here aim to weight the contributions of NO2 and SO2 observations equally in our inversion, even though current measurement capabilities do not naturally provide a balanced view of these two species; otherwise, assimilation of SO2 does very little to impact the inversion when assimilating NO2 unless we were to perform many more iterations of the optimization algorithm than is computationally feasible, as the inversion will first only correct for the terms that make the largest contributions to the cost function. To accomplish this, we use the ratio of number of observations of NO2 to SO2 (α, shown in Table S1) to scale the observation terms in the cost function and check the reduction of NO2 and SO2 prediction error in Table 3. SO2 prediction errors have generally larger reductions than those of NO2, especially during spring, summer and autumn when SO2 observation errors are smaller than 1016 molec/cm2 (red solid line in Figure 4). The domain‐wide and nation‐wide summer peak of absolute NO2 observation errors in the NASA standard product are mainly caused by the larger absolute uncertainties in the summer in remote regions, although in major economic regions NO2 absolute uncertainties have winter peaks in the NASA standard product. The larger absolute NO2 uncertainties in the remote regions lead to smaller constraints on NO emissions in these locations in the summer. However, the seasonalities of absolute NO2 uncertainties in the product used in this study may be different from those of other products. For instance, Boersma et al. (2018) show a winter peak of OMI NO2 uncertainties over China. Monthly variations in the magnitude of the cost function reduction for each species are related to balance of observation and parameter errors. These balances all involve uncertainties in prior emissions and observations, which are not considered when we scale α and γ using the number of observations. However, the relative changes of these two terms are of the same magnitude, suggesting that the weighting strategy using α still helps balance the role of the two species in the optimization.

Table 3

Cost function Reductions in NO2 and SO2 Observation Terms in 4D‐Var Inversions in 2010, Reported as the Percent Difference of the Species‐Specific Prediction Error Terms in the First and Converged Iteration

	January	February	March	April	May	June
SO2 (%)	41.7	48.5	47.6	82.5	19.8	20.1
NO2 (%)	21.9	40.0	31.3	12.8	8.6	5.9

Figure 4

Monthly mean uncertainties in Ozone Monitoring Instrument NO2 (blue) and SO2 (red) observations over the entire East Asia domain (solid line), over China (dotted line), and over India (dashed line) in 2010.

The equivalent prior error after scaling by γ is 5.7%. This suggests that in regions where emissions most directly impact the cost function, a 40% uncertainty estimate may be too large. In January 2010, 99.4% (single species inversion) and 91.4% (joint inversion) of SO2 scaling factors and 94.9% (single species inversion) and 96.1% (joint inversion) of NO scaling factors change less than three times the standard deviation (5.7%), consistent with the relationship between standard deviation and data coverage in a normal distribution.

Changes in emissions when assimilating NO2 and SO2 individually and jointly are shown in Figure 5. Incorporation of NO2 observations leads to up to a 242% increase of SO2 emissions over NCP, Guizhou Province, and South Korea (first row). Posterior NO emissions from the joint inversion also increase in the NCP by up to 49%. These increases in the joint posterior SO2 and NO emissions happen in regions where prior SO2 simulations are lower than OMI observations, as shown in Figure 6. In these high‐NO regions, increases of NO emissions leads to decreases of OH and increases of SO2 concentration to match OMI SO2 observations.

Figure 5

Difference of posterior and prior emissions when assimilating NO2 and SO2 observations individually (left column) and jointly (center column) and differences in posteriors from joint and single species inversions (right column), all for January 2010.

Figure 6

Differences in GEOS‐Chem simulation (using prior emissions) and Ozone Monitoring Instrument SO2 slant column density in January 2010.

Joint posterior SO2 emissions also decrease by up to 99% more than single species posterior in India and Sichuan, Henan, and Hubei Province of China. In remote regions, increase of NO leads to increase of OH, which decreases SO2 column. Therefore, NO emissions increase in most rural areas in India, and Henan and Hubei provinces in China to match the lower OMI SO2 column than prior model simulation (Figure 6).

Posterior SO2 emissions barely change in countries other than China and India in single species inversion (Qu et al., 2019) and also in this joint inversion (less than 4%). Therefore, we focus only on these two countries in this work. Seasonalities of 4D‐Var posteriors and bottom‐up emissions in China and India are compared in Figure 7. In China, the joint SO2 posterior has larger seasonal variation than the single species posterior. The annual budgets of SO2 emissions in China decrease from 12.4 Tg S in the prior to 11.8 Tg S (single species) and 11.6 Tg S (joint) in the posteriors. In India, both single species and joint posterior SO2 emissions have similar seasonality, which is different from the prior. The annual budget of posterior SO2 emissions in India is 3.2 Tg S (single species) and 2.5 Tg S (joint), which also decreases compared to the prior (4.4 Tg S). The seasonality of NO emissions does not change significantly in China and India. The annual budgets of NO emissions are 1–3% larger in the joint posterior (China 7.9 Tg N; India 3.3 Tg N) than the single species posterior (China 7.8 Tg N; India 3.2 Tg N), and both are within 4% of the prior (China 8.1 Tg N; India 3.3 Tg N).

Figure 7

Seasonality of SO2 and NO emissions in China (left) and India (right) for 2010 from 4D‐Var. OMI = Ozone Monitoring Instrument.

The error weighted difference between predicted and simulated NO2 and SO2 slant columns (forcing term) is lower in the summer than in the winter, which therefore leads to smaller reductions of the cost function in the summer. The largest reduction of NO2 forcing in the posterior simulation occurs mostly in eastern China and northern India. Monthly mean reductions of NO2 forcing in the posterior simulation over the entire domain go up to 14.6% in 2010. The largest SO2 forcings are in India, which are mainly positive in the prior simulation. Over China, prior SO2 forcing is mostly negative over North China plain and positive in the rest of the regions. The largest reductions of SO2 forcing are over all of India and the eastern part of China. Monthly SO2 forcings have reduced by up to 29.6% in the posterior in the entire domain. The mismatch of the posterior simulation and observations are caused by prior constraints and observation errors.

The posterior NO and SO2 concentrations are generally larger over the North China Plain but smaller over the rest of the region (Figure 8). Changes in O3 and OH concentration in the joint posterior simulation (Figure 8) are generally in the same direction. Monthly changes of surface O3 concentrations are between −35.6% and 142.7%; changes of OH concentrations are between −25.7% and 100.8%. Their changes are larger in the fall and the winter than in the spring and summer. Changes in CO concentrations are within 1%. Changes of surface O3 and OH concentration in the joint inversion are mainly driven by changes in NO emissions. If only optimizing SO2 emissions using SO2 observations, monthly changes in surface O3 concentrations are less than 1.1 ppbv and lead to increases of OH concentration (within 21%) in India and south China.

Figure 8

Changes of O3 and OH concentrations before and after jointly optimizing NO and SO2 emissions in January, April, July, and October of 2010.

Trends of NO and SO2 Emissions

The hybrid joint inversion method is applied to recover the trend of NO and SO2 emissions in China and India. The spatial distribution of emission changes using the single‐species mass balance and joint mass balance approaches are shown in Figure S3. The higher joint posterior SO2 emissions in the NCP and lower values in Henan and urban areas in India are consistent with the differences from 4D‐Var results (Figure 5). The interannual variations of joint and single species posteriors are more similar than their spatial distribution, as shown in Figure 9.

Figure 9

Annual budget of NO and SO2 emissions in China (left column) and India (right column). The magenta lines are the sum of anthropogenic emissions from EDGAR v4.3.2 and natural emissions used in GEOS‐Chem; the green line is the sum of anthropogenic emissions from MEIC (Wang et al., 2018) and natural emissions. We also compare our emissions with top‐down SO2 emissions from Li, McLinden, et al. (2017) shown in pink and four top‐down NO emissions from Ding, Miyazaki, et al. (2017), that is, MIROC model based on ensemble Kalman filter and Ozone Monitoring Instrument data (yellow), CHASER model based on ensemble Kalman filter and OMI data (cyan), DECSO model based on extended Kalman filter and OMI data (orange), and DECSO model based on extended Kalman filter and GOME data (purple).

In Figure 9, we also evaluate the trend of the hybrid posterior emissions by comparing them with bottom‐up emissions trends from Multi‐resolution Emission Inventory for China (MEIC) and Emissions Database for Global Atmospheric Research (EDGAR) (v4.3.2). Posterior SO2 emissions in China start to decrease in 2008 and drop by 16.5% (joint) and 16.2% (single species) till 2009, an amount of decrease which is consistent with the 11.8% decrease in the MEIC inventory, although MEIC SO2 starts to decrease in 2007 instead of 2008. This decrease in SO2 emissions is not reflected in the EDGAR inventory (magenta line). On the other hand, posterior SO2 emissions show qualitatively consistent increases with the EDGAR inventory in India. However, posterior growth rates (18.7% for single and 19.2% for joint) are about a quarter of that in EDGAR (76.2%).

The turning points in the trends of posterior SO2 emissions in 2007, 2009, and 2011 in China are consistent with top‐down estimates from Li, McLinden et al. (2017; bright yellow line in Figure 9) and Koukouli et al. (2018; estimated over a different domain, not shown in Figure 9), although differences in the magnitude of our posterior SO2 emissions are up to 29.6% compared with Li, McLinden, et al. (2017) and 33.7% compared with Koukouli et al. (2018; emissions in this study are summed over 15–50°N and 102–132°E for this comparison). Our top‐down estimate of 884 Gg S (joint posterior) and 887 Gg S (single species posterior) in China in April 2008 is 4.8% and 5.2% higher, respectively, than the 843 Gg S estimates in Wang et al. (2016). The growth rate of top‐down SO2 emissions in India in this study is 77.4% lower than that in Li, McLinden, et al. (2017). However, the analysis in Li, McLinden, et al. (2017) is based on trends in 47 large point sources of SO2 emissions, whereas our trends are based on national total emissions. For comparison, if we consider only grid cells that have emissions larger than 106 kg S per grid in the prior inventory, our posterior SO2 growth rate would be 29.1% larger than the growth rate using all posterior emissions in India, which would more closely agree with Li, McLinden, et al. (2017).

NO emissions in China increase from 2005 to 2011 by 24.9% (joint inversion) and 25.0% (single species inversion) and decrease from 2011 to 2012 by 1.8% (joint) and 1.6% (single species). In comparison, the MEIC and EDGAR bottom‐up inventories have larger growth (43.2% and 30.0%) from 2005 to 2011. Joint NO emissions' growth (15.9%) in India is 56.2% larger than that of single species posterior (10.2%) from 2005 to 2012, with nearly doubled growth rate from 2009 to 2012 (6.5% in joint one, 3.3% in single species one), but both are smaller than the 59.9% growth in EDGAR inventory. The magnitude of posterior NO emissions from the single species inversion is 5.8–13.2% larger than the result from the same system in Qu et al. (2017) due to use here of a different version of the OMI NO2 product (v2.1) and different cloud fraction filtering criteria (correction factor, CF < 0.2 in this work but only CF = 0 is used in Qu et al., 2017). Although the magnitude of posterior NO emissions have good consistency with the MEIC inventory (including natural sources, from GEOS‐Chem), especially after 2008, posterior NO emissions are likely underestimated due to the usage of spatially coarse prior information in the OMI retrievals, which lead to systematic low bias in NO2 retrievals (Goldberg et al., 2018; Laughner et al., 2016). The omission of halogen chemistry in the version of GEOS‐Chem used here possibly causes higher simulated troposphere NO2 burden than simulation with halogen reactions (e.g., about 4.2% in 2005 in Sherwen et al., 2016) and also contributes to underestimates in posterior NO emissions.

However, the annual budgets of top‐down NO emissions from joint and single species inversions are higher than four other top‐down estimates in Ding, Miyazaki, et al. (2017) by 12.1–82.9% (joint posterior) and 9.2–76.8% (single species posterior). These emissions are derived using MIROC‐Chem, CHASER, and DECSO chemical transport model based on ensemble Kalman filter and extended Kalman filter algorithms. Posterior NO emissions in this study have a peak in 2011 consistent with those from MIROC and CHASER, and also consistent growth rate within the range of 22.7% to 28.7% from 2005 to 2011.

To compare the trend of hybrid joint posterior emissions and the trend of 4D‐Var joint posterior emissions, we performed 4D‐Var inversion for each May from 2005 to 2012, since the trend in this month is most representative of the interannual variation of NO and SO2 emissions in the hybrid posterior emissions. While the sign of the changes in posterior emissions are mostly consistent between the two different inversion methods, the magnitudes are slightly different. As shown in Figure S4, the changes in NO emissions are smaller in the 4D‐Var posteriors than in the hybrid posteriors; the opposite is generally found for SO2. Consequently, the correlation coefficients between 4D‐Var and hybrid posterior emissions are relatively higher (0.95 for NO in China, 0.95 for SO2 in China, 0.69 for NO in India, and 0.97 for SO2 in India), but the absolute differences between emission changes compared to 2005 level are slightly larger, that is, within 10.6% for NO in China, 10.9% for SO2 in China, 25.0% for NO in India, and 20.7% for SO2 in India.

Evaluation With Surface Measurements

Surface Concentrations in China

We next evaluate the hybrid posterior emissions by comparing monthly mean GEOS‐Chem simulated concentrations with monthly mean surface measurements. The locations of measurements sites and comparisons to GEOS‐Chem prior and posterior simulations are shown in Figures S5 and S6. Posterior SO2 simulations have improved correlation with surface measurements in every month of 2010, and the joint posterior emissions have the highest correlation in 9 months of the year (top right panel of Figure 10). The NMB and NMSE of the SO2 joint posterior reduced by up to 26% and 18% compared to single species posterior in most months between November and February. These improvements are related to the relative larger uncertainties in OMI SO2 retrievals in these four months (Figure 4), when incorporation of NO2 observations helps improve SO2 emission estimates. This behavior is also corroborated by the pseudo observation test in section 3.1.2, where the larger the noise in the SO2 observation, the better the performance of joint inversion than SO2 single species inversion to improve SO2 emissions. On an annual basis, the NMSE of surface SO2 concentration reduce from 0.77 (the prior) to 0.58 and 0.61 in the singe and joint posteriors, respectively. Spatial distributions of annual mean SO2 concentrations in GEOS‐Chem and at monitoring sites are compared in section S2. Joint posterior SO2 surface concentrations have smaller NMSEs over cities than simulations using the prior and single species posterior SO2 emissions. Over all provincial capitals in China, the NMSE of the annual mean surface SO2 concentration reduces from 0.66 to 0.40 (single posterior) and 0.36 (joint posterior); over provincial capitals in East China, the NMSE reduces from 0.30 to 0.19 (single posterior) and 0.15 (joint posterior).

Figure 10

Monthly mean surface SO2 concentration and statistics of monthly mean SO2 sampled at locations of monitoring sites in GEOS‐Chem simulations and surface measurements in 2010 from 608 sites, which have measurement throughout the year. The magenta vertical error bars in the top left panel show the standard deviations of the surface measurements. NMSE = normalized mean square error.

The NO2 measurement techniques employed in China's monitoring sites detect NO concentration after converting NO2 and some unknown fraction of NOz to NO. We compare monthly GEOS‐Chem surface NO and NO2 concentration with reported NO2 from the monitoring sites in Figure S7, where simulated NO2 are mostly biased low and simulated NO mostly biased high. Following Lamsal et al. (2008), we use a correction factor to account for the interference in the NO2 measurements by the molybdenum converter. Given that we have three sets of model simulations using prior, hybrid single posterior and hybrid joint posterior emissions, we calculated the correction factor for each GEOS‐Chem simulation and divided the simulated NO2 concentrations by the corresponding correction factor to convert simulated NO2 to the measured species. The monthly mean correction factors are between 0.79 and 0.83 in 2010, and the annual mean correction factors are between 0.83 and 0.85 from 2005 to 2012 in China.

Another factor contributing to the low bias in the simulations is the model resolution. As pointed out in Valin et al. (2011), model simulations with resolution greater than 12 km may not accurately represent nonlinear effects of NO2 loss rates. Coarse resolution representation of large source of NO2 will be biased low due to enhanced OH concentration. We therefore focus more on the seasonality and interannual variation instead of the absolute magnitude of surface NO2 concentration when compared to measurements.

In Figure 11, surface measurements are higher than model simulations in all months. This can be caused by the NOz contributions to the measurements and biases inherent in comparing model simulations at the 0.5° × 0.667° resolution with point measurements (mostly in cities) for a short‐lived species. Still, the joint posterior has reduced NMSE from May to October and reduced absolute value of NMB from May to September compared to the single species posterior, corresponding to the months when OMI NO2 observations have relatively larger error (Figure 4) and benefit from synergistic change of OH and O3 concentration in the model by assimilating OMI SO2 observations. These conclusions from NO and SO2 emissions constrained by OMI observations are consistent with results from pseudo observation tests in sections 3.1.2 and 3.4.

Figure 11

Comparison of monthly mean surface NO2 concentration in GEOS‐Chem simulations sampled at locations of monitoring sites with surface NO2 measurements from 614 sites in China, 2010. The magenta vertical error bars in the top left panel show the standard deviations of the surface measurements. NMSE = normalized mean square error.

Trends from 2005 to 2012 of surface SO2 and NO2 concentrations in China are similar in the joint and single species posterior simulation (Figure 12). The average of surface NO2 and SO2 concentrations from 2005 to 2012 are lower in the simulations than the measurements, i.e., measured surface SO2 (NO2) is 44.6 μg/m3 (35.7 μg/m3), 41.0 μg/m3 (22.1 μg/m3) in the prior, 34.6 μg/m3 (21.9 μg/m3) in the joint posterior, and 34.9 μg/m3 (22.4 μg/m3) in the single species posterior. Surface SO2 concentrations decrease by 28.3% (in situ), 10.9% (joint) and 14.9% (single) over the period, consistent with the trend in SO2 emissions. The low biases of simulated NO2 concentrations when compared with measurements are mainly from West China, where posterior simulations are 53.0–63.5% lower than measurements. In comparison, simulated NO2 concentration in East China are 3.5–17.2% lower than measurements. Posterior surface NO2 concentration increase by 19.3% (joint) and 19.8% (single) from 2005 to 2011 and decrease by 2.7% (joint) and 2.6% (single) from 2011 to 2012. The surface measurements only increase by 5.5% from 2005 to 2011, reflecting potential issues with model deposition and vertical mixing or consistent measurements of NOz species across monitoring sites. Still, all simulated surface NO2 and SO2 concentrations are within one standard deviation of surface measurements. Percent changes of surface SO2 and NO2 concentration compared to 2010 are shown in Figure S8.

Figure 12

Annual mean surface SO2 and NO2 concentration in China. There are 1,118 monitoring sites providing measurements between 2005 and 2012, but only 272 (283) of them have SO2 (NO2) data in every month from 2005 to 2012. We only use sites that have data cover the whole period in the in situ trend calculation. Model simulations are sampled at grid cell that has in situ measurements. The magenta vertical error bars are the standard deviations of the surface measurements.

To evaluate the trend of 4D‐Var and hybrid posterior simulations, we compare the simulated surface concentrations using posterior emissions in May of each year estimated from these two methods (described in the last paragraph of section 4.2) with measurements in China. Compared with the changes in the measurements from 2005 to 2012 (Figure S9, +11.0% for NO2 and −23.1% for SO2), the changes in 4D‐Var posterior simulations (+22.2% for NO2 and −16.7% for SO2) are more consistent than the changes in hybrid posterior simulations (+25.1% for NO2 and −2.7% for SO2). The 4D‐Var posterior simulations also have better correlations (R = 0.79 for NO2 and R = 0.83 for SO2) with the surface measurements than the hybrid posterior simulations (R = 0.78 for NO2 and R = 0.75 for SO2). However, the hybrid method facilitates long‐term emission estimates to be calculated much more efficiently than the 4D‐Var approach, with the sacrifice of accuracy.

Surface Concentration in India

The number of NO2 and SO2 monitoring sites from CPCB in India is much less than in China. Most sites only have data for a few years or months during the 2005–2012 period. Since the long‐term trend of these data may only reflect spatial variability in SO2 concentrations across different locations where measurements are available in different years instead of interannual variation, we only use SO2 concentrations in 2010 for the evaluation. Prior simulations of surface SO2 concentrations are higher than the measurements in all months of 2010; posterior simulations have reduced bias and NMSE than the prior except for March to July (Figure 13). Joint posterior surface SO2 concentrations are lower than single species posteriors by 10–65% except for July. The large NMSE for the joint posterior results in June is caused by few measurements in that month and the fact that the prior simulation of surface SO2 concentrations is lowest in this month.

Figure 13

Comparison of GEOS‐Chem surface SO2 concentration with surface measurement in India in 2010 (from 14 sites). Posterior emissions are from the hybrid inversion. The magenta vertical error bars in the top left panel are the standard deviations of the surface measurements. The numbers inside the first panel stand for total number of daily measurements from all sites in each month. NMSE = normalized mean square error.

The annual averages of simulated surface NO2 concentrations (Figure 14) are 65–67% (prior: 12.0 ppb, single species posterior: 11.7 ppb, joint posterior: 11.5 ppb) smaller than the surface measurements (34.4 ppb). The lower values of simulated NO2 concentrations are related to comparison of concentrations from the 0.5° × 0.667° resolution model to only a single measurement site per grid cell, which may be located close to localized sources. In China, concentrations above 15 ppb are more widespread and thus the in situ measurements appear to better represented by the model at this resolution.

Figure 14

Comparison of GEOS‐Chem simulations with surface NO2 measurements in India in 2010 (from 14 sites). Posterior emissions are from the hybrid inversion. The magenta vertical error bars show the standard deviations of the surface measurements. The numbers inside the figure stand for total number of daily measurements from all 14 sites in each month.

Discussion and Conclusions

We combine the GEOS‐Chem adjoint 4D‐Var joint inversion framework with a newly developed joint mass balance method to derive long‐term NO and SO2 emissions in East Asia using simultaneous observations of NO2 and SO2 from OMI. In pseudo observation tests, this hybrid joint posterior has a larger decrease of NMSE (5.4–11.5%) than using a joint mass balance alone. Joint inversion has the potential to improve NO and SO2 emission estimates through synergistic change of OH and O3 concentration; however, the inclusion of both NO2 and SO2 observations does not always improve model performance when noises in observations and prior emission are large. This impact from observation noise is further evident when comparing model simulations using posterior emissions with surface measurement in China and India—the joint inversion posterior has reduced NMSE for SO2 (NO) in winter (summer) months when uncertainties in SO2 (NO2) observations are large and benefit from the inclusion of NO2 (SO2) observations.

Percentage changes of joint posterior emissions from the single species posterior emissions from 4D‐Var and mass balance inversions range up to 242% in January in individual grid cells. Still, annual emissions budget, seasonality, and interannual variations of joint and single species posterior emissions are similar. In 2010, posterior SO2 emissions in China (11.8 Tg S from single species and 11.6 Tg S from joint) and India (3.2 Tg S from single species and 2.5 Tg S from joint) are both lower than the prior (12.4 Tg S in China and 4.4 Tg S in India); posterior NO emissions also reduce or stay the same in China (single species 7.8 Tg N, joint 7.9 Tg N) and India (single species 3.2 Tg N, joint 3.3 Tg N) than the prior (China 8.1 Tg N, India 3.3 Tg N). Joint posterior SO2 emissions are lower and NO are higher than single species posterior emissions, since NO emissions increase in NO limited regions to help increase SO2 concentrations. The lower values in the posterior than prior could be related to model resolution errors in representing NO2 and SO2 columns at the satellite footprint scale (Valin et al., 2011) and the underestimates of NO2 peak column concentration due to coarse resolution in the retrieval prior information (Goldberg et al., 2018). In China, NO emissions peak in summer due to increase of natural sources (Qu et al., 2017), whereas SO2 emissions peak in winter. Posterior SO2 emissions peak in the summer in India, whereas NO emissions peak in spring. Both joint and single species SO2 posterior emissions decrease by 16.2% (single species) and 16.5% (joint) in China since 2008; NO emissions increase by 25.0% (single species) and 24.9% (joint) up through 2011. NO and SO2 emissions in India increase by 10.2–19.2% throughout the period. Hybrid posterior uncertainties are not estimated because they depend upon uncertainties in the 4D‐Var posterior emissions in 2010 that require additional computations (beyond the scope of this work) to quantify (Bousserez et al., 2015). Instead, we consider comparisons with other independent emission inventories and in situ measurements to provide evaluation of the accuracy of our posterior emissions.

The trend of surface SO2 concentrations in the posterior simulations does not precisely track those of surface measurements, suggesting possible errors in assimilated meteorology, modeled deposition, or inconsistent measurement methods, which warrant further investigation. The seasonality of surface NO concentrations from joint posterior simulations has better correlation with surface measurements than those from the prior or single species posterior simulations.

This study demonstrates the benefits of synergistic change of OH and O3 concentration when performing emission inversions for rapidly oxidized species. However, only chemical interactions of NO and SO2 are incorporated in inversions. Correlations in emission errors of NO and SO2 will be addressed in future work. The performance of the joint inversion depends on the relative uncertainties of SO2 and NO2 observations, which in theory would allow the data assimilation system to de‐emphasize the importance of highly uncertain observations. However, precise quantifications of measurement uncertainties are not known. While this paper focuses on improvements to just simulations of SO2 and NO2 concentrations, the impacts of the joint SO2 and NO2 inversion on other species such as aerosols and O3 suggest that future studies may benefit from expanded multispecies, multicomponent data assimilation. Works such as Miyazaki et al. (2017) have included, beyond constraints on NO emissions, direct assimilation of remotely sensed O3. HCHO assimilation would also potentially change O3 concentrations (e.g., Hamer et al., 2015; Kaiser et al., 2018) and thus help enforce consistent NO and SO2 emissions; joint assimilation of AOD could also potentially improve NO and SO2 emission estimates, as suggested by previous studies wherein AOD assimilation leads to improved agreement of model estimated trace gas concentrations with retrievals from OMI (Xu et al., 2013). Retrievals from TROPOspheric Monitoring Instrument and future geostationary satellites (e.g., Geo‐stationary Environmental Monitoring Spectrometer, and Tropospheric Emissions: Monitoring of Pollution, Sentinel 4) could provide a higher density of NO2 and SO2 observations to improve future joint data assimilation.

Supporting Information S1

Click here for additional data file.