Satellite Data Reveals a Common Combustion Emission Pathway for Major Cities in China

Extensive fossil fuel combustion in rapidly-developing cities severely affects air quality and public health. We report observational evidence of decadal changes in the efficiency and cleanness of bulk combustion over large cities in mainland China. We combine air quality retrievals from mature satellite instruments across 2005-2014 to estimate the trends in enhancement ratios of CO and SO2 to NO2 (ΔCO/ΔNO2 and ΔSO2/ΔNO2) over these cities and infer emergent bulk combustion properties. Our 15 results show a robust coherent progression of declining-to-growing ΔCO/ΔNO2 (-5.4±0.7% to +8.3±3.1%) and slowly-declining ΔSO2/ΔNO2 (-6.0±1.0% to -3.4±1.0%) from Shenyang, Beijing, Shanghai, to Shenzhen relative to 2005, which is not evident in the trends of emission ratios reported in Representative Concentration Pathway (RCP8.5) inventory. This progression is likely due to a shift towards cleaner combustion from industrial and residential sectors in Shanghai and Shenzhen, which is presently 20 obfuscated by China’s still relatively higher dependence on coal. Such progression is well-correlated with economic development and traces a common emission pathway that resembles evolution of air pollution in more developed cities. Our results highlight the utility of augmenting observing and modeling capabilities by exploiting enhancement ratios in constraining the time variation of emission ratios in current inventories. The ability to monitor combustion efficiency and effectiveness of pollution control 25 becomes increasingly important in assessing sustainable control strategies as cities and/or countries continue to socioeconomically develop.


Introduction
Urban agglomeration, particularly megacities (i.e., cities with >10 million inhabitants) are expected to continue growing (in size and number) over the coming decades (Jalkanen, 2012;World Bank, 30 2015). This is especially problematic since it is in these large cities where anthropogenic activities are most intense, accompanied by immense energy consumption mainly in the form of fossil fuel combustion (Mage et al., 1996;Kennedy et al., 2015). These lead to enhanced emissions of air pollutants, greenhouse gases, and waste energy, largely impacting air quality (AQ), climate, and ecosystems (Baklanov et al., 2016, Lelieveld et al., 2015. At present, estimates of city-to-national-scale emissions from fossil fuel 35 combustion remain uncertain, especially in rapidly-developing regions where combustion is still poorly characterized due to the lack of detailed information on energy use, combustion practices, and pollution control strategies (Streets et al., 2013;Creutzig et al., 2015). This is also confounded by larger uncertainties on other sources of pollution that may be associated with urbanization (e.g., deforestation, agriculture, and fires). These alone preclude us to accurately assess the changes in atmospheric 5 composition due to anthropogenic activities at scales that are relevant to AQ, energy, and environmental policy (National Academies of Sciences, Engineering, and Medicine, 2016). Such is the case for cities in China even with the scientific attention the country has received in the past decades. As China grew into the world's second largest economy, its rapid development resulted to substantial emissions (Richter et al., 2005) and more frequent occurrences of most severe pollution 10 events in many of its megacities, most notably Beijing (Guo et al., 2014). These affect not only local AQ and public health but are reported to impact hemispheric-to-global atmospheric environment (Lin et al., 2014;Verstraeten et al., 2015). Along with the growth of these cities is a growing body of evidence of decreasing emissions and associated pollution levels in some cities in China, pointing to important changes in AQ as a result of development, AQ management, and regional-to-national socioeconomic 15 initiatives embodied within its Five-Year Plans (FYP) (Reuter et al., 2014;Krotkov et al., 2016;van der A et al., 2017;Sun et al., 2018;Koukouli et al., 2018). However, these changes in AQ as a result of efforts to control air pollution are still obfuscated at present by the increase in combustion activities along with uncertainties in bottom-up emission inventories and diversity in economic structure and growth across cities (Wang and Hao, 2012;Mi et al., 2017). Monitoring these reductions at city scale remains to be a 20 challenge especially when narrowly viewed within the context of a single pollutant, and more so when attributing them to a particular emission sector. Fossil fuel emissions from an evolving megacity follow a pattern that can be potentially monitored and refined by combining observational constraints on combustion activity (abundance of combustion products) with efficiency and effectiveness of pollution control strategies or 'cleanness' (enhancement 25 ratios of these products) (Silva et al., 2013;Hassler et al., 2016;Silva and Arellano, 2017;Tang et al., 2018), alongside information on the state of socio-economic development (e.g., gross domestic product (GDP) or income) and a priori estimates from bottom-up emission inventories. In particular, the 'cleanness' of combustion of a known fossil fuel type can be determined stoichiometrically by measuring the relative abundance of intermediate products such as carbon monoxide (CO), nitrogen oxides (NO X ), 30 sulfur dioxides (SO 2 ), and soot particles with final products like carbon dioxide (CO 2 ). Please see Methods section for more details. Most of these products are currently monitored as criteria pollutants by surface measurement networks and as tracers of pollution by satellite remote sensing (Streets et al., 2013;Duncan et al., 2014). In fact, these combustion products are revealed in space as very distinct bulk enhancements over a megacity metropolitan location in marked spatial contrast with the city's surroundings (Bechle et 35 al., 2011;Lamsal et al., 2013). At a scale of a megacity being monitored from space, these enhancements are analogous to smoke plumes coming from a stationary smokestack. And so, observations of these megacity plumes enable us to monitor bulk anthropogenic activity and transboundary pollution. They have also been used in recent years to refine the spatiotemporal distribution of emissions (Lamsal et al., 2013;Hakkarainen et al., 2016), to indicate bulk combustion efficiency, inter-megacity differences and 40 fire phase (Silva et al., 2013;Silva and Arellano, 2017;Tang and Arellano, 2017), and to infer fossil fuel CO 2 emissions (Konovalov et al., 2016) among others. From an annual to decadal standpoint, it is reasonable to interpret the long-term changes in spatial covariations between these observed pollutant enhancements within the megacity to reflect dominant shifts in bulk combustion characteristics (e.g., changes in fuel mixture and technology practice), which can then be indicative of an emission pathway for a given megacity (e.g., Parrish et al., 2002;Parrish, 2006;Russell et al., 2012;Silva et al., 2013;Hassler et al., 2016;Silva and Arellano, 2017). Data sampling and collocation issues, as well as retrieval 5 information content and chemical nonlinearities between these pollutants, do not quite manifest at decadal scales more than emission changes, especially when treated as a smokestack in the analysis. In this study, our goal is to uncover space-based evidence of dominant shifts in the cleanness of bulk combustion of a large cities across the current decade (through these ratios), associate these shifts to particular sectors, and identify a common emission pathway across these cities. Along the same line to 10 studies on environmental Kuznets curves (EKC, Stern, 2004) and human development (Lamb et al., 2014), we attempt to connect this pathway to economic growth by finding a power law relationship between the ratios observed for each major city in China and the city's GDP per capita. As cities in China grow, emissions from fossil fuel combustion evolve accordingly depending on the rate and type of socioeconomic development, technological innovation, and environmental policies (Chan and Yao, 2008;15 Bechle et al., 2011;Zhang et al., 2012;Wang et al., 2012;He and Wang, 2012;Luo et al., 2014;Koukouli et al., 2018;Sun et al., 2018). This evolution however cannot be reflected at shorter time scales. As a basis for comparison, pollution controls adopted in developed countries like United States and Europe, which followed a progression from first controlling SO 2 , CO, and then NO X (Crippa et al., 2016), reflect some aspects of decadal-scale sustainable development that can be brought to light in the case of China. 20 We analyze the emergent patterns of the 'cleanness' of bulk combustion in the past decade (2005-2014) based on enhancement ratios between intermediate products of combustion (∆CO/∆NO & and ∆SO & /∆NO & ) observed within each megacity and urban agglomeration in China. We use gridded monthly-averaged satellite retrievals of total columns of CO from Measurement of Pollution In The Troposphere (MOPITT), tropospheric columns of NO 2 from Ozone Monitoring Instrument (OMI), and 25 planetary boundary layer (PBL) columns of SO 2 from OMI to derive monthly estimates of these ratios using spatial regression analysis and subsequently derive estimates of the decadal trends of these ratios using time series analysis. We then compare these trend estimates to inferred trends from a couple of model-derived abundance ratios and several emission ratios from current bottom-up emission inventories including estimates based on the Representative Concentration Pathways scenario (RCP8.5) (Riahi et al., 30 2011). We also conducted a simple inverse analysis to update the contribution of major emission sectors in RCP8.5 to fit our estimates of decadal changes in enhancement ratios. Section 2 describes data and methods used in this study. Results and discussions are presented in Section 3. Section 4 is summary and implication of this study.

Study Region
We considered all 31 provincial capitals and five special cities (Beijing, Shanghai, Shenzhen, Tianjin, and Chongqing) in mainland China for our analysis. These cities comprise the main urban Atmos. Chem. Phys. Discuss., https://doi.org/10.5194/acp-2018-1121 Manuscript under review for journal Atmos. Chem. Phys. Discussion started: 4 December 2018 c Author(s) 2018. CC BY 4.0 License. agglomerations in the country (see Figure 1 for coverage). For purposes of finding long-term emergent patterns on its emission characteristics, we focused our analysis to 12 representative urban agglomerations. These 12 cities cover the four economic regions of China (i.e., East Coast: Beijing, Tianjin, Shanghai, Guangzhou, Shenzhen; Central China: Wuhan, Northeast China: Harbin, Shenyang; and Western China: Chengdu, Chongqing, Xian, Hohhot). Based on prior information from RCP8.5 and National Bureau of 5 Statistics of China (http://data.stats.gov.cn), these cities already exhibit largely diverse pollution and economic development attributes illustrated in Figure 1 as differences in magnitude, sectoral, and temporal distribution of emissions and GDP per capita for 2005 to 2014 between these cities. Our goal is to assess whether the long-term patterns that are seen in these a priori emission estimates are consistent with observations. We also considered Los Angeles and other large cities in the United States (New York  10 City, Chicago, Houston, Phoenix, Boston, Seattle, and Miami) for comparison.

Data
The main datasets used in this study are summarized in Table 1. This includes multiple satellite retrievals, representative emission inventories, and a couple of model simulations and chemical reanalysis.

Satellite Retrievals and Data Processing 15
We use the NASA Terra Measurement of Pollution In The Tropophere (MOPITT) version 6, Level 2, multispectral (Thermal Infrared/Near Infrared) retrievals of carbon monoxide (CO) total columns for CO (Deeter et al., 2014), tropospheric column retrievals from NASA Aura/ Dutch Ozone Monitoring Instrument NO & (DOMINO) v2.0 for NO & (Boersma et al., 2011), and Ozone Monitoring Instrument (OMI) Planetary Boundary Layer (PBL) SO & , version 3, Level 2 (Krotkov et al., 2006). We collected 20 daily MOPITT CO, OMI NO 2 , and OMI SO 2 retrievals that are available within a 2˚ by 2˚ area around each city center. This radius was selected to cover the extent of each city based on NO 2 footprints (Bechle et al., 2011;Lamsal et al., 2013) and geopolitical maps of city boundaries. We gridded each set of retrievals into 0.1-degree cells commensurate to the finest retrieval resolution among MOPITT and OMI. We then averaged them across each month to minimize spatiotemporal collocation issues. The gridded 25 data for each city domain is average across the month to minimize spatiotemporal collocation issues (see Table 1 for differences in sampling of MOPITT and OMI). As a result, there are 400 points for each species (CO, SO 2 , NO 2 ) per city and month. We note that these retrievals have been used in the past to study decadal changes for individual (or a pair of) pollutants but not to derive enhancement ratios (Krotkov et al., 2016). While CO retrieved from thermal infrared (TIR) radiances are mostly sensitive to 30 free tropospheric CO, it has also been reported to be capable of observing lower tropospheric CO, especially when retrieved from TIR and near infrared (NIR) radiances (Deeter et al., 2014). We recognize however that retrievals of SO 2 from OMI have been reported to exhibit low sensitivity to weak SO 2 signals, in particular to less than 30 to 70 kTon per year of point source emissions (Krotkov et al., 2016). While our spatial and temporal smoothing, along with anchoring our SO 2 analysis with NO 2 data (please see 35 later description of our regression analysis), should help in enhancing the SO 2 signal from cities with low SO 2 emissions, these SO 2 retrievals are useful as large SO 2 abundances are still observed across the majority of cities in China (Krotkov et al., 2016). We also used CO retrievals from the Infrared Atmospheric Sounding Interferometer (IASI), Level 2 (De Wachter et al., 2012), and tropospheric column NO 2 from FP7 QA4ECV OMI, v1 (Boersma et al., 2017) to verify consistency in our trend estimates.

Emission Inventories and Model Simulations
Multiple bottom-up emission inventories for CO, NO 2 and SO 2 are analyzed, namely Emission 5 Database for Global Atmospheric Research (EDGAR, Crippa et al., 2016), Representative Concentration Pathways (RCP8. 5, Riahi et al., 2011), Regional Emission inventory in ASia (REAS) version 2.1 (Kurokawa et al., 2013), and Hemispheric Transport of Air Pollution (HTAP, Janssens-Maenhout et al., 2015). We also use top-down emission estimates of CO and NO 2 from the Tropospheric Chemical Reanalysis (TCR) based on CHASER-LETKF assimilation system (Miyazaki et al., 2017). Since these 10 emission inventories have different spatial resolutions (see details in Table 1) and are available in the form of fluxes (units in kg/m 2 /s), we can upscale/downscale them by simply regridding into 0.1˚ by 0.1˚ cells similar to our approach for satellite data to facilitate comparison. We then consider all cells within the 2˚ by 2˚ area around the city center. For annual emissions, we only take the sum of all cells within the geopolitical boundary of the city (see Figure 2). All of the cities extend to less than the 2˚ by 2˚ area that 15 we set as our city domain. We also use model data for CO and NO 2 from the Community Atmosphere Model with Chemistry (CAM-chem; Gaubert et al., 2016) and TCR to derive CO and NO 2 abundance ratios associated with the bottom-up emissions used in these models (i.e., RCP in CAM-Chem and EDGAR in CHASER). The associated retrieval averaging kernels and prior information are applied to the daily-averaged model CO 20 and NO 2 vertical profiles of mixing ratios from CAM-chem and CHASER, along with appropriate spatial interpolation and/or partial column integrations. Since the spatial resolution (about 2˚~3˚) of CAM-chem and CHASER outputs that we analyzed are far coarser than 0.1˚, we only considered the associated abundance ratio rather than deriving enhancement ratio across the month where non-stationarity and nonlinearity issues are more likely to exist. 25

Deriving Enhancement Ratios using Spatial Regression Analysis
For each city, we regress the gridded monthly-average CO and SO 2 to NO 2 to calculate monthly enhancement ratios (∆CO/∆NO & and ∆SO & /∆NO & ). We use NO 2 as our control variable as NO 2 has the shortest lifetime (hours) among these products. Except for lightning, NO X is mostly produced from hightemperature anthropogenic combustion processes and because of its short lifetime it is observed as 30 distinctly and spatiotemporally local surface enhancements with relatively very low background concentrations. Along with the availability of NO 2 retrievals from satellites at fine spatial scale and over long period, NO 2 allows us to effectively identify intra-megacity combustion activities and define the urban extent (Bechle et al., 2011;Lamsal et al., 2013;Hakkarainen et al., 2016). In other words, NO & is a good proxy for combustion activity. We use a reduce major axis regression (Smith, 2009) to estimate 35 the slopes (∆ ∆ ) representing enhancement ratio across the spatial extent of the megacity and intercept ( *+ ) for CO and SO 2 representing the background levels when there is no combustion (within the megacity and free-tropospheric contribution). This follows the approach introduced by (Fujita et al. (1992)  and Parrish et al. (2002). However, we note that we use the spatial covariations of these species relative to NO 2 rather than their temporal covariations as in previous studies. Please see Section 3 for implications of this approach. We only consider statistically significant and positive slopes as we are focusing on sources and not sinks of these combustion products. These monthly ratios are then average across the year for analysis and archived for time series (decadal) analysis (see Section 3). Note that they can be 5 considered to be comparable to emission ratios when observations are taken at or near the source and if they are normalized to account for air mass variations (Fujita et al., 1992;Parrish et al., 2002;Parrish, 2006;Hassler et al., 2016). Here, we treat emissions of these species across the entire extent of the megacity as a point source and normalized all ratios to year 2005 values.
It is important to note that we view each large city as a big smokestack that emits an aggregate of 10 combustion products that can then be observed by satellite remote sensing as column-integrated quantities. The spatial (0.1-degree) covariation of these aggregate within the 2-degree radius is interpreted as bulk characteristic of spatially heterogeneous combustion sources within the megacity. Monthly enhancement ratios are hence interpreted as the linear sensitivity in CO or SO 2 to intra-megacity spatial variations in combustion activity as defined by NO 2 . We emphasize that these enhancement ratios are not derived using 15 time covariations but spatial covariations to minimize potential non-stationarities (e.g., differences in lifetimes between species) and influence of free-tropospheric signatures in MOPITT CO, which should be reflected as part of a larger scale contribution to CO *+ in this analysis given that we anchor the regression on OMI NO 2 . Possible confounding factors such as biogenic sources of CO in a megacity is also minimized in our analysis by treating CO data only when NO 2 is observed since NO 2 is not largely 20 co-emitted from CO biogenic sources. Although spatial and temporal smoothing can minimize the effect of lightning (NO X ) and fires (NO X and CO) since they are emitted intermittently relative to anthropogenic combustion, our findings must be interpreted to represent changes in bulk combustion cleanness over a megacity rather than specific combustion cleanness.

Time Series Analysis and Curve Fitting 25
The focus of this work is to study the long-term changes in the spatial covariations of these monthly-averaged CO and SO 2 to NO 2 as expressed in terms of enhancement ratios. We hypothesized that at decadal scale the changes in covariations reflect the dominant changes in megacity emission characteristics. We use two approaches to calculate the decadal trend in our normalized estimates of these ratios. For linear trend analyses, we use the Robust Regression Using Iteratively Reweighted Least-30 Squares (Holland and Welsch, 1977). This minimizes the influence of outliers relative to traditional leastsquares fit especially when the relationship is not fully linear. We also use another trend analysis algorithm in our subsequent inverse analysis. Instead of using the annual mean values and estimate the linear trend across 2005-2014, we estimate the associated decadal trends in ∆CO/∆NO & and ∆SO & /∆NO & using the seasonal trend decomposition with LOESS (locally weighted scatterplot smoothing) or STL 35 algorithm (Cleveland et al., 1990). This algorithm separates the seasonal, inter-annual, and decadal contributions of monthly ratios. We use the smoothing windows for the decadal, inter-annual, and seasonal trends of 121 months, 25 months, and 5 months, respectively based on analysis of CO decadal trends in Jiang et al. (2018). As in Gaubert et al. (2017), we tested several other windows and found consistent temporal patterns across cities. For non-linear curve fitting, we use robust least square 40 regressions with Least Absolute Residuals (LAR) method (within the cftool function in MATLAB) to fit a power law function to the annual-mean ratios and GDP per capita. This method also minimizes the influence of extreme values on the fit.

Inverse Analysis
We conduct an inverse analysis of the long-term trends in monthly enhancement ratios to further 5 expound our findings by associating the overall changes to sectoral changes. In this case, we are interested in finding the decadal contribution of the time series (2005-2014) of monthly statistically-significant enhancement ratios derived from our previous regression and time series analysis. We decomposed the a priori estimate of monthly emission ratio of CO to NO X (and SO 2 to NO X ) from RCP8.5 as a product of: a) ratio of effective emission factors for each of the four sectors namely energy, industry, transport, and 10 others; and b) fractional contribution of NO 2 emissions from each sector to the total NO 2 emissions for all four sectors. We then use a two-step Monte-Carlo-based Bayesian inversion method to estimate effective emission factors and fractional contribution of NO 2 emissions from each sector. Please refer to Appendix A for a short derivation of this decomposition, and Appendix B for details in the inverse analysis. 15

Observed Patterns of Enhancement Ratios in Chinese and U.S. Cities
Shown in Figure (Chan and Yao, 2008). In particular, combustion-related activities in Shenyang, Beijing, 20 Shanghai, and Shenzhen can be characterized to follow a progression from heavy to light manufacturing, export processing, and service industries. For this analysis, Shenyang, Beijing, Shanghai, and Shenzhen represent the progression across the 12 select cities of increasing GDP per capita along with decreasing to increasing ∆CO/∆NO & (-5.4±0.7% to +8.3±3.1%) and decreasing rate of ∆SO & /∆NO & reductions (-6.0±1.0% to -3.4±1.0%) relative to 2005 ( Figure 3 and Table 2). This pattern in enhancement ratios is not 25 evident in the rate of change of CO, SO 2 , and NO 2 column abundance, for which we find increasing rate of decrease in CO (-0.1±0.3% to -1.0±0.2%) and SO 2 (-1.9±0.9% to -5.5±1.1%) abundance along with decreasing rate of increase in NO 2 abundance from Shenyang (+5.2±1.4%) to Shenzhen (1.8±0.7%) ( Table 2) consistent with previous studies of these species. In fact, we find a decreasing-to-increasing pattern in the derived enhancements of CO due to combustion (i.e., ∆ ./0* = CO − CO *+ across these 30 four levels of development. This temporal pattern is reasonably robust as these CO enhancements are dominantly from combustion-related processes that co-emit NO 2 by our study design. Although naturally-produced CO and NO 2 like biogenic CO and lightning NO X introduce a strong seasonality on these ratios even within the megacity, we find that when we average the monthly ratios using only the months corresponding to a 35 particular season (i.e., more fires and lightning during the summer), we still find a similar temporal pattern (albeit different in magnitude) in derived ∆CO/∆NO & (Figure 2). We have minimized the influence of inter-annual variations due to meteorology (e.g., changes in air mass) by analyzing molar ratios (e.g., mole CO/mole NO 2 ) rather than absolute molar concentrations (e.g., mole CO/mole air). Normalizing these ratios to 2005 values should have also minimized the impact of the differences in the magnitude of these ratios between these cities. The impact of meteorology on inferred decadal trends through variations in columnar abundance is more evident when absolute magnitudes of single species are analyzed. In 5 addition, potential drifts of biases in time caused by systematic errors in the instrument and/or retrieval algorithm cannot account for the differences in the temporal pattern that we find across these cities. Such biases should be commonly reflected in all cities yet we see differences between cities. In fact, we find very similar progression pattern when we use the Infrared Atmospheric Sounding Interferometer (IASI) CO retrievals (De Wachter et al., 2012) instead of MOPITT or OMI QA4ECV (Boersma et al., 2017) 10 instead of OMI DOMINO. Interestingly, we find that the increasing enhancement ratio of CO to NO 2 in Shenzhen (and to a lesser extent in Shanghai) remarkably resembles the relative changes in CO to NO 2 ratios in more developed megacities (Los Angeles and New York) and several urban agglomerations in the United States (see Figure 3a for Los Angeles and Table 2 and Figure S1 for all other select cities). More importantly, the increasing pattern that we see in Los Angeles (~ +7±1%) relative to 2005 is 15 generally consistent to the increasing trend (~ +4%) after 2007 of ground-based CO to NO X enhancement ratio in Los Angeles as reported by Hassler et al. (2016). While it is a common understanding that modernization brings about larger energy use coupled with higher economic productivity but poorer environmental quality (i.e., increasing abundance of pollutants), the changes in lifestyle concomitant with human development results in a shift to fewer activities (including increase use of renewable energy) 20 along with more efficient and cleaner combustion and changes in fuel types (coal to natural gas) (Mazur and Rosa, 1974). This eventually leads to increases in relative sensitivities of CO and SO 2 to NO 2 . Along the same line as previous studies suggesting emissions of CO, SO 2 , NO 2 , and their ratios can be indicators of modernization to some extent (Krotkov et al., 2006;Russell et al., 2012;Luo et al., 2014;Hassler et al., 2016) Unlike ∆ ./0* , ∆ & ./0* in all four Chinese cities still show a decreasing trend relative to 2005 while ∆ & ./0* in Los Angeles show an increasing pattern consistent with its ∆ ./0* . On one hand, there is a striking difference in absolute magnitudes in SO 2 abundance between these cities (as has been reported), reflecting large-scale differences in combustion practice. Yet, the low SO 2 abundance in Los Angeles makes it also difficult to detect possibly large SO 2 point sources (Krotkov et al., 2016). Enhanced SO 2  35 signal can still be detected as the spatial first-order derivatives of SO 2 with NO 2 at megacity-scale should not be largely (non-linearly) influenced by its absolute magnitude. We find that there is a tighter correspondence between SO 2 and NO 2 abundance in Chinese cities than in U.S. cities, which might suggest differences in fuel use as SO 2 is mainly produced within a megacity from burning of sulfurcontaining fossil fuel (mostly coal, oil, and natural gas) and to a smaller extent from industrial processes 40 cities is due to continuing heavier reliance of these cities (and China) on coal burning relative to United States (Wang and Hao, 2012;Qi et al., 2016;Yang et al., 2016). In 2005, it was estimated that coal accounts for about 69% and 23% of the total primary energy consumption in China and U.S., respectively. While there are on-going activities regulating coal-related emissions, including usage of low-sulfur coals, installation of flue gas desulfurization (FGD) facilities, and closing of small units, coal consumption in 5 China remains to increase in the past decade (Qi et al., 2016;Yang et al., 2016). In terms of mass, it has increased by 70% from 2005 to 2014 (Korsbakken et al., 2016). On the other hand, the use of coal in U.S. has been found to be slightly decreasing along with previous adoption of SO 2 control technologies (Taylor et al., 2005).

Inconsistencies with A Priori Estimates 10
The satellite-based ∆CO/∆NO & patterns are inconsistent with emission-and model-based ratios (Figures 3 and S1, Table 2). As previously introduced, estimates of the ratios of emissions can be related to observed ratios of enhancements when these observations are taken at or near the source. In this case, we assume that a megacity is a big smokestack emitting mostly combustion-related pollutants (i.e., CO, NO 2 , and SO 2 ) that can be observed from space with MOPITT and OMI. In addition, NO 2 is considered 15 to be the dominant form of NO X that can be observed at this scale. From a global atmospheric chemistry modeling (CTMs) perspective, the associated abundance over megacities is represented as one to four discrete vertical column(s) assuming spatial resolution of these CTMs of one to two degrees. While recognizing the associated month-to-month variability in ∆CO/∆NO & and expected differences on how these ratios should be compared, the trends in emission ratios relative to 2005 of CO to NO X from bottom-20 up emission inventories (EDGAR4.2 and RCP8.5) and top-down emission estimates (CHASER, Miyazaki et al., 2017) do not appear to follow the progression from Shenyang to Shenzhen of decreasing to increasing ∆CO/∆NO & relative to 2005 ( Figure 3). This is also true for the ratios of CO to NO 2 abundance from CAM-Chem and CHASER CTMs, which are mostly consistent (except in Los Angeles) with the trends of their associated emission ratios (i.e., CAM-Chem and CHASER emissions are based 25 on RCP8.5 and EDGARv4.2 inventories, respectively). The a posteriori emission ratios in Beijing from Miyazaki et al. (2017), which uses CHASER-LETKF to assimilate MOPITT CO and OMI NO 2 retrievals among other retrievals, also appear to initially follow the emission ratios from EDGAR. Furthermore, the ratios of SO 2 to NO X emissions from RCP8.5 follow the trend of ∆SO & /∆NO & in Chinese cities but tend to diverge in Los Angeles, whereas the emission ratios from EDGAR exhibit a lack of trend in China and 30 Los Angeles. A closer look at linear trends of the ratios for each sector in RCP8.5 ( Figure S2) reveals inconsistencies in the trends which cannot be addressed by simple scaling of activity levels in bottom-up inventories (Zheng et al., 2018). All these differences underscore the need to reduce uncertainties in representing time-varying emission activity and emission factors in CTM inputs. There is also a need to quantify errors in model physics and dynamics in transforming emissions to abundance, as well as in data 35 assimilation and inverse methods in integrating observations into models including representativeness of these retrievals. We highlight here the need to improve not only the accuracy but also the consistency of AQ predictions across pollutants in megacities. Initial results from an improved set of multi-species data assimilation runs using CHASER-LETKF show better agreements with the trends in ∆CO/∆NO & (Miyazaki et al., 2017). Such improvements highlight an under-explored utility of available observational 40 constraints on the changes in emission ratios. We emphasize here that while these differences are expected and have been previously reported, our findings highlight the need to focus on improving model treatments of the dynamic nature of emission factors in these megacities.

Combustion Emission Pathway for Chinese Cities
We conducted an inverse analysis of the ratios shown in Figure 3 to further expound on these 5 patterns by associating them to sectoral changes and identifying a common combustion emission pathway across these four levels of development. Please see details of the matrix-vector product and inversion methodology in Section 2.5 and Appendix B. The result of this inversion is a set of a posteriori time series estimates of sectoral CO to NO X and SO 2 to NO X ratios, such that the corresponding time series estimates of the total CO to NO X and SO 2 to NO X ratios match the decadal trends of ∆CO/∆NO & and 10 ∆SO & /∆NO & inferred from these satellite retrievals. Again, we note that we use the STL-inferred decadal trend as the data to fit (not the monthly-mean ratios nor the linear trend in Figure 3), as this is the most appropriate data for analyzing long-term changes in emission sectors. The results of our inverse analysis are presented in Figure 4. This figure consists of five 2-D line plots of a posteriori (solid) and a priori (dashed) time series of SO 2 to NO X emission ratios (ESO 2 /ENO X ) 15 in y-axis versus corresponding values of CO to NO X emission ratios (ECO/ENO X ) in x-axis. The five plots correspond to the annual total (center panel, Figure 4a) and sectoral emission ratios (four side panels, Figures 4b to 4e) of each of the four cities selected in Figure 3. The time series, which is normalized to 2005 values, starts at the origin (1,1) and ends at the arrow tip of the line. Each 2-D plot also contains an inset showing the corresponding emission trajectory for Los Angeles. The center panel of Figure 4 is 20 similar to Figure 3 but plotted jointly and with the a posteriori time series of emission ratios now corresponding to the time series of enhancement ratios (i.e., STL-inferred decadal trend). We find that the progression in combustion characteristics across these four cities is clearly evident from this diagram and very consistent with the linear trends in Figure 3. In Shenyang, both ESO 2 /ENO X and ECO/ENO X are decreasing relative to 2005 at a faster rate (as represented by the length of the line) than in Beijing. 25 On the other hand, we see a clear shift in Shanghai and most notably in Shenzhen to a slightly decreasing ESO 2 /ENO X and increasing ECO/ENO X leading their emission trajectories toward a different state of 'combustion cleanness'. The combustion emission ratios in Los Angeles (and other cities in U.S.) lies however at a different state than Shanghai and Shenzhen. In particular, we find ESO 2 /ENO X and ECO/ENO X in Los Angeles to be both linearly increasing relative to 2005 values. And so, there exists a 30 progression of decreasing-to-increasing sensitivities of CO and SO 2 to NO 2 from Shenyang to Shenzhen to Los Angeles (gray semi-circular trace in Figure 4a) relative to 2005 that appears to be related to socioeconomic development consistent with the current understanding of human development pathways (Lamb et al., 2014), which in our case may be a consequence of air quality management practice and improved efficiency in China (Sun et al., 2018; van der A et al., 2017) and U.S. (Hassler et al., 2016;35 Russell et al., 2012). Altogether, this leads us to suggest a common combustion emission pathway for the megacities in mainland China that begins with a reduction in SO 2 , followed by CO, and continues with a reduction in NO X and potentially on volatile organic compounds (VOCs) later on. To illustrate, we still see increases in NO X abundance in Shenyang although CO and SO 2 are already decreasing, whereas in Shenzhen, we see NO X starting to decrease (at a faster rate) along with decreasing CO and SO 2 abundance. The rate at which SO 2 , CO, and NO 2 are decreasing is not at a level that is observed in Los Angeles. And so, while the satellite data reveals a combustion emission pathway in these Chinese megacities, these cities are yet to reach conditions that is at par with megacities in more developed cities in U.S. and Europe. It is worth noting that the a priori estimates from RCP8.5 do not follow this pathway, even for Los Angeles, suggesting inconsistencies and necessary updates on temporal changes in emission factors, 5 effectiveness in pollution control technologies, and/or more information on fuel use mixtures in this emission inventory. It also appears that the pathway represented in RCP is similar to all cities and more resembling the emission pathway for Beijing. Furthermore, the traces in sectoral emission ratios from RCP8.5 all point to decreasing ratios relative to 2005 and are primarily driven by the energy (transportation) sector, which constitute more than 10 one-third of NO X emissions in Chinese (U.S.) cities (Figure 4b to 4e). Our inversion results to slight adjustments in Chinese energy emission pathway towards little to no changes in CO to NO X emission ratios ( Figure 4b). Adjustments from the transportation sector are also small in terms of direction and slower in terms of its rate of change relative to 2005 RCP values (Figure 4c). This is certainly not the case in Los Angeles where CO to NO X and SO 2 to NO X ratios follow quite the opposite pathway of increasing 15 ratios from the energy sector and increasing CO to NO X with no change in SO 2 to NO X from the transportation sector. This is expected in United States because of cleaner fuel standards (REF).
Significant shifts on these ratios relative to 2005 are clearly evident from the industry and other (i.e., agriculture, residential, and waste) sectors in the cities in China (Figure 4d and 4e). Shanghai and most notably Shenzhen show a shift to increasing CO to NO X with slightly decreasing SO 2 to NO X that are not 20 reflected in RCP8.5. The emission ratios from industry and other (mostly residential) sectors need to be adjusted significantly in our inversion to match the shifts in observed ∆CO/∆NO & and ∆SO & /∆NO & in these two cities. As earlier mentioned, tertiary (service) industries including export processing activities are dominant in Shanghai and Shenzhen than in Shenyang and Beijing. The shift in recent years to increasing CO to NO X reflects a larger rate of decrease in NO X levels than CO from the industrial and 25 residential sectors of these cities. While a more detailed investigation is warranted to narrowly identify the activities and/or policies driving this shift (van der A et al., 2017), it is clear that changes in combustion activity alone cannot account for these shifts and that updates on emission factors for these sectors in RCP8.5 are needed. We find that these findings are robust across a suite of error assumptions in the inverse analysis. This update applies all the more to all sectors in RCP emissions for Los Angeles. Again, 30 this is well supported by studies like Hassler et al. (2016). where they reported increasing CO to NO X enhancement ratio after 2007 in Los Angeles along with a 45% decline of NO X emissions based on their fuel-based inventory. This is in contrast to decreasing RCP8.5-based MACCity emission ratios that they also reported for Los Angeles. This increase in enhancement ratios (similar to this work) is attributed to a combination of factors such as the decrease in NO X from freight traffic activity during U.S. recession 35 and implementation of new NO X emission control technologies and regulations to meet Tier two emission standards on U.S. light-duty vehicles. They also noted that differences in the trends of ∆CO/∆NO & are still observed even between cities from developed countries like U.S. and Europe, as these cities differ in terms of transportation practices and lifestyles (e.g., increase in light duty diesel vehicles). It is also now conceivable that ∆CO/∆NO & can be further influenced by shifts in relative importance of emission sectors 40 (e.g., VOCs in petrochemical and pharmaceutical industries) as activity decreases with efficiency, pollution is controlled, and lifestyle changes whenever cities evolve . A recent study (Jiang et al., 2018) revealing an over-estimation in the decrease of USEPA NO X emissions based on OMI NO 2 and MOPITT CO retrievals with USEPA ground station measurements of NO 2 , also suggests potential changes in 'bulk' combustion characteristics in urban regions of the United States. Along with these studies, our results suggest that regional to global emission inventories, which are used as input to predictive models of atmospheric composition, have to reflect: a) the evolution of air pollution 5 for a given city (sectoral shifts) and b) the differences in combustion practices from city to city, in order to capture these observed magnitudes and variations in enhancement ratios.

Socioeconomic Dependence of Urban Enhancement Ratios in China
Here, we attempt to connect these emission pathways to the larger pattern of economic growth across the 31 capital cities and five special cities in mainland China. We find in particular a power law 10 relationship between the observed annual-mean ∆CO/∆NO & (and ∆SO & /∆NO & ) and GDP per capita. This is not to derive an overall EKC for China, as this in fact requires a very long record of environmental quality, but specifically to investigate how economic development shapes how 'clean' the bulk combustion in Chinese cities would be. These enhancement ratios complement abundance and/or emissions of pollutants as traditional measures of air pollution. Unlike Figure 3 and 4, our focus is to 15 illustrate the larger dependence of enhancement ratios on GDP per capita. As discussed in the Methods section, we relate the enhancement ratio of a megacity to the ratio of the product of emission factor ( 678.986 ) and effectiveness of control technology (1 − 678.986 ) for CO and NO X species in the case of ∆CO/∆NO & for example. We use a robust least-squares regression with least absolute residuals method to fit a curve of the form: = < , where is ∆CO/∆NO & or ∆SO & /∆NO & and is GDP per capita. Our 20 results are presented in Figure 5a and 5b for ∆CO/∆NO & and ∆SO & /∆NO & , respectively. The 12 cities considered in our analysis of emission pathways are marked with colors corresponding to its level of urban development described in previous section. Note that the magnitudes of enhancement ratios derived from this work is a factor of 10 higher than ratios derived from ground-based networks. We attribute this discrepancy to differences in air mass and volume, representativeness, and vertical sensitivity between 25 abundance retrieved as total or tropospheric columns and in-situ and point samples in units of mixing ratios. Nevertheless, we find a strong power law relationship with GDP per capita having coefficients (R 2 =0.98) of negative two-thirds and negative one-half for ∆CO/∆NO & and ∆SO & /∆NO & , respectively. Likely, the coefficients in ∆SO & /∆NO & will converge to that in ∆CO/∆NO & as changes in fuel type and SO 2 controls should decrease SO 2 abundance. While each city is unique and that the evolution of air 30 pollution may be different from city to city, there also exist a clear signature of urbanization at national level that reflects the influence of economic growth on the cleanness of bulk combustion. Similar power law relationships (albeit different coefficients) have been reported in studies of urban growth and development (Bechle et al., 2011;Lamsal et al., 2013;Bettencourt et al., 2013), energy flows (Creutzig et al., 2015) and carbon emissions (Fragkias et al., 2013). Our results suggest that enhancement ratios 35 scale with GDP per capita, with lower GDP per capita like Shenyang and other cities (gray dots) having higher enhancement ratios, while Shenzhen and other cities (yellow dots) with highest GDP per capita in China lie among cities with the lowest enhancement ratios. As we have shown in Figure 4 (and Table 2), the ratios in Shenzhen tend to increase with time (and GDP) but this increase has its limits and appears to be dwarfed by cities with highest enhancement ratios. We note, however, that identifying a mechanistic 40 Atmos. Chem. Phys. Discuss., https://doi.org/10.5194/acp-2018-1121 Manuscript under review for journal Atmos. Chem. Phys. Discussion started: 4 December 2018 c Author(s) 2018. CC BY 4.0 License. rationale of these negative scaling coefficients is beyond the scope of this work and hence is not proposed. A unified relationship cannot also be established across countries as there are obvious differences in socioeconomic and air pollution conditions in China and U.S. that cannot be accounted for ( Figure 6). Nevertheless, we suggest incorporating this observable along with estimates of emissions to future scaling studies, especially as we move past RCPs and toward recent developments in building more realistic 5 emission scenarios that integrate socioeconomic and environmental development pathways like the Shared Socioeconomic Pathways (SSPs; O'Neil et al., 2014).

Summary and Implications
The main goal of this work is to provide observational evidence from Earth observing satellites of emission pathways of combustion-related air pollutants as a result of urban growth in economically 10 developing countries like China. A new observational perspective on monitoring one of the major consequences of urbanization is introduced, not to replace existing observing capabilities but to further exploit the information that is already available. Following the pioneering work by Parrish et al. (2002), the sensitivities of intermediate products of combustion can be derived from existing satellite retrievals of air quality (AQ) to inform changes in bulk combustion characteristics (and consequently emissions) of 15 a megacity. This is especially relevant as the number of megacities continue to grow in the coming decades, mostly at locations that lack sufficient AQ monitoring capabilities. Enhancement ratios of CO to NO 2 and SO 2 to NO 2 over megacities in mainland China that are derived from MOPITT and OMI satellite instruments show a coherent long-term progression in recent years of decreasing to increasing ratios relative to 2005 that is well correlated with economic development. These trace a common emission 20 pathway that resembles the evolution of air pollution in more developed cities in the United States which is characterized by transitions in energy use and subsequent implementation of pollution control and regulation. Although we find cleaner combustion as cities in China develop consistent with their Five Year Plans, this is presently obfuscated by increasing fuel use particularly its heavy reliance on coal. We propose the use of these enhancement ratios derived from existing satellite retrievals to complement 25 existing surface AQ networks, including carbon-related satellite observing systems in further constraining combustion efficiency and effectiveness of control technologies and policies. Augmenting existing capabilities (Saeki et al., 2017) is particularly relevant especially with the aid of big data informatics and machine learning as well as the advent of activities focusing specifically on tracking fossil fuel emissions like the CO 2 Human Emissions project (https://www.che-project.eu). While we recognize the current 30 limitations of these retrievals (e.g., collocation, sensitivity), our findings appear to be robust across retrievals and methods and supported by previous studies using these retrievals in a different way (Krotkov et al., 2016;Jiang et al., 2018) or ground measurements (Hassler et al., 2016). We strongly suggest that the capability to monitor relatively long-term changes in atmospheric composition has to be supported and continued with complementary new satellite and field missions and deployments (Streets 35 et al., 2013; National Academies of Sciences, Engineering, and Medicine, 2016). The relative importance of monitoring combustion efficiency and effectiveness of pollution control increases as a city and country continue to socioeconomically develop and become sustainable. Despite past and present studies (Mazur and Rosa, 1974;Lamb et al., 2014), it is only in most recent years that we have developed comprehensive and integrated monitoring and prediction systems which paved new measures of air pollution and new developments in emission scenarios like SSPs. For China, more detailed information on energy use and improved emission inventories are increasingly becoming available for assessment (Li et al., 2017;Zhong et al., 2017). As we also recognize some of the challenges to quantify socioeconomic variables such as the impact of international trade on air pollution (Lin et al., 5 2014), economic structural upgrading (Mi et al., 2017), greater utilization of renewable energy, and even metrics of performance (Ramaswami et al., 2013), from a physical science perspective, our results strongly support these new developments. We find inconsistencies between the long-term spatiotemporal patterns of emission ratios from RCP8.5 and model predictions of abundance ratios and the corresponding patterns derived from observed enhancement ratios. Scientific improvements in representing the evolution 10 of air pollution (Lewis, 2018) and emission pathways (Mitchell et al., 2017) can be made by considering observationally-constrained time-varying emission factors and confronting emissions and physical models with available data not only for their accuracy but also for their consistency in representing both carbon and AQ-related combustion products. 15 Data availability. The raw data used in this study are available online (links to satellite data and emission inventories can be found in Table 1; Socioeconomic Data: Annual GDP and population are directly taken from China Statistical Yearbook compiled by National Bureau of Statistics of China (http://www.stats.gov.cn/)). Model outputs and reanalysis data are available upon request from the authors. 20 Acknowledgments. This study is supported by NASA ACMAP Grant NNX17AG39G. K.M's reanalysis is supported by JSPS KAKENHI Grant 15K05296 and 18H01285. We acknowledge MOPITT, IASI, and OMI retrieval teams for CO, NO 2 , SO 2 , data, respectively. We also thank EDGAR, HTAP, REAS, RCPs data teams for the emission inventories. All the satellite data and emission inventories are available to the public online. We thank Kevin Bowman, Cenlin He, and Sam Silva for insightful discussions. 25 Author contributions. The initial idea was provided by AFA. and WT. CHASER-LETKF experiments were performed and provided by KM. CAM-chem modeling experiments were performed and provided by BG. Data analyses and results interpretation were performed by WT. and AFA. HW provided key expert guidance on MOPITT CO. The manuscript was written by AFA. and WT. 30

Appendix A. Combustion Emission Ratios and their Decomposition
In a combustion process using a hydrocarbon fuel, CO and elemental carbon (e.g., soot or BC) are 35 produced when combustion is incomplete; otherwise carbon in the fuel is oxidized to CO 2 (Eq. 1). In addition, NO and NO 2 are produced from the oxidation of nitrogen from the fuel itself and from decomposition of N 2 in air at high temperatures (Flagan and Seinfeld, 2012). Sulfur dioxide (SO 2 ) is also produced when the fuel used in the combustion process contains sulfur (such is the case for low-grade fuels). 40 Atmos. Chem. Phys. Discuss., https://doi.org/10.5194/acp-2018-1121 Manuscript under review for journal Atmos. Chem. Phys. where ? is the total mass of emissions for species , ?,6 is its associated emission factor for a specific source/sector , 6 is the activity level of the source. ?,6 corresponds to effectiveness of control 10 measure and ?,6 = ?,6 • 1 − ?,6 is the effective emission factor. When we take the ratio of emissions (Eq. 2) of co-emitted species and , this ratio can be expressed as the sum of the products of the ratio of effective emission factors ( ?,K,6 QQR ) and the fractional contribution of emission sector ( ?,6 ) (Eq. A3). 20 We decomposed the a priori estimate of monthly emission ratio of CO to NO X (and SO 2 to NO X ) from RCP8.5 as a product of: a) ratio of effective emission factors for each of the four sectors namely energy, industry, transport, and others ( ?,K,6 QQR ); and b) fractional contribution of NO 2 emissions from each sector to the total NO 2 emissions ( ?,6 ) for all four sectors s ( T : energy, & :industry, U :transport, V :others). In matrix-vector form, this can be expressed as: 25 decadal trends of the RCP8.5 CO to NO X and SO 2 to NO X emission ratios using the decadal trends of ∆CO/∆NO & and ∆SO & /∆NO & as observational data ( ). We use the decadal trend of enhancement ratios of CO to NO 2 and SO 2 to NO 2 (derived using STL), calculate their annual averages and normalized to 2005 values, and then take these as our observational (fitting) data. Our goal is to estimate and given subject to the following constraints: a) errors in and are 10% and 25% of their values, b) errors in 5 is 5% of its value, c) error covariances of and are uncorrelated and diagonal ( , ) and d) sum of is unity. Since this is an under-determined inverse problem, we apply prior information on and using the RCP emissions ( , ). We conduct our inverse analysis into two-step: 1) estimate the most likely that results to estimates of best fitting the decadal trend, 2) estimate using the new estimate of . For Step 1, first, we draw n=10,000 samples of assuming its errors are normally distributed with 10 mean to be its prior and covariance to be the diagonal of its squared errors. Second, we use the maximum a posteriori (MAP) solution to the Bayesian problem to estimate for every sample. i.e., We draw a new sample if any of the elements in is negative. Third, we take the mean of 100 15 samples resulting to the lowest root-mean-square errors relative to the data. We use this mean as our new estimate of ( ). For

Appendix B. Inverse Analysis
Step 2, we apply the same MAP solution using and = to estimate and . Similar to a Kalman filter, we cycle this procedure for each year starting from 2006 to 2014. We use the new estimates of , , and for a given year as priors for the succeeding cycle with fix inflation on the covariance of 1.25 to minimize filter divergence. We note that the additional constraints (positive , 20 sum of is unity) minimizes the underdeterminacy of the problem. This is supported by post-inverse analysis diagnostics (i.e., averaging kernels) showing that elements of are resolved by the trend data.