Satellite data reveal 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. In order to estimate the trends in enhancement ratios of CO and SO₂ to NO₂ (ΔCO∕ΔNO₂ and ΔSO₂∕ΔNO₂) and infer emergent bulk combustion properties over these cities, we combine air quality retrievals from widely used satellite instruments across 2005–2014. We present results for four Chinese cities (Shenyang, Beijing, Shanghai, and Shenzhen) representing four levels of urban development. Our results show a robust coherent progression of declining to growing ΔCO∕ΔNO₂ relative to 2005 ($-\mathrm{5.4}\pm \mathrm{0.7}$ to $+\mathrm{8.3}\pm \mathrm{3.1}\mathit{\%}$ yr⁻¹) and slowly declining ΔSO₂∕ΔNO₂ ($-\mathrm{6.0}\pm \mathrm{1.0}$ to $-\mathrm{3.4}\pm \mathrm{1.0}$ % yr⁻¹) across the four cities. The coherent progression we find 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 that is not yet seen in Shenyang and Beijing. This overall trend is presently 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 in emission ratios in current inventories. As cities and/or countries continue to socioeconomically develop, the ability to monitor combustion efficiency and effectiveness of pollution control becomes increasingly important in assessing sustainable control strategies.

1 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, 2015). Anthropogenic activities are most intense in megacities, 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 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 prevent us from accurately assessing the changes in atmospheric 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 in substantial emissions (Richter et al., 2005) and more frequent occurrences of severe pollution events in many of its megacities, most notably Beijing (Guo et al., 2014). These affect not only local AQ and public health but are also reported to impact the 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. This points to important changes in AQ as a result of development, AQ management, and regional to national socioeconomic 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 the city scale remains a 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 ratios of these products) (Silva and Arellano, 2017; Silva et al., 2013; Hassler et al., 2016; Tang et al., 2018, 2019), alongside information on the state of socioeconomic 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), sulfur dioxides (SO₂), and soot particles with final products like carbon dioxide (CO₂). Please see the 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 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; Ding et al., 2017); to indicate bulk combustion efficiency, inter-megacity differences, and fire phase (Silva et al., 2013; Silva and Arellano, 2017; Tang and Arellano, 2017); and to infer fossil fuel CO₂ 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 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 large cities across the recent decade (through these ratios), associate these shifts with particular sectors, and identify a common emission pathway across these cities. Along the same line of 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; Bechle et al., 2011; Zhang et al., 2012; Wang and Hao, 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 timescales. As a basis for comparison, pollution controls adopted in developed countries like the United States and Europe, which followed a progression from first controlling SO₂, 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.

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₂ from Ozone Monitoring Instrument (OMI), and planetary boundary layer (PBL) columns of SO₂ from OMI to derive monthly estimates of these ratios. We conduct 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., 2011). We also conduct 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 Sect. 3. Section 4 is summary and implication of this study.

2 Data and methods

2.1 Study region

We consider 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 agglomerations in the country (see Fig. 1 for coverage). For purposes of finding long-term emergent patterns on its emission characteristics, we focus 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 the National Bureau of Statistics of China (http://data.stats.gov.cn, last access: 22 March 2019), these cities already exhibit largely diverse pollution and economic development attributes illustrated in Fig. 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 consider Los Angeles and other large cities in the United States (New York City, Chicago, Houston, Phoenix, Boston, Seattle, and Miami) for comparison.

Table 1 List of satellite products and emission inventories used in this study. All these datasets are re-gridded into $\mathrm{0.1}{}^{\circ }\times \mathrm{0.1}{}^{\circ }$ if not originally at this resolution. This version of CHASER–LETKF does not provide emissions of SO₂.

^* Last access date: 22 March 2019.

Figure 1 Time series (2005–2014) of RCP8.5 combustion-related emissions of NO_x (first quad), CO (second quad), and SO₂ (third quad) all in units of grams per year per square meter and GDP per capita (fourth quad) in units of CNY 10⁵ per capita per year for each of the 12 select major cities (red dots) in mainland China. The scales of each quadrant are indicated in the legend (lower left of the map). The total emissions for each combustion product are broken down into four major sectors: energy, industry, land transport, and others, which is the sum of agriculture, residential and commercial, and waste treatment and disposal). The GDP per capita is also broken down into primary (direct use of natural resources), secondary (industry and manufacturing), and tertiary (service) sectors. Each blue dot corresponds to one of the 36 designated provincial capital and special cities in mainland China.

2.2 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.

2.2.1 Satellite retrievals and data processing

We use the NASA Terra Measurement of Pollution In The Troposphere (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 the Ozone Monitoring Instrument (OMI) planetary boundary layer (PBL) SO₂, version 3, Level 2 (Krotkov et al., 2006). We collect daily MOPITT CO, OMI NO₂, and OMI SO₂ retrievals that are available within a $\mathrm{2}{}^{\circ }\times \mathrm{2}{}^{\circ }$ area around each city center. This radius is selected to cover the extent of each city based on NO₂ footprints (Bechle et al., 2011; Lamsal et al., 2013) and geopolitical maps of city boundaries. We grid each set of retrievals into $\mathrm{0.1}{}^{\circ }\times \mathrm{0.1}{}^{\circ }$ grids that commensurate to the finest retrieval resolution among MOPITT and OMI. We then average them across each 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₂, NO₂) 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 (e.g., Krotkov et al., 2016). While CO retrieved from thermal infrared (TIR) radiances is mostly sensitive to free tropospheric CO, it has also been reported to be capable of observing lower tropospheric CO, especially when retrieved jointly from TIR and near-infrared (NIR) radiances (Worden et al., 2010; Deeter et al., 2014). We recognize, however, that retrievals of SO₂ from OMI have been reported to exhibit low sensitivity to weak SO₂ signals, in particular to less than 30 to 70 kTon yr⁻¹ of point source emissions (Krotkov et al., 2016). While our spatial and temporal smoothing, along with anchoring our SO₂ analysis with NO₂ data (please see later description of our regression analysis), should help in enhancing the SO₂ signal from cities with low SO₂ emissions, these SO₂ retrievals are useful as large SO₂ abundances are still observed across the majority of cities in China (Krotkov et al., 2016). We also use CO retrievals from the Infrared Atmospheric Sounding Interferometer (IASI), Level 2 (De Wachter et al., 2012), and tropospheric column NO₂ from FP7 QA4ECV OMI v1 (Boersma et al., 2017) to verify consistency in our trend estimates.

Table 2 Summary of percent rate of change for select cities in China and the United States. Numbers that follow the ± sign are standard errors.

We note that using $\mathrm{2}{}^{\circ }\times \mathrm{2}{}^{\circ }$ area to represent cities does lead to slight overlap over Guangzhou and Shenzhen, Beijing, and Tianjin. This does not affect our analyses of emission inventories because we apply geopolitical maps of city boundaries to calculate emissions for each city (see Sect. 2.2.2). This does have an impact on our analyses of satellite observations because we use all the grids in the $\mathrm{2}{}^{\circ }\times \mathrm{2}{}^{\circ }$ area to conduct the spatial regression. However, we do not expect the overlap to significantly change our results because (1) the overlapped area is relatively small; (2) the overlapped cities are sometimes considered together as a whole region because of their similarities and connections (for example, the Jing-Jin-Ji megalopolis and the Pearl River Delta); and (3) the overlapped cities are in the same classes with similar patterns based on our analyses (i.e., Beijing and Tianjin are both in class 2, while Guangzhou and Shenzhen are both in class 4; Table 2).

2.2.2 Emission inventories and model simulations

Multiple bottom-up emission inventories for CO, NO₂, and SO₂ are analyzed, namely the Emission 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₂ from the Tropospheric Chemical Reanalysis (TCR) based on the CHASER–LETKF assimilation system (Miyazaki et al., 2017). Since these emission inventories have different spatial resolutions (see details in Table 1) and are available in the form of fluxes (units in kilograms per square meter per second), we upscale or 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 Fig. 2). All of the cities extend to an area less than the 2^∘ by 2^∘ that we set as our city domain.

Figure 2 Spatial regression analysis of satellite retrievals of CO and SO₂ to NO₂ by season (blue: March–May, MAM; red: June–August, JJA; green: September–November, SON; orange: December–February, DJF). The left column shows an example of scatter plots and linear regression for Beijing (a, b, c) and Los Angeles (d, e, f). The center column corresponds to the changes across 2005 to 2014 in the ratios calculated for a given season. The rightmost column panels show the city domain ($\mathrm{2}{}^{\circ }\times \mathrm{2}{}^{\circ }$) with the geopolitical extent of the city of Beijing and Los Angeles.

We also use model data for CO and NO₂ from the Community Atmosphere Model with Chemistry (CAM-chem; Gaubert et al., 2016) and TCR to derive CO and NO₂ 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 and NO₂ 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 analyze are far coarser than 0.1^∘, we only consider the associated abundance ratio rather than deriving enhancement ratio across the months when non-stationarity and nonlinearity issues are more likely to exist.

2.3 Deriving enhancement ratios using spatial regression analysis

For each city, we regress the gridded monthly-average CO and SO₂ to NO₂ to calculate monthly enhancement ratios (ΔCO∕ΔNO₂ and ΔSO₂∕ΔNO₂). We use NO₂ as our control variable as NO₂ has the shortest lifetime (hours) among these products. Except for lightning, NO_x is mostly produced from high-temperature anthropogenic combustion processes. And because of its short lifetime, it is observed as distinctly and spatiotemporally local surface enhancements, with relatively very low background concentrations. Along with the availability of NO₂ retrievals from satellites at fine spatial scales and over long periods, NO₂ 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 reduced major axis regression (Smith, 2009) to estimate the slopes (Δy∕Δx) representing enhancement ratio across the spatial extent of the megacity and intercept (y^bg) for CO and SO₂ 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₂ rather than their temporal covariations as in previous studies. Please see Sect. 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 averaged across the year for analysis and archived for time series (decadal) analysis (see Sect. 3). Note that they can be 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 normalize all ratios to 2005 values.

It is important to note that we view each large city as a big smokestack that emits an aggregate of combustion products that can then be observed by satellite remote sensing as column-integrated quantities. The spatial (0.1^∘) covariation in this aggregate within the 2^∘ radius is interpreted as a bulk characteristic of spatially heterogeneous combustion sources within the megacity. Monthly enhancement ratios are hence interpreted as the linear sensitivity in CO or SO₂ to intra-megacity spatial variations in combustion activity as defined by NO₂. We emphasize that these enhancement ratios are not derived using 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^bg in this analysis given that we anchor the regression on OMI NO₂. Possible confounding factors such as biogenic sources of CO in a megacity are also minimized in our analysis by treating CO data only when NO₂ is observed since NO₂ is not largely 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.

2.4 Time series analysis and curve fitting

The focus of this work is to study the long-term changes in the spatial covariations of these monthly-averaged ratios of CO and SO₂ to NO₂, as expressed in terms of enhancement ratios. We hypothesize that at a 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 robust regression using iteratively reweighted least squares (Holland and Welsch, 1977). This minimizes the influence of outliers relative to a traditional least-squares 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 estimating 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 seasonal-trend decomposition using LOESS (STL) algorithm (Cleveland et al., 1990). This algorithm separates the seasonal, interannual, and decadal contributions of monthly ratios. We use the smoothing windows for the decadal, interannual, and seasonal trends of 121, 25, and 5 months, respectively, based on analysis of CO decadal trends in Jiang et al. (2018). As in Gaubert et al. (2017), we test several other windows and found consistent temporal patterns across cities. For nonlinear curve fitting, we use robust least-square regressions with the 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.

2.5 Inverse analysis

We conduct an inverse analysis of the long-term trends in monthly enhancement ratios to further expound our findings by associating the overall changes with 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 that are derived from our previous regression and time series analysis. We decompose the a priori estimate of the monthly emission ratio of CO to NO_x (and SO₂ 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) and (b) fractional contribution of NO₂ emissions from each sector to the total NO₂ 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₂ emissions from each sector. Please refer to Appendix A for a short derivation of this decomposition and Appendix B for details of the inverse analysis.

3 Results and discussions

3.1 Observed patterns of enhancement ratios in Chinese and US cities

In this subsection, we present observed patterns of enhancement ratios in Chinese and US cities. We firstly show spatial regression analysis of satellite retrievals of the ratio of CO and SO₂ to NO₂ by season (taking Beijing and Los Angeles for demonstration) in Fig. 2. Although naturally produced CO and NO₂ like biogenic CO and lightning NO_x introduce a strong seasonality in these ratios even within the megacity, we find that when we average the monthly ratios using only the months corresponding to a 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₂ and ΔSO₂∕ΔNO₂ (see Fig. 2). This is reasonable as these CO as well as SO₂ enhancements are dominantly from combustion-related processes that co-emit NO₂ by our study design, pointing to the robustness of analyzing annual-mean ΔCO∕ΔNO₂ and ΔSO₂∕ΔNO₂.

Figure 3 Changes in annual-mean enhancement ratios (black) from MOPITT and OMI retrievals of ratios of CO to NO₂ (top) and SO₂ to NO₂ (bottom) for select cities in China and the US relative to 2005. Its associated emission ratios (ECO∕ENO_x and ESO₂∕ENO_x) from RCP8.5 (red) and EDGAR4.2 (blue) and the top-down estimate from CHASER (orange) and model-simulated abundance ratios from CHASER (purple) and CAM-Chem (green) chemistry transport models are superimposed. Gray areas are 90 % confidence intervals of the linear fit (black lines). The four Chinese cities represent the four classes/levels of urban development across 12 selected cities in China.

Shown in Fig. 3 are linear trends of annual-mean ΔCO∕ΔNO₂ and ΔSO₂∕ΔNO₂ relative to 2005 values in four Chinese cities. These cities are representative of a certain level of urban development across mainland China. The four levels in this study are defined using broad clustering between the average GDP per capita per year and the rate of change in ΔCO∕ΔNO₂ that is derived from satellite observations. This is shown in Table 2, where a general rule resulting from this analysis would be a classification mainly based on GDP per capita per year, except for Harbin and Wuhan. Combustion-related activities in Shenyang, Beijing, Shanghai, and Shenzhen can be characterized to follow a progression from heavy to light manufacturing, export processing, and service industries (Chan and Yao, 2008). 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₂ ($-\mathrm{5.4}\pm \mathrm{0.7}$ to $+\mathrm{8.3}\pm \mathrm{3.1}$ % yr⁻¹) and decreasing rate of ΔSO₂∕ΔNO₂ reductions ($-\mathrm{6.0}\pm \mathrm{1.0}$ to $-\mathrm{3.4}\pm \mathrm{1.0}$ % yr⁻¹) relative to 2005 (Fig. 3 and Table 2). This pattern in enhancement ratios is not evident in the rate of change in CO, SO₂, and NO₂ column abundance, for which we find increasing rate of decrease in CO ($-\mathrm{0.1}\pm \mathrm{0.3}$ to $-\mathrm{1.0}\pm \mathrm{0.2}\mathit{\%}$ yr⁻¹) and SO₂ ($-\mathrm{1.9}\pm \mathrm{0.9}$ to $-\mathrm{5.5}\pm \mathrm{1.1}\mathit{\%}$ yr⁻¹) abundance, along with a decreasing rate of increase in NO₂ abundance from Shenyang ($+\mathrm{5.2}\pm \mathrm{1.4}$ % yr⁻¹) to Shenzhen (1.8±0.7 % yr⁻¹) (Table 2). This is 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., $\mathrm{\Delta }{\mathrm{CO}}_{\mathrm{comb}}=\langle \mathrm{CO}-{\mathrm{CO}}^{\mathrm{bg}}\rangle$), across these four levels of development.

We have minimized the influence of interannual variations due to meteorology (e.g., changes in air mass) by analyzing molar ratios (e.g., moleCO∕moleNO₂) rather than absolute molar concentrations (e.g., moleCO∕moleair; Parrish et al., 2002, 2006). As the co-emitted species (i.e., CO, SO₂, and NO₂) are subject to the same meteorological conditions (affecting transport, dilution, and lifetime), their enhancement ratios are expected to be less sensitive to meteorology compared to the absolute molar concentrations. This is supported by the fact that decadal ΔCO∕ΔNO₂ as well as ΔSO₂∕ΔNO₂ for different seasons have similar trends (Fig. 2). Previous studies have also proven that the ratios compared to the concentrations themselves are relatively immune to changing meteorological conditions, and can provide insights into the magnitude and temporal trends of the emissions (Parrish et al., 2002, 2006, 2009; Silva et al., 2013; Hassler et al., 2016). In addition, they can be directly compared to the corresponding emission ratios under certain circumstances. However, we note that even though the ratios derived from satellite observations are relatively less sensitive to meteorology, the methodology cannot eliminate all the impacts from meteorology. The enhancement ratios may be impacted by the meteorological conditions because lifetimes of different air pollutants may respond to meteorological conditions differently. Nevertheless, we believe such an impact should not influence our main conclusions for the following two reasons. (1) Our analysis focuses on decadal trends instead of short-term trends. As shown by a previous study, meteorology also plays an important role on relatively short timescales, but meteorology probably plays a lesser role in the longer-term trends (Krotkov et al., 2016). (2) The satellite retrieval samples are taken over the megacities (right above strong emission sources) instead of downwind of the pollution sources, making them more representative of megacity sources.

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 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 a 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) instead of OMI DOMINO. Interestingly, we find that the increasing enhancement ratio of CO to NO₂ in Shenzhen (and to a lesser extent in Shanghai) remarkably resembles the relative changes in CO-to-NO₂ ratios in more developed megacities (Los Angeles and New York) and several urban agglomerations in the United States (see Fig. 3e for Los Angeles and Table 2 and Fig. S1 in the Supplement for all other select cities). More importantly, the increasing pattern that we see in Los Angeles ($\sim +\mathrm{7}\pm \mathrm{1}$ % yr⁻¹) relative to 2005 is generally consistent with the increasing trend ($\sim +$4 % yr⁻¹) after 2007 of a ground-based CO-to-NO_x enhancement ratio in Los Angeles as reported by Hassler et al. (2016). 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). However, the changes in lifestyle concomitant with human development result in a shift to fewer activities (including increased use of renewable energy), 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₂ to NO₂. Similar to previous studies suggesting emissions of CO, SO₂, and NO₂ 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), our finding on this progression in ΔCO∕ΔNO₂ serves as satellite-based evidence of a dominant shift in the cleanness of bulk combustion in a more economically developed city within a developing country like China.

Conversely, there is no clear difference in the observed enhancement ratios (ΔSO₂∕ΔNO₂) and derived enhancements of SO₂ due to combustion (ΔSO_2comb) between cities. The sensitivity of SO₂ to NO₂ relative to 2005 in Shenzhen does not follow the increasing pattern in Los Angeles (Fig. 3b). Unlike ΔCO_comb, ΔSO_2comb in all four Chinese cities still shows a decreasing trend relative to 2005 while ΔSO_2comb in Los Angeles shows an increasing pattern consistent with its ΔCO_comb. There is a striking difference in absolute magnitudes in SO₂ abundance between these cities (as has been reported), reflecting large-scale differences in combustion practice. Yet, the low SO₂ abundance in Los Angeles also makes it difficult to detect possibly large SO₂ point sources (Krotkov et al., 2016). Enhanced SO₂ signal can still be detected as the spatial first-order derivatives of SO₂ with NO₂ at a megacity scale should not be largely (nonlinearly) influenced by its absolute magnitude. We find that there is a tighter correspondence between SO₂ and NO₂ abundance in Chinese cities than in US cities. This might suggest differences in fuel use as SO₂ is mainly produced within a megacity from burning of sulfur-containing fossil fuel (mostly coal, oil, and natural gas) and to a smaller extent from industrial processes (e.g., smelting). Here, we postulate that the absence of an apparent shift in ΔSO₂∕ΔNO₂ across the four Chinese cities is due to continuing heavier reliance of these cities (and China) on coal burning relative to the United States (Wang and Hao, 2012; Bhattacharya et al., 2015; Qi et al., 2016; Yang et al., 2016; Sun et al., 2018; Zheng et al., 2018). In terms of the sectoral share, the majority of NO_x emissions over the Los Angeles basin are from transport according to a recent fuel-based inventory (Hassler et al., 2016), whereas fossil fuel combustion (from power generation and industry) is the most dominant NO_x source in China (Sun et al., 2018). In terms of the energy share, it was estimated that coal accounted for about 69 % and 23 % of the total primary energy consumption in China and the US in 2005, respectively. Actions including usage of low-sulfur coals, installation of flue gas desulfurization (FGD) facilities, and closing of small units have been taken to reduce coal-related emissions in China. The aforementioned procedure to reduce SO₂ in China is most likely to be the dominant driving factor of the declining ΔSO₂∕ΔNO₂ (Li et al., 2018; Zheng et al., 2018). While there are ongoing activities regulating coal-related emissions, coal consumption in China has continued 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). Conversely, the use of coal in the US has been found to be slightly decreasing along with previous adoption of SO₂ control technologies (Taylor et al., 2005). In addition, previous studies have reported recent reduction in NO_x emissions over China since 2011 based on satellite observations and emission inventories (Liu et al., 2016; van der A. et al., 2017). The installation of selective catalytic reduction (SCR) equipment at power plants and new emissions standards for vehicles both contribute to the NO_x emission reduction (Liu et al., 2016; van der A. et al., 2017; Wu et al., 2017). Conversely, based on our analysis of decadal trends (2005–2014), only NO₂ over Shenzhen decreased overall in the decade, while 10-year average changes in NO₂ over Shenyang, Beijing, and Shanghai were overall positive (Table 2). Intra-decadal changes as reported in Liu et al. (2016) (from increasing to decreasing NO_x emissions around 2011) do not contradict the derived 10-year trend in this work, especially over Shenyang and Beijing where NO_x emissions were still rapidly increasing during the first half of the decade (2005–2011). The changes in SO₂ emissions and NO₂ emissions together contribute to the trends of ΔSO₂∕ΔNO₂ that we found. Positive ΔNO₂ and negative ΔSO₂ produce negative ΔSO₂∕ΔNO₂ over the three cities; while negative ΔSO₂ and negative ΔNO₂ (albeit smaller in magnitude) still produce negative ΔSO₂∕ΔNO₂ but of smaller magnitude over Shenzhen than ΔSO₂∕ΔNO₂ over the other cities (Table 2). This indicates a stronger influence of the changes in SO₂ emissions (as reflected in ΔSO₂) in the decreasing trends of these ratios.

3.2 Inconsistencies with a priori estimates

The satellite-based ΔCO∕ΔNO₂ patterns are inconsistent with emission- and model-based ratios (Figs. 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₂, and SO₂) that can be observed from space with MOPITT and OMI. In addition, NO₂ is considered to be the dominant form of NO_x that can be observed at this scale. From a global atmospheric chemistry modeling (chemistry transport models, CTMs) perspective, the associated abundance over megacities is represented as one to four discrete vertical column(s) assuming spatial resolution of these CTMs of 1 to 2^∘. While recognizing the associated month-to-month variability in ΔCO∕ΔNO₂ and expected differences in how these ratios should be compared, the trends in emission ratios relative to 2005 of CO to NO_x from bottom-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 (i.e., decreasing to increasing ΔCO∕ΔNO₂ relative to 2005 from Shenyang to Shenzhen; Fig. 3). This is also true for the ratios of CO to NO₂ 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 on the 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₂ retrievals among other retrievals, also appear to initially follow the emission ratios from EDGAR. Furthermore, the ratios of SO₂ 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 Los Angeles. A closer look at linear trends of the ratios for each sector in RCP8.5 (Fig. S2) reveals inconsistencies in the trends, which cannot be addressed by simple scaling 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 assimilation and inverse methods in integrating observations into models including representativeness of these retrievals. Here we highlight 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 agreement with the trends in ΔCO∕ΔNO₂ (Miyazaki et al., 2017). Such improvements highlight an under-explored utility of available observational 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.

3.3 Combustion emission pathway for Chinese cities

We define the combustion emission pathway as a trajectory in time of the overall changes in emissions due to combustion with respect to socioeconomic development (e.g., Riahi et al., 2011; Steinberger et al., 2012; Li et al., 2016; Marangoni et al., 2017). In this section, we identify a common combustion emission pathway across these four levels of development and associate them with sectoral changes through inverse analysis. We will briefly describe the inverse analysis of the ratios in Sect. 3.3.1, present our findings on combustion emission pathway in Sect. 3.3.2, and elucidate the driving factors by means of time traces in sectoral emission ratios in Sect. 3.3.3.

3.3.1 Inverse analysis of the ratios

We conduct an inverse analysis of the ratios shown in Fig. 3 to further expound on these patterns, by associating them with sectoral changes. Please see details of the matrix–vector product and inversion methodology in Sect. 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₂-to-NO_x ratios, such that the corresponding time series estimates of the total CO-to-NO_x and SO₂-to-NO_x ratios match the decadal trends of ΔCO∕ΔNO₂ and Δ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 or the linear trend in Fig. 3), as these are the most appropriate data for analyzing long-term changes in emission sectors.

Figure 4 Joint traces of the annual changes in a priori (dotted line) and a posteriori (solid line) estimates of ECO∕ENO_x (x axis) and ESO₂∕ENO_x (y axis) relative to 2005 for the four selected Chinese cities (Shenyang: green; Beijing: blue; Shanghai: orange; Shenzhen: purple) representing four levels of urban development. These traces are presented as line arrows (with the origin at x=1, y=1 and endpoint corresponding to the years 2005 and 2014, respectively) for total emission ratios (a) and four sectoral ratios (b–e). The other sector is the sum of mostly residential and commercial areas along with agriculture and waste treatment and disposal. The inset for each panel represents the associated traces for Los Angeles, which is added as basis for comparison. In each panel, the different quadrants (lower left, lower right, and upper right) correspond to decreasing ECO∕ENO_x and ESO₂∕ENO_x, increasing ECO∕ENO_x but decreasing ESO₂∕ENO_x, and increasing ECO∕ENO_x and ESO₂∕ENO_x relative to 2005, respectively. The gray semicircular arrow in panel (a) represents our suggested common combustion emission pathway for Chinese cities.

3.3.2 Combustion emission pathway

The results of our inverse analysis are presented in Fig. 4. This figure consists of five 2-D line plots of a posteriori (solid) and a priori (dashed) time series of SO₂-to-NO_x emission ratios (ESO₂∕ENO_x) on the y axis versus corresponding values of CO-to-NO_x emission ratios (ECO₂∕ENO_x) on the x axis. The five plots correspond to the annual total (Fig. 4a) and sectoral emission ratios (Fig. 4b to e) of each of the four cities selected in Fig. 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. Figure 4a is similar to Fig. 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 Fig. 3. In Shenyang, both ESO₂∕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. Conversely, we see a clear shift in Shanghai and most notably in Shenzhen to a slightly decreasing ESO₂∕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 US) are, however, different than in Shanghai and Shenzhen. In particular, we find ESO₂∕ENO_x and ECO₂∕ENO_x in Los Angeles to be both linearly increasing relative to 2005 values. And so, there is a progression of decreasing to increasing sensitivities of CO and SO₂ to NO₂ from Shenyang to Shenzhen to Los Angeles (gray semicircular trace in Fig. 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). In this case, it 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 the US (Hassler et al., 2016; 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₂, followed by CO, and continues with a reduction in NO_x and potentially in volatile organic compounds (VOCs) later on. To illustrate, we still see increases in NO_x abundance in Shenyang although CO and SO₂ are already decreasing, whereas in Shenzhen, we see NO_x starting to decrease (at a faster rate) along with decreasing CO and SO₂ abundance. The rate at which SO₂, CO, and NO₂ are decreasing is not at a level that is observed in Los Angeles. And so, while the satellite data reveal a combustion emission pathway in these Chinese megacities, these cities have yet to reach conditions that are at par with megacities in more developed cities in the US 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, 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 resembles the emission pathway for Beijing more.

3.3.3 Traces in sectoral emission ratios

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 constitutes more than one-third of NO_x emissions in Chinese (US) cities (Fig. 4b to e). Our inversion results in slight adjustments in the Chinese energy emission pathway towards little to no changes in CO-to-NO_x emission ratios (Fig. 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 (Fig. 4c). This is certainly not the case in Los Angeles where CO-to-NO_x and SO₂-to-NO_x ratios follow quite the opposite pathway of increasing ratios from the energy sector and increasing CO to NO_x, with no change in SO₂ to NO_x from the transportation sector. This is expected in the United States because of cleaner fuel standards (Shindell et al., 2011; Zhang et al., 2012; Kheirbek et al., 2014; Yang et al., 2016; Paulot et al., 2017). 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 (Fig. 4d and e). Shanghai and most notably Shenzhen show a shift to increasing CO-to-NO_x ratio with a slightly decreasing SO₂-to-NO_x ratio that is not 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 mentioned earlier, tertiary (service) industries including export processing activities are more dominant in Shanghai and Shenzhen than in Shenyang. The shift in recent years to an increasing CO-to-NO_x ratio reflects a larger rate of decrease in NO_x levels than CO from the industrial and 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, this is well-supported by studies like Hassler et al. (2016), where they reported an 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 the US recession and implementation of new NO_x emission control technologies and regulations to meet Tier 2 emission standards on US light-duty vehicles. They also noted that differences in the trends of ΔCO∕ΔNO₂ are still observed even between cities from developed countries like the US 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 (e.g., VOCs in petrochemical and pharmaceutical industries) as activity decreases with efficiency, pollution is controlled, and lifestyle changes whenever cities evolve (McDonald et al., 2018). A recent study (Jiang et al., 2018) revealing an overestimation in the decrease in U.S. EPA NO_x emissions based on OMI NO₂ and MOPITT CO retrievals with U.S. EPA ground station measurements of NO₂ 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 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.

Figure 5 Annual-mean enhancement ratios (in units of moles per mole) of CO to NO₂ (a) and SO₂ to NO₂ (b) for all 36 provincial capitals and cities (2005 to 2014) as a function of its corresponding annual GDP per capita (in units of CNY per capita per year). The 12 select cities analyzed in this study are plotted in color, where each color represents four increasing levels or classes of urban development (e.g., Shenyang: Class 1; Beijing: Class 2; Shanghai: Class 3; Shenzhen: Class 4). The rest of the 36 cities are plotted in gray. Superimposed on panels (a) and (b) is a fitted curve (black dashed line) based on the power-law relationship of the data, which is indicated in the plot by its corresponding equation.

Figure 6 Weekly (a) and seasonal (b) cycle of the satellite-based enhancement ratios averaged for the 12 cities in China (red) and for eight cities in the US (blue). The error bars stand for standard deviation across cities.

3.4 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 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 Figs. 3 and 4, our focus is to 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 (EF_species) and effectiveness of control technology (1−CE_species) for CO and NO_x species in the case of ΔCO∕ΔNO₂ for example. We use a robust least-squares regression with a least absolute residuals method to fit a curve of the form y=ax^k, where y is ΔCO∕ΔNO₂ or ΔSO₂∕ΔNO₂ and x is GDP per capita. Our results are presented in Fig. 5a and b 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 the previous section. Note that the magnitudes of enhancement ratios derived from this work are 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 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 k coefficients (R²=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₂ controls should decrease SO₂ abundance. While each city is unique and the evolution of air pollution may be different from city to city, there is also a clear signature of urbanization at the 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, 2013), energy flows (Creutzig et al., 2015), and carbon emissions (Fragkias et al., 2013). Our results suggest that enhancement ratios 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 the highest GDP per capita in China are among the cities with the lowest enhancement ratios. As we have shown in Fig. 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 the highest enhancement ratios. We note, however, that identifying a mechanistic 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 the US that cannot be accounted for (Fig. 6). Nevertheless, we suggest incorporating this observation along with estimates of emissions into future scaling studies, especially as we move past RCPs and toward recent developments in building more realistic emission scenarios that integrate socioeconomic and environmental development pathways like the Shared Socioeconomic Pathways (SSPs; O'Neil et al., 2014).

4 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 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 a megacity. This is especially relevant as the number of megacities continues to grow in the coming decades, mostly at locations that lack sufficient AQ monitoring capabilities. Enhancement ratios of CO to NO₂ and SO₂ to NO₂ 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. This is well-correlated with economic development. These trace a common emission 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 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 and Patra, 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₂ Human Emissions project; https://www.che-project.eu, last access: 22 March 2019). While we recognize the current limitations of these retrievals (e.g., collocation, sensitivity), our findings appear to be robust across retrievals and methods and are 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 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 the 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 is 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., 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 of air pollution (Lewis, 2018) and emission pathways (Mitchell et al., 2018) can be made by (1) considering observationally constrained time-varying emission factors and (2) 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.

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 produced when combustion is incomplete; otherwise carbon in the fuel is oxidized to CO₂ (Eq. 1). In addition, NO and NO₂ are produced from the oxidation of nitrogen from the fuel itself and from decomposition of N₂ in air at high temperatures (Flagan and Seinfeld, 2012). Sulfur dioxide (SO₂) is also produced when the fuel used in the combustion process contains sulfur (such is the case for low-grade fuels).

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{\mathrm{C}}_{x_{\mathrm{1}}}{\mathrm{H}}_{x_{\mathrm{2}}}{\mathrm{O}}_{x_{\mathrm{3}}}{\mathrm{N}}_{x_{\mathrm{4}}}{\mathrm{S}}_{x_{\mathrm{5}}}+n_{\mathrm{1}}\left(\mathrm{1}+e\right)\left({\mathrm{O}}_{\mathrm{2}}+\mathrm{3.76}{\mathrm{N}}_{\mathrm{2}}\right)\\ {\displaystyle }& {\displaystyle }\to n_{\mathrm{2}}{\mathrm{CO}}_{\mathrm{2}}+n_{\mathrm{3}}{\mathrm{H}}_{\mathrm{2}}\mathrm{O}+n_{\mathrm{4}}{\mathrm{O}}_{\mathrm{2}}+n_{\mathrm{5}}{\mathrm{N}}_{\mathrm{2}}+n_{\mathrm{6}}\mathrm{CO}\\ \text{(A1)}& {\displaystyle }& {\displaystyle }+n_{\mathrm{7}}\mathrm{NO}+n_{\mathrm{8}}{\mathrm{NO}}_{\mathrm{2}}+n_{\mathrm{9}}{\mathrm{SO}}_{\mathrm{2}}+n_{\mathrm{10}}\mathrm{C}+\mathrm{\dots }\end{array}$$

Emissions of these intermediate products are typically expressed as

$$\begin{array}{}\text{(A2)}& {\displaystyle }{E}_x=\sum _s\left[A_{\mathrm{s}}\cdot {\mathrm{EF}}_{x,s}\cdot \left(\mathrm{1}-{\mathrm{CE}}_{x,s}\right)\right]=\sum _s\left[A_{\mathrm{s}}\cdot {\mathrm{EEF}}_{x,s}\right],\end{array}$$

where E_x is the total mass of emissions for species x, EF_x,s is its associated emission factor for a specific source/sector s, and A_s is the activity level of the source. CE_x,s corresponds to effectiveness of control measure and ${\mathrm{EEF}}_{x,s}={\mathrm{EF}}_{x,s}\cdot \left(\mathrm{1}-{\mathrm{CE}}_{x,s}\right)$ is the effective emission factor. When we take the ratio of emissions (Eq. 2) of co-emitted species x and y,

$$\begin{array}{}\text{(A3)}& {\displaystyle \frac{E_y}{E_x}}={\displaystyle \frac{\sum _s\left[A_{\mathrm{s}}\cdot {\mathrm{EEF}}_{y,s}\right]}{\sum _s\left[A_{\mathrm{s}}\cdot {\mathrm{EEF}}_{x,s}\right]}}=\sum _s\left({\displaystyle \frac{{\mathrm{EEF}}_{y,s}}{{\mathrm{EEF}}_{x,s}}}\right)\left({\displaystyle \frac{E_{x,s}}{E_{x,\mathrm{total}}}}\right),\end{array}$$

this ratio can be expressed as the sum of the products of the ratio of effective emission factors ($R_{x,y,s}^{\mathrm{EEF}}$) and the fractional contribution of emission sector f for species x (f_x,s) (Eq. A3).

Appendix B: Inverse analysis

We decomposed the a priori estimate of monthly emission ratio of CO to NO_x (and SO₂ to NO_x) from RCP8.5 as a product of (a) the ratio of effective emission factors for each of the four sectors, namely energy, industry, transport, and others ($R_{x,y,s}^{\mathrm{EEF}}$), and (b) fractional contribution of NO₂ emissions from each sector to the total NO₂ emissions (f_x,s) for all four sectors s (s₁: energy; s₂: industry; s₃: transport; s₄: others). In matrix–vector form, this can be expressed as

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\left[\begin{array}{c}\mathrm{ECO}/{\mathrm{ENO}}_x\\ {\mathrm{ESO}}_{\mathrm{2}}/{\mathrm{ENO}}_x\end{array}\right]\\ {\displaystyle }& {\displaystyle }=\left[\begin{array}{cccc}{R}_{\mathrm{CO}/{\mathrm{NO}}_x,s_{\mathrm{1}}}^{\mathrm{EEF}}& R_{\mathrm{CO}/{\mathrm{NO}}_x,s_{\mathrm{2}}}^{\mathrm{EEF}}& R_{\mathrm{CO}/{\mathrm{NO}}_x,s_{\mathrm{3}}}^{\mathrm{EEF}}& R_{\mathrm{CO}/{\mathrm{NO}}_x,s_{\mathrm{4}}}^{\mathrm{EEF}}\\ R_{{\mathrm{SO}}_{\mathrm{2}}/{\mathrm{NO}}_x,s_{\mathrm{1}}}^{\mathrm{EEF}}& R_{{\mathrm{SO}}_{\mathrm{2}}/{\mathrm{NO}}_x,s_{\mathrm{2}}}^{\mathrm{EEF}}& R_{{\mathrm{SO}}_{\mathrm{2}}/{\mathrm{NO}}_x,s_{\mathrm{3}}}^{\mathrm{EEF}}& R_{{\mathrm{SO}}_{\mathrm{2}}/{\mathrm{NO}}_x,s_{\mathrm{4}}}^{\mathrm{EEF}}\end{array}\right]\\ \text{(B1)}& {\displaystyle }& {\displaystyle }\left[\begin{array}{c}{f}_{{\mathrm{NO}}_x,s_{\mathrm{1}}}\\ f_{{\mathrm{NO}}_x,s_{\mathrm{2}}}\\ f_{{\mathrm{NO}}_x,s_{\mathrm{3}}}\\ f_{{\mathrm{NO}}_x,s_{\mathrm{4}}}\end{array}\right]\end{array}$$

or

$$\begin{array}{}\text{(B2)}& {\displaystyle }\mathit{y}=\mathbf{H}\mathit{x}.\end{array}$$

We use a two-step Monte Carlo-based Bayesian inversion method to estimate both H and x of the following cities: Shenyang, Beijing, Shanghai, Shenzhen, and Los Angeles. We focus our analysis on the decadal trends of the RCP8.5 CO-to-NO_x and SO₂-to-NO_x emission ratios using the decadal trends of ΔCO∕ΔNO₂ and ΔSO₂∕ΔNO₂ as observational data (y). We use the decadal trend of enhancement ratios of CO to NO₂ and SO₂ to NO₂ (derived using STL), calculate their annual averages normalized to 2005 values, and then take these as our observational (fitting) data. Our goal is to estimate H and x given y subject to the following constraints: (a) errors in H and x are 10 % and 25 % of their values, (b) errors in y are 5 % of its value, (c) error covariances of y and x are uncorrelated and diagonal (S_e, S_a), and (d) the sum of x is unity. Since this is an under-determined inverse problem, we apply prior information on H and x using the RCP emissions (H_a, x_a). We conduct our inverse analysis in two steps: (1) estimate the most likely H that results in estimates of x best fitting the decadal trend, (2) estimate x using the new estimate of H. For Step 1, first, we draw n=10 000 samples of H assuming its errors are normally distributed, with mean as its prior and covariance as the diagonal of its squared errors. Second, we use the maximum a posteriori (MAP) solution to the Bayesian problem to estimate x for every sample, i.e.,

$$\begin{array}{ll}{\displaystyle }\widehat{\mathit{x}}& {\displaystyle }={\mathit{x}}_a+{\left({\mathbf{H}}_a^T{\mathbf{S}}_e^{-\mathrm{1}}{\mathbf{H}}_a+{\mathbf{S}}_a^{-\mathrm{1}}\right)}^{-\mathrm{1}}{\mathbf{H}}_a^T{\mathbf{S}}_e^{-\mathrm{1}}\left(\mathit{y}-{\mathbf{H}}_a{\mathit{x}}_a\right)\widehat{\mathbf{S}}\\ \text{(B3)}& {\displaystyle }& {\displaystyle }={\left({\mathbf{H}}_a^T{\mathbf{S}}_e^{-\mathrm{1}}{\mathbf{H}}_a+{\mathbf{S}}_a^{-\mathrm{1}}\right)}^{-\mathrm{1}}.\end{array}$$

We draw a new sample if any of the elements in $\widehat{\mathit{x}}$ are negative. Third, we take the mean of 100 H samples resulting in the lowest root-mean-square errors relative to the data. We use this mean as our new estimate of H ($\widehat{\mathbf{H}}$). For Step 2, we apply the same MAP solution using x_a and ${\mathbf{H}}_a=\widehat{\mathbf{H}}$ to estimate $\widehat{\mathit{x}}$ and $\widehat{\mathbf{S}}$. Similar to a Kalman filter, we cycle this procedure for each year starting from 2006 to 2014. We use the new estimates of $\widehat{\mathit{x}}$, $\widehat{\mathbf{H}}$, and $\widehat{\mathbf{S}}$ for a given year as priors for the succeeding cycle with fixed inflation on the covariance of 1.25 to minimize filter divergence. We note that the additional constraints (positive $\widehat{\mathit{x}}$, sum of $\widehat{\mathit{x}}$ is unity) minimize the indeterminacy of the problem. This is supported by post-inverse analysis diagnostics (i.e., averaging kernels) showing that elements of $\widehat{\mathit{x}}$ are resolved by the trend data. Since H is drawn based on Monte Carlo sampling, we do not have a diagnostic for the relative contributions of the prior and the data to H. We chose the mean across 100 H values resulting in estimates of $\mathbf{H}\widehat{\mathit{x}}$ with the lowest RMSEs relative to the data. The changes in $\widehat{\mathbf{H}}$ relative to H_a can be explored in the sectoral changes shown in Fig. 4. This is especially the case for Shanghai and Shenzhen where the change in H is larger than the change in x.