Space-based quantification of per capita CO2 emissions from cities

Urban areas are currently responsible for ∼70% of the global energy-related carbon dioxide (CO2) emissions, and rapid ongoing global urbanization is increasing the number and size of cities. Thus, understanding city-scale CO2 emissions and how they vary between cities with different urban densities is a critical task. While the relationship between CO2 emissions and population density has been explored widely in prior studies, their conclusions were sensitive to inconsistent definitions of urban boundaries and the reliance upon CO2 emission inventories that implicitly assumed population relationships. Here we provide the first independent estimates of direct per capita CO2 emissions (Epc) from spaceborne atmospheric CO2 measurements from the Orbiting Carbon Observatory-2 (OCO-2) for a total 20 cities across multiple continents. The analysis accounts for the influence of meteorology on the satellite observations with an atmospheric model. The resultant upwind source region sampled by the satellite serves as an objective urban extent for aggregating emissions and population densities. Thus, we are able to detect emission ‘hotspots’ on a per capita basis from a few cities, subject to sampling restrictions from OCO-2. Our results suggest that Epc declines as population densities increase, albeit the decrease in Epc is partially limited by the positive correlation between Epc and per capita gross domestic product. As additional CO2-observing satellites are launched in the coming years, our space-based approach to understanding CO2 emissions from cities has significant potential in tracking and evaluating the future trajectory of urban growth and informing the effects of carbon reduction plans.


Introduction
As concerns mount among the global community for stabilizing climate and limiting warming to within 1.5°C from pre-industrial levels (Allen et al 2019), cities have attracted attention as an arena where carbon emissions have the potential for significant reductions. While ongoing urbanization brings challenges to human welfare (Asrar et al 2019), urbanization is also an opportunity to transition cities to be more 'sustainable'. As such, many mayors around the world have made pledges to reduce carbon emissions from their cities (e.g. C40, Watts et al 2015). Against this societal and scientific backdrop, urgent needs exist for studies of carbon emissions and how such emissions vary with socioeconomic drivers. Given the highly dynamic and heterogeneous nature of urban systems, empirical relationships between urban CO 2 emissions and other characteristics about the cities may help inform the effectiveness of carbon reduction actions.
As two of the most widely studied urban characteristics, total population count and population density are closely linked to the intensity of socioeconomic activities and resulting CO 2 emissions (Fragkias et al 2013, Gudipudi et al 2016. Compared to population size, population density contains information about the spatial allocation of individuals, with implications Original content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. for urban planning. Furthermore, per capita emission accounting is commonly used as the basis with which emission reductions obligations as part of climate change agreements (Aldy 2006).
Numerous papers have attempted to elucidate city-scale CO 2 emissions and how they covary with population (and density) with both empirical fits and theoretical explanations (Bettencourt et al 2007, Bettencourt 2013, Fragkias et al 2013, Oliveira et al 2014, Gately et al 2015, Gudipudi et al 2016, Ribeiro et al 2019. For example, a 'sublinear' emission-density relationship implies a slower increase in total emission (i.e. decrease in per capita emission, E pc ) relative to the increase in population density. In addition to population (density), competing factors may affect the efficiency of population density in modifying E pc . To address the impact on emissions from competing factors, Fragkias et al (2013) and Gately et al (2015) examined their emission-population relationship with additional factors such as per capita personal income, job density, and population growth rate. Recent efforts have been made towards overcoming the possible confounding effects from urban area, GDP, and energy consumption, using sophisticated frameworks further inspired by the economic theory of production or the Kaya Identity (Gudipudi et al 2019, Ribeiro et al 2019.
However, prior studies examining per capita CO 2 estimates and how they covary with urban characteristics are subject to two main weaknesses: (1) Definitions of urban extent differ among studies.
Consequently, inconsistent urban definitions lead to discrepancies in per capita urban emissions and divergent scaling relations and conclusions on 'whether larger cities are greener' (Fragkias et al 2013, Oliveira et al 2014. (2) Lack of observational constraints. Most prior studies infer CO 2 emissions from self-reported data or emission inventories. Differing methods and purposes may result in reported emissions being inconsistent between local, regional, and national scales (Guan et al 2012). In addition, population data are often used in inventory development. Thus, analyses of relationships between inventory-based emissions and population could simply reproduce the inventories' underlying assumptions.
As global-scale CO 2 measurements have become available in recent years, studies have emerged to leverage spaceborne CO 2 measurements for studying urban emissions and urban characteristics. For example, column CO 2 (XCO 2 ) enhancements due to urban emissions derived from OCO-2 have been related to urban area sizes and self-reporting fossil fuel emissions over dozens of global cities (Labzovskii et al 2019). However, this study did not explicitly derive CO 2 emissions from the satellite measurements nor account for atmospheric transport.
In this work, we provide the first independent observational constraint on per capita emissions and their relationship with population density and per capita GDP, by using spaceborne CO 2 measurements from OCO-2 , launched by NASA in 2014. OCO-2 retrieves the vertically integrated atmospheric column CO 2 (XCO 2 ) concentrations and has already demonstrated capability for informing the global carbon cycle (Liu et al 2017) and quantifying fine-scale CO 2 emissions from power plants (Nassar et al 2017) and a megacity (Hedelius et al 2018).

Data and methodology
By using OCO-2 measurements (section 2.1) and an atmospheric transport model (section 2.2), our goal is to extract the urban XCO 2 signal and use it to quantify city-wide per capita CO 2 emissions (E pc ) in different cities, with an attempt to examine how the urban emissions scale with variables such as population density and GDP (sections 2.3-2.4).

Target cities and OCO-2 observations
We select 20 cities across the midlatitudes over multiple continents (black boxes in figure 1(a)) based on favorable data quality and sampling patterns from OCO-2, minimal interference from vegetation, and relatively high population density according to the Gridded Population of the World (GPWv4; Center for International Earth Science Information Network CIESIN, Columbia University 2017). As an indicator of potential interference from vegetation, contiguous Solar-Induced Fluorescence (SIF) product (Zhang et al 2018) during north hemispheric non-growing season in 2016 is used to confirm minimal vegetation activities in target cities (figure 1(a)). We also examined the residual biogenic impact on our emission esimtates (SI appendix, Text D2 is available online at stacks.iop.org/ERL/15/035004/mmedia).
Only satellite tracks whose upwind regions overlapping with the urban land during the non-growing season between September 2014 and April 2019 were analyzed for every city. The numbers of qualified satellite tracks span from 6 to 10 for northern hemispheric cities and are smaller for southern hemispheric cities (table S1). In the end, 134 satellite tracks were analyzed for emission estimates. For every track, we selected 50-90 screened satellite soundings with quality flags of zero from version 9 of OCO-2 lite files (OCO-2 Science Team/Michael Gunson, Annmarie Eldering 2018) to interpret track-specific urban XCO 2 signals (further explained in section 2.3). More details on the criteria for selecting cities and satellite tracks are described in SI appendix, text A.

Atmospheric transport model and input datasets
To calculate CO 2 emissions from the observed urban XCO 2 signal and account for wind patterns affecting the satellite measurements, we adopted the atmospheric column version of the Stochastic Time-Inverted Lagrangian Transport (X-STILT, Wu et al 2018) model. 100 air parcels are released per vertical level in the atmospheric columns ('column receptors') corresponding to selected satellite soundings locate. To transport air parcels in the model, we used the Eta Data Assimilation System (EDAS) at 40 km grid spacing for the US cities and the global data assimilation system at 0.5°grid spacing (Rolph et al 2017) for the remaining cities. Errors in the wind fields and their impacts on per capita emissions are carefully investigated (section 2.4).
Next, source-receptor sensitivities, also referred to as 'footprints' (Lin et al 2003, are generated at the horizontal resolution of 1 km×1 km with a spatial coverage of 20°-lat×20°-lon around each city center. 'Column footprints' (X Foot , gray (a) Maps of 1 km×1 km population density (PPS, capita km −2 ; in red) from GPWv4, with daily mean clear-sky contiguous SIF averaged over the north hemispheric non-growing reason in 2016 (mW m −2 nm −1 Sr −1 ; in green). Only areas with high PPS (>500 capita km −2 ) are plotted to better reveal population centers. Spatial domains of 2°×3°(black rectangles) are used to select OCO-2 soundings and to place X-STILT receptors. (b) Maps of latitudinal integrated source-receptor sensitivities (X Foot , ppm-°N/(μmol m −2 s −1 ); in gray) and PPS-weighted X Foot (ppm-°N/(μmol capita −1 s −1 ); in purple) in log 10 scale from X-STILT with XCO 2 measurements from OCO-2 (ppm; blue to red) over Las Vegas on 02/08/2018 2000 UTC. (c) Latitude series of screened XCO 2 (quality flags QF=0, gray triangles) with average values in each latitude bin (black triangles). Background XCO 2 (dashed green line) is defined as the mean observed XCO 2 over modeled background latitude ranges. The latitude range (light red) affected by the urban plume is interpolated from X-STILT. shading in figure 1(b)) are matrices that link potential upwind source regions to downwind atmospheric columns with units of ppm/(μmol m −2 s −1 ), with proper weighting using OCO-2 vertical weighting functions. Therefore, the upwind source regions of the satellite soundings are elucidated by X-STILT, which is also the region for examining the overlap with gridded population density maps (purple shading in figure 1(b)). Only population densities within source regions are taken into account to match with the corresponding urban signal observed by the satellite. Our objective approach of inferring source regions from the atmospheric model avoids having to impose, a priori, city boundaries.
The main dataset used to derive population density is the Gridded Population of the World (GPWv4, Center for International Earth Science Information Network CIESIN, Columbia University 2017). Land-Scan 2017 High Resolution Global Population Dataset (Dobson et al 2000) was used to evaluate possible bias in GPWv4 and its influence on scaling results (SI appendix, Text D1). Both products provide global population data at 1 km grid spacing. We adopted gridded GDP per capita (Kummu et al 2018), to assess how CO 2 emissions covary with socioeconomic variables.
2.3. Calculations of satellite-based urban signal, background XCO 2 , E pc, and effective PPS We solve for an intensive quantity-the per capita CO 2 emissions (E pc )-from all direct emitters over the source region affecting the satellite-observed XCO 2 , as elucidated by X-STILT. Instead of attempting to calculate spatially explicit emissions, we treat the entire city as a whole and estimate an overall emission value per satellite track per city. The city-wide calculation is more robust against errors in the atmospheric model and their impact on emission estimates (Wu et al 2018).
E pc per track per city is calculated as the ratio between an urban signal and a population density (PPS)-weighted column footprint as showed in equation (1): where x, y are indices for the spatial grid, l is the latitude of every receptor/sounding. Spatial column footprints X Foot (x, y, l) are generated at every 50-90 column receptors for a given track. The integral in the numerator is a latitudinally integrated XCO 2.ff enhancement due to fossil fuel combustion, or simply an urban XCO 2 signal. The entire summation in the denominator describes the overlapping region between X Foot and PPS, which is revealed by the purple shading in figure 1(b). More specifically, XCO 2.ff along the satellite swath are calculated as differences between a track-specific background and observed total XCO 2 . We followed the forward-simulation approach presented in Wu et al (2018) to estimate an urban-enhanced latitude range (from l 1 to l 2 ; red ribbon in figure 1(c)) and background value over background region outside the latitude range. As not all source areas indicated by X Foot are associated with high population density, gridded PPS from GPWv4 is convolved with latitudinally integrated X Foot . Lastly, a linear regression curve is fitted whose slope serves as the one overall E pc estimate per city (figure S1).
To examine the relationship between E pc and PPS, a representative PPS value is required for every city. Instead of simply taking the average of PPS within a prescribed urban area, we calculate an effective population density (PPS eff ) as the footprint-normalized PPS: where XF tot is the sum of footprints excluding those over lands or water bodies whose PPS are zero and X Foot serves as a spatial weighting factor. In other words, population densities within the urban core are often weighted more, with more leverage on influencing the satellite-observed XCO 2 . Therefore, PPS eff vary among different satellite tracks under different meteorological conditions. Besides PPS eff , we also calculated average PPS over the source region (PPS upwind ) to scale against E pc as a sensitivity test (SI appendix, text D1). Similar to the calculation of PPS eff , effective per capita GDP is also calculated for every track and city with spatial footprint serving as an areal weighting function.

Uncertainty estimates for E pc
The overall fractional error of E pc is comprised of fractional errors in both observed urban signals and the atmospheric model (namely, from horizontal transport), according to E pc calculation in equation (1). Vertical transport in X-STILT typically has a smaller uncertainty impact than horizontal transport (Wu et al 2018) and thus excluded from this study. Both observational and modeling errors are first quantified for each column receptor and then aggregated to track-level errors, after accounting for error covariances. Finally, by assuming statistical independence of errors between tracks, track-level errors decay as a function of N , where N represents the numbers of tracks per city. Technical details of error quantifications are available in SI appendix, text C2.

OCO-2/X-STILT E pc estimates
The average urban signal observed by OCO-2 varies from 0.21 to 1.56 ppm (Lat-avg. column in table S1). Our estimates of E pc range from 1.22 to 44.30 tCO 2 yr −1 . Estimates especially for low-emission cities fall within the same approximate range as national-level estimates from World Bank (World Bank 2014) or the Emissions Database for Global Atmospheric Research (EDGAR, Janssens-Maenhout et al 2019), even given differences in reporting years. Relative uncertainties of E pc fluctuate from 16.0% to 46.9% for different cities (table S1). These uncertainties stem from a combination of measurement noise and errors in retrieving satellite XCO 2 , background XCO 2 , and the winds used to drive X-STILT.
A few cities stand out especially due to their high E pc . Yinchuan is one of the four major coal production cities in China (Liu and Cai 2018). Due to rich coal resources and intensive coal production (Liu and Cai 2018), E pc for Yinchuan is about five times higher than China's national-level estimates (table S1). In contrast, two other cities in China, Lanzhou and Xi'an, have comparable or slightly higher E pc as the national average level. Middle Eastern cities are well-known for their natural gas and oil productions, while high E pc of Johannesburg may be attributed to its coal-fired power plants as also seen from NO x instruments like TROPOMI (Reuter et al 2019).
To further understand the unique emission structures for individual cities and plausible causes for the extremely high E pc , we adopt sectoral breakdown of emissions from an emission inventory (EDGARv5.0, Crippa et al 2019, https://data.europa.eu/doi/10. 2904/JRC_DATASET_EDGAR). These sectoral emissions, independent from OCO-2 based emissions, are utilized to calculate the XCO 2 percentage from every sector over all fossil fuel CO 2 (FFCO 2 sectors (ppm/ ppm; in %). It is noteworthy that ∼12.9%-81.5% of the EDGAR-modeled XCO 2 enhancements are caused by power industries for different cities, and cities with high E pc have high shares from their power industries ( figure S2(a)). This finding underscores the importance of power plants in overall urban emissions (Singer et al 2014). As direct emissions from power industries are often unrelated to urban population (density) while emissions from buildings are expected to be related, we adopted the power-to-building XCO 2 enhancement ratio ( figure S2(b)) as an indicator of the relative importance of power industries. This indicator identified Riyadh, Yinchuan, and Johannesburg as cities with heavy power industries given their high power-to-building concentration ratios.

Possible factors that drive E pc
We next turn to the scaling relationship between E pc and the effective population density (PPS eff , equation (2))an overall population density representing each city. E pc generally declines as PPS eff rises, even after uncertainties in E pc are accounted for via Monte Carlo simulation. 5000 slopes of the E pc -PPS eff relation in logarithmic space from the Monte Carlo analysis (blue lines in figure 2) have a mean value of -1.06 with a standard deviation of 0.083. Since cities with heavy power industries (defined in section 3.1) may skew the E pc -PPS eff relation, we also report regression slopes with mean of -0.97 (gray lines in figure 2) after removing such cities from consideration.
In addition to PPS, we examine the relationship between E pc and per capita GDP (GDP pc ), which may shed light on E pc 's dependency on city-wide living standards. E pc increases linearly with GDP pc in logarithmic space, excluding a few industrial cities (figure S3). As a simple attempt to disentangle emission effects from GDP pc , we carried out multiple regression fits using both PPS and GDP pc in logarithmic space with E pc 's uncertainties accounted via Monte Carlo simulations. For all 20 cities, coefficients of GDP pc are all positive with a mean of 0.59, while coefficients of PPS become less negative than the previous fit using just PPS as the explanatory variable (-1.06 to -0.83), implying that the decrease in E pc is partially limited by the positive correlation between E pc and GDP pc .
Possibly to a lesser extent, E pc -PPS relationship could be affected by other factors, e.g. emission seasonality caused by lack of summertime tracks, potential biases in GPWv4 and PPS eff calculations, and minor residual vegetation influences gradient. For instance, E pc using LandScan may differ from those using GPWv4, especially for Doha and Riyadh (5th versus 6th columns in table S1). However, sensitivity analyses suggest that the inverse E pc -PPS relation is robust (SI appendix, text D). Lastly, population dynamics (e.g. transient versus resident population) may potentially affect the track-level E pc and scaling slopes, but has not been studied in this paper due to unavailability of such datasets to our knowledge.

Production-and consumption-based emission accounting
These are ultimately two different accounting methods for anthropogenic carbon emissions, also reflected in emission inventories that are either consumptionbased Our estimates represent direct emissions from local emitters from a production-based perspective, in contrast to the concept of 'carbon footprint' or 'Scope 3 emission' (Moran et al 2018) which is calculated based on consumption (e.g. material, fuel, travel). Therefore, comparisons between estimates of direct emission and carbon footprint help identify various cities' roles as potential net carbon exporters/ importers (Moran et al 2018). For example, Johannesburg, Riyadh, Urumqi, and Yinchuan have much larger direct emissions than carbon footprints, whereas Madrid and Seoul have smaller direct estimates (table S1), which likely implies the former three cities being net 'carbon exporters' and the latter two cities as 'importers'.
The scaling relationship obtained in this study is estimated based on a 'production-based/direct' emissions perspective, as informed from the satelliteobserved CO 2 plumes that attribute emissions to sources at locations where the CO 2 is emitted to the atmosphere. However, urban residents rely on food and energy produced elsewhere. It is also possible that the scaling relationship could differ for the alternative consumption-based perspective, since residents of large and dense cities can consume goods and services that result in emissions elsewhere Kammen 2014, Moran et al 2018). For example, it has been suggested that the carbon footprint per household (consumption perspective) and its variation with population density differs between urban cores, suburbs, and rural areas (Jones and Kammen 2014). However, we did not examine the population density dependence on CO 2 emissions at sub-urban scale or for subsets of population densities, given limited sample size and increasing complexity for calculating E pc over suburbs and rural places from lower signal-to-noise ratio in XCO 2 and potential interference from biospheric fluxes.
While consumption-based emission accounting brings unique insights to international mitigation efforts, it may suffer from higher uncertainties due to emission reallocations (Peters 2008). Nonetheless, we stress that both the direct emissions perspective adopted in our analysis and the consumption-based emissions perspective are equally important and relevant, from both research and mitigation standpoints.

Discussion
Studies focusing on urban agglomeration and carbon emissions have increased in number within the past few decades. Some have related total carbon emissions to total population (Fragkias et al 2013, Oliveira et al 2014, Gudipudi et al 2019 while others used per capita emissions versus population density (Gately et al 2015, Gudipudi et al 2016. The conversion between these two scaling relations requires an assumption of urban area when converting population density to total population or E pc to total emissions (SI appendix, text B). The purpose and advantage of our top-down approach is to avoid having to explicitly outline the spatial extent of the urban area for collecting emissions or socioeconomic data by relying on the atmosphere to sample the emissions upwind of satellite tracks ( figure 1(b)) and focusing on an intensive quantity like E pc . Therefore, we discuss our results in light of Figure 2. Urban scaling relations between OCO-2/X-STILT derived per capita emissions (E pc , tCO 2 capita −1 yr −1 ) and mean effective population density (PPS eff , capita km −2 ) for 20 cities in log space. Three cities with heavy power industries are shown with blue points, while the other cities are in black. Error bars represent the total error impacts from observed signals and horizontal transport within X-STILT. Linear regression slopes and y-intercepts are calculated using ordinary least squares (OLS) regression method, based on the best-estimates and standard deviations of E pc via 5000 Monte Carlo simulations. Regression lines with the mean values (-1.06 versus -0.97) of 5000 slopes and y-intercepts for all 20 or 17 less industrial cities are plotted in light blue or black lines.
previously-published papers reporting the relationship between E pc and population density. A few decades ago, an inverse relation was first discovered between per capita annual gasoline use and population density for dozens of large cities around the globe (Newman and Kenworthy 1989). A similar decline in per capita on-road CO 2 emission with increased population density was found for US cities (Gately et al 2015). Outside the transportation sector, residential and commercial CO 2 emissions also scale sublinearly with population density over US cities (Gudipudi et al 2016).
More recently, an analysis using ground-based CO 2 measurements derived a nonlinear relation between the total daytime fossil fuel CO 2 (FFCO 2 ) fluxes and PPS eff across five local regions in Salt Lake City . Total FFCO 2 fluxes first increase with PPS eff but then stabilize for PPS eff >200 individual km −2 , suggesting denser neighborhoods even within the same urban area are associated with lower per capita emissions. Linear, sublinear and superlinear relations are associated respectively with individual human needs; infrastructure; and creation of information, wealth and resources (Bettencourt et al 2007). Thus, we speculate that our derived sublinear relation can likely be explained by denser cities' capabilities in enhancing productivities and reducing per capita emissions for sectors like transportation, where infrastructures and resources are shared and commuting distance and cost are largely reduced (Bettencourt et al 2007, Glaeser and Kahn 2010). The use of renewable or nuclear energy could affect E pc by replacing traditional energy sources like coal or natural gas. However, the quantitative impact on E pc and the scaling relation from switching to clean energy is difficult to estimate with the current dataset.
We acknowledge that the city and track numbers sampled in this study, while covering multiple continents, remain small. The limited sample size renders difficult detailed comparisons within or across nations. Therefore, we only attempted to focus on the overall scaling relation without further breakdowns. Although we tried to account for competing factors' impact on E pc -PPS scaling relationship, more comprehensive studies using more observational data are needed. Nevertheless, the limited number of cities analyzed here arise from the small number of suitable tracks, the choice to filter out seasons and areas with potential interference from vegetation, and the fact that OCO-2 was not explicitly designed for urban monitoring but for large-scale, global observations of atmospheric CO 2 .

Conclusion
Despite low sample sizes, our study offers a unique viewpoint at the scale of urban carbon emissions using independent space-based observations, rather than the erstwhile approach of relying on emission inventories. In the upcoming few years, rapid growth in satellite observations is expected. For instance, the 'snapshot area mapping' (SAM) mode from OCO-3 (Eldering et al 2019) would yield more frequent observations with a wider spatial coverage around hundreds of cities. OCO-3's SAM mode will carry out additional XCO 2 observations distributed densely over a region spanning approximately 80 km by 80 km, centered over the urban core of cities from around the world. Thus, researchers will have the opportunity to derive CO 2 emissions from different parts of an urban area and relate them with population statistics or other socio-economic indicators at sub-urban scales. The investigation into the temporal trends of carbon emissions, population densities, or their resultant scaling relations would also be meaningful due to ongoing carbon mitigation efforts.