Southern California megacity CO₂, CH₄, and CO flux estimates using ground- and space-based remote sensing and a Lagrangian model

We estimate the overall CO₂, CH₄, and CO flux from the South Coast Air Basin using an inversion that couples Total Carbon Column Observing Network (TCCON) and Orbiting Carbon Observatory-2 (OCO-2) observations, with the Hybrid Single Particle Lagrangian Integrated Trajectory (HYSPLIT) model and the Open-source Data Inventory for Anthropogenic CO₂ (ODIAC). Using TCCON data we estimate the direct net CO₂ flux from the SoCAB to be 104 ± 26 Tg CO₂ yr⁻¹ for the study period of July 2013–August 2016. We obtain a slightly higher estimate of 120 ± 30 Tg CO₂ yr⁻¹ using OCO-2 data. These CO₂ emission estimates are on the low end of previous work. Our net CH₄ (360 ± 90 Gg CH₄ yr⁻¹) flux estimate is in agreement with central values from previous top-down studies going back to 2010 (342–440 Gg CH₄ yr⁻¹). CO emissions are estimated at 487 ± 122 Gg CO yr⁻¹, much lower than previous top-down estimates (1440 Gg CO yr⁻¹). Given the decreasing emissions of CO, this finding is not unexpected. We perform sensitivity tests to estimate how much errors in the prior, errors in the covariance, different inversion schemes, or a coarser dynamical model influence the emission estimates. Overall, the uncertainty is estimated to be 25 %, with the largest contribution from the dynamical model. Lessons learned here may help in future inversions of satellite data over urban areas.

1 Introduction

About 43 % of global anthropogenic carbon dioxide (CO₂) emissions come directly from urban areas, and urban final energy use accounts for about 76 % of CO₂ emissions (Seto and Dhakal, 2014). Associations of cities that recognize their significant emissions of CO₂ to the atmosphere – such as the C40 Cities Climate Leadership Group (C40) – seek to reduce their greenhouse gas (GHG) emissions and develop local resilience to changing climate. There is a need to track long-term anthropogenic GHG emissions from urban areas to aid urban planners and ensure commitments are met.

Bottom-up (BU) inventories (e.g., of CO₂) can be derived by accounting for various emission activities such as transportation, electricity generation, industry, and heating. BU inventories have some inherent uncertainty due to imperfect emission models, which are largely based on extrapolation of controlled studies and rely on assumptions of fuel consumption, and due to disagreements in downscaling methods (Duren and Miller, 2012; Sargent et al., 2018). Uncertainties in how emissions are calculated and in the underlying activity data used to construct inventories make them susceptible to systematic biases by nature (Oda et al., 2017). On the national level, 2σ uncertainties range from 4.0 % to 17.5 % for the 10 largest emitters (Oda et al., 2018). Uncertainties on the grid cell level are unique to the disaggregation method, but may be in the range of 4 %–190 % (2σ) (Andres et al., 2016). Top-down (TD) emission estimate methods rely on measurements of gases along with models of atmospheric transport, which have their own inherent uncertainties. Measures of emissions and emission changes are generally more reliable when TD and BU methods are in agreement (Duren and Miller, 2012).

Tracking emissions from a TD perspective requires observations. Various networks, such as the Total Carbon Column Observing Network (TCCON) and the National Oceanic and Atmospheric Administration (NOAA) Earth System Research Laboratory (ESRL) in situ CO₂ network can aid in long-term measurements but are too sparse to track emissions from more than a few cities. Some urban areas have ground-based networks (Lauvaux et al., 2016; Mitchell et al., 2018; Sargent et al., 2018; Shusterman et al., 2016; Verhulst et al., 2017). Significant progress has been made in minimizing the cost, deployment time, and data delivery from these networks. However, they still require a significant number of personnel hours and are difficult to scale up to more than a few dozen areas for long-term observations. Urban observation networks can provide finer spatial and temporal details on emission sources, but space-based observations are likely the only way to track emissions TD for more than a few dozen cities.

Within the past 10 years, two satellites have been shown to have high-precision (better than 1 ppm) small-footprint (<100 km²) CO₂ observing capabilities, including the Greenhouse Gases Observing Satellite (GOSAT, in orbit 2009) and the Orbiting Carbon Observatory-2 (OCO-2, in orbit 2014). Several other satellites are planned or are already in orbit with this same potential. Combined, OCO-2 and GOSAT can cover about 1 % of the Earth's surface every 3 days, and though this is only a small fraction, it is unprecedented. Other missions such as TanSat (in orbit 2016), GAS onboard FY-3D (in orbit 2017), GOSAT-2 (in orbit 2018), OCO-3 (expected 2019), and GeoCARB (expected 2023) may further bolster coverage. Space-based observations of methane (CH₄) have been made from GOSAT and the TROPOspheric Monitoring Instrument (TROPOMI, in orbit 2017) and will be made from the GOSAT-2 and planned GeoCARB missions. Carbon monoxide (CO) is measured using Measurements of Pollution in the Troposphere (MOPITT, in orbit 1999) and TROPOMI and will be from GOSAT-2. There is presently a lack of studies that have assimilated satellite trace-gas abundance data into inversion schemes to determine urban emissions.

We test trajectory-based inversion schemes to see if they can reproduce known emissions (from inventories and previous studies) from the California South Coast Air Basin (SoCAB). Our goal is not to apportion spatially, but rather to come up with a single number for the total flux and an estimate of uncertainty. Fluxes from this urban area (pop. ∼ 16.3 million) have been studied extensively, and it provides a test bed to evaluate methods. We discuss the components used to build our inversion in Sect. 2. Typical urban enhancements are described in Sect. 3. Fluxes of CO₂, CO, and CH₄ using TCCON data and of CO₂ using OCO-2 data are discussed in Sect. 4 along with sources of uncertainty. In Sect. 5 we discuss emission ratios, which can also be used to evaluate our flux results. We conclude by summarizing uncertainty and mentioning expansions and areas of improvement in Sect. 6.

2 Data sources and methods

2.1 Observations of column-averaged dry-air mole fractions

We use observations of column-averaged dry-air mole fraction (denoted X_gas) to tie model abundances to fluxes. Column averages are calculated by dividing the retrieved amount of the gas of interest (molecules cm⁻²) by the retrieved total column of dry air (molecules cm⁻²). X_gas values are less sensitive to changes in surface pressure and water vapor than total column amounts in units of molecules per square centimeter (Wunch et al., 2015).

Data are obtained from the TCCON and OCO-2. We use TCCON data from the California Institute of Technology (Caltech) site in Pasadena, California (Wennberg et al., 2014), and from the NASA Armstrong Flight Research Center (AFRC) site near Lancaster, California (Iraci et al., 2014). Values of $X_{{\mathrm{CO}}_{\mathrm{2}}}$, X_CO, and $X_{{\mathrm{CH}}_{\mathrm{4}}}$ were generated using the operational GGG2014 algorithm (Wunch et al., 2015). The Caltech site (lat 34.136, long −118.127, 240 m a.s.l.) is located in an urban environment within the SoCAB. As the name implies, the SoCAB is a basin surrounded by mountains, except towards the southwest, which boarders the Pacific Ocean. AFRC (lat 34.960, long −117.881, 700 m a.s.l.) is located outside the basin ∼100 km to the north in a much more sparsely populated area. Because of the lower population density, the AFRC is often considered a “background” site. However, depending on airflow patterns, recent emissions from the SoCAB may be observed at the AFRC so we use the term background loosely to indicate where lower concentrations are typically observed. Coincident data from both sites are available from July 2013 to August 2016, after which the AFRC instrument was relocated. In total, there are 5355 paired hourly averaged observations on 783 days.

OCO-2 data are available starting September 2014 when the instrument began its nominal operational mission (OCO-2 Science Team et al., 2017). Here, we use $X_{{\mathrm{CO}}_{\mathrm{2}}}$ data generated using the NASA Atmospheric CO₂ Observations from Space (ACOS) version 8r algorithm (O'Dell et al., 2018). We also do a partial analysis on v7r data for comparison with past studies that used these data with a focus on the SoCAB (Hedelius et al., 2017a; Schwandner et al., 2017). Because OCO-2 is in a sun-synchronous orbit with an equatorial crossing time of around 13:00 local solar time, all observations are in the early afternoon. OCO-2 has 8 longitudinal pixels, with a footprint of ∼3 km² each. To reduce over-weighting target mode observations, OCO-2 data are gridded to $\mathrm{0.01}{}^{\circ }\times \mathrm{0.01}{}^{\circ }$. Before filtering there are 6098 pre-averaged OCO-2 observations on 29 different overpass days when the AFRC TCCON site also collected background observations.

In Appendix A we describe filtering, background subtraction, boundary conditions, and our accounting for averaging kernels. In short, we determine enhancements of various gases (ΔX_gas) by finding the difference between observations within the basin (either the Caltech TCCON or OCO-2) compared with the AFRC TCCON site.

2.2 A priori flux estimates

Our flux estimate involves scaling the a priori spatial inventory, or subregions of the prior up or down to reduce the measurement–model mismatch. More important than the total prior absolute flux is the distribution of sources. EDGAR (Emissions Database for Global Atmospheric Research; EC-JRC/PBL, 2009) and FFDAS v2.0 (Fossil Fuel Data Assimilation System; Asefi-Najafabady et al., 2014) are available globally at a 0.1^∘ resolution. We use the year 2016 version of the Open-source Data Inventory for Anthropogenic CO₂ (ODIAC2016), which is available globally at a resolution of 30 arcsec from 2000 to 2015 (Oda and Maksyutov, 2011, 2015; Oda et al., 2018). We also compare total SoCAB emissions from the 2015 version of ODIAC (ODIAC2015), which is based on a projection of the Carbon Dioxide Information Analysis Center (CDIAC) country total emissions. ODIAC has a monthly variation and compared to the annual average seasonal flux ratios are 1.06 (DJF), 0.97 (MAM), 1.00 (JJA), and 0.97 (SON). We assume that 2015 emissions are identical to those in 2016. A generic temporal hourly scaling factor product (TIMES – Temporal Improvements for Modeling Emissions by Scaling) available at a $\mathrm{0.25}{}^{\circ }\times \mathrm{0.25}{}^{\circ }$ resolution can be applied to spatial inventories such as ODIAC to improve temporal emissions (Nassar et al., 2013). However, TIMES has a single peak for midday emissions, which is inconsistent with morning and afternoon rush hour periods in the SoCAB. We instead use the Hestia-LA v1.0 weekly profile reported by Hedelius et al. (2017a, Fig. 2 therein), which has both morning and afternoon rush hour peaks. We use ODIAC over the domain longitude −121.5 to −114.5 and latitude 30.5 to 37.5. Hestia-LA v2.5 is expected to be an even more accurate spatiotemporal inventory for the SoCAB (Gurney, 2018; Gurney et al., 2012). As a sensitivity test we also derive a flux based on Hestia-LA 2.5 over the region in which it is available, and Vulcan 3.0 is used for the rest of the area within the US. These were gridded to the same scale as the ODIAC.

This same prior is used for CO, but total emissions are 1 % of CO₂ emissions on a molar basis (0.6 % of mass) based on the results of Wunch et al. (2009). Figure 1 shows the ODIAC2016 prior for 1 month.

Figure 1 A priori flux maps for CO₂ (a) and CH₄ (b) for select months. The same spatiotemporal prior for CO₂ (ODIAC2016) was used for CO, but scaled to 1 % on a per mole basis. The methane prior was created based on point sources, total emissions, and the population distribution. The black lines are coastlines and the geopolitical boundaries of the SoCAB. Blue lines are county borders.

A detailed CH₄ inventory is also available for the SoCAB, which we do not use because it would be difficult to scale (Carranza et al., 2018). For the US the Harvard–EPA (Environmental Protection Agency) inventory is already available at $\mathrm{0.1}{}^{\circ }\times \mathrm{0.1}{}^{\circ }$ resolution (Maasakkers et al., 2016), and globally the EDGAR inventory is available at $\mathrm{0.1}{}^{\circ }\times \mathrm{0.1}{}^{\circ }$ resolution (EC-JRC/PBL, 2009). We make our own 30 arcsec × 30 arcsec methane prior using landfills, nightlights, expected total emissions, and the US Harvard–EPA inventory (Maasakkers et al., 2016) shown in Fig. 1. Due to a lack of information outside the US on point sources, such as landfills, our methane prior is also not scalable beyond a national level.

For our methane prior we first distribute emissions from landfills as point sources (available 2010–2015, https://ghgdata.epa.gov/ghgp/main.do, last access: 17 August 2017) and use 2015 emissions for 2016. Emissions from the Puente Hills landfill were doubled because the EPA estimate (average of 13.6 Gg CH₄ yr⁻¹) is low compared to previous estimates of 34 Gg CH₄ yr⁻¹ (Peischl et al., 2013). After doubling Puente Hills emissions, EPA total SoCAB (144 Gg CH₄ yr⁻¹) and Olinda Alpha (13.5 Gg CH₄ yr⁻¹) landfill emissions are similar enough to other studies (Peischl et al., 2013) that we do not double emissions from other landfills in the SoCAB. Chino dairy emissions were added in as a $\sim \mathrm{0.1}{}^{\circ }\times \mathrm{0.1}{}^{\circ }$ source (Chen et al., 2016; Viatte et al., 2017). Outside of the SoCAB, CH₄ manure and enteric fermentation were added from the $\mathrm{0.1}{}^{\circ }\times \mathrm{0.1}{}^{\circ }$ Harvard–EPA inventory (Maasakkers et al., 2016). SoCAB emissions are assumed to sum to 400 Gg CH₄ yr⁻¹ based on the work of Wunch et al. (2016), and the rest of the emissions were distributed based on population, which was assumed to correspond with the January 2017 Suomi NPP nightlights (15 arcsec). An average monthly trend was included based on results of Wong et al. (2016), and emissions were assumed to be constant on a monthly timescale. Because the Aliso Canyon leak effectively doubled the SoCAB CH₄ emissions for its duration from 23 October 2015 to 11 February 2016 (Conley et al., 2016), it was also added as a point source.

We use various publicly available statistics to get a sense of annual CO₂ emissions from the SoCAB. Literature estimates range from 99 Tg CO₂ yr⁻¹ (Fischer et al., 2017) to 211 Tg CO₂ yr⁻¹ (Wunch et al., 2009). Table 1 lists statistics for the SoCAB. We assume the nonresidential natural gas (NG) use is for industry or power already accounted for in the EPA inventory. Because most of the food consumed in the SoCAB is grown outside the basin, such as in the Midwestern US and Central Valley (CV), there is a CO₂ return flux to the croplands from both human respiration and food waste. In the US, 60 million metric tons (MMT) of food are lost annually at the retail and consumer levels compared with 129 MMT consumed (Dou et al., 2016), roughly one-third of all food calories (not counting inedible food-related biomass). Presumably, most food waste decomposition would be accounted for in EPA landfill emissions. However, CO₂ emissions from food waste could be underestimated if food waste is composted, if there were unaccounted for methanotrophs, or if aerobic respiration is significantly underestimated (e.g., from rapid decomposition while still exposed to oxygen), which would decrease the CH₄ : CO₂ emission ratio commonly assumed to be unity for managed landfills on a per mole basis (RTI, 2010). Thus, we add 30 % to human respiration emissions of 917 g CO₂ d⁻¹ person⁻¹ (Prairie and Duarte, 2007) for food waste losses. We assume the flux from vegetation is balanced (i.e., no net change in plant biomass or soil carbon) within the basin. This choice is because of uncertainty as to whether there is a net uptake of CO₂ by the biosphere in the SoCAB (Park et al., 2018) or if the excess CO₂ in the atmosphere from the biosphere (Newman et al., 2016) is due to more respiration than photosynthetic uptake. We estimate the uncertainty due to the biosphere is less than ±10 %. Based on these various statistics we estimate a BU net flux of the order of 110 Tg CO₂ yr⁻¹ from the SoCAB.

Table 1 Statistics for the SoCAB.

Most of these values are approximations. ^a https://www.aqmd.gov/nav/about/jurisdiction, last access: 12 November 2018. ^b http://www.dot.ca.gov/hq/tsip/hpms/datalibrary.php, last access: 12 November 2018. ^c Assuming 95 % of kilometers light duty vehicles with 9.1 km per liter (KPL) fuel efficiency and 5 % trucks with 2.5 KPL (https://www.fhwa.dot.gov/policyinformation/statistics/2013/, last access: 12 November 2018, VM-1). ^d Vehicle kilometers and fuel emissions are independent estimates. ^e http://www.cdtfa.ca.gov/taxes-and-fees/spftrpts.htm, last access: 12 November 2018. ^f Based on emissions of 1.3×917 g CO₂ d⁻¹ person⁻¹ (Prairie and Duarte, 2007). ^g http://www.ecdms.energy.ca.gov/gasbycounty.aspx, last access: 12 November 2018. ^h https://www.epa.gov/sites/production/files/2015-07/documents/emission-factors_2014.pdf, last access: 12 November 2018. ⁱ Emissions within or near geographical SoCAB boundaries only.

2.3 Dynamical models

A dynamical model is needed in conjunction with the a priori flux estimates to generate forward model X_gas enhancements. Our model uses Lagrangian trajectories driven by existing archived forecast or reanalysis datasets. An advantage of archived model data is there is no need to run a Eulerian model first, and they are more accessible to a broader community. However, taking existing results without model evaluation may propagate hidden errors and biases, which could influence flux results. Archived data usually have coarser spatiotemporal resolutions than custom models and cover larger domains than the area of interest. Custom runs allow models to be parameterized differently and nudged to reduce the measurement–model mismatch for the regions of interest.

We use the North American Mesoscale Forecast System (NAM) at 12 km resolution (3 h temporal) from the NOAA data archive as the primary model source. NAM is run with a non-hydrostatic version of the WRF at its core with a Mellor–Yamada–Janjić planetary boundary layer (PBL) scheme (Coniglio et al., 2013). Estimates of model error are described in Appendix B. Though NAM data are only available over North America, other archived models are available at lower resolution with global coverage (e.g., the Global Data Assimilation System (GDAS) 0.5 ^∘, 3 h product). The NOAA ESRL recently began publicly releasing 3 km, 1 h archived data from the High-Resolution Rapid Refresh (HRRR) model that covers the US (Benjamin et al., 2016). This product holds the potential to improve flux estimates at smaller scales.

We use HYSPLIT-4 (Stein et al., 2015) with the three archived NOAA data products described above. Our base method is to use mean 48 h back trajectories with NAM 12 km for the lowest 20 % of the atmosphere, which we assume is the only part of the atmosphere enhanced with local emissions at the measurement site. Trajectories are equally spaced in pressure every 0.3 % of the column. By comparison, the GDAS model takes 0.71±0.18 (1σ) times as long to run, and the HRRR model takes 33.3±7.1 (1σ) times as long. Because HRRR takes substantially longer, we only run it for a subset of months – July and October 2015 and January and April 2016. Other studies (Fischer et al., 2017; Janardanan et al., 2016) used multiple particles released at each level. We assume that over the multiyear time series the ensemble of mean trajectories is, on average, representative of the upwind influences on the receptor sites without the additional turbulence term.

Figure 2 shows back trajectories for one layer that end at the observation sites at 14:00 (UTC−7). Trajectories from multiple vertical levels are combined to determine residence times or footprints as described in Appendix C. There are three major origins for air at the Caltech site. The primary source is from over the ocean and over downtown Los Angeles (southwest). The second major source is from the Mojave Desert (northeast), and the third source is from the CV (northwest; see Fig. E1 in Appendix).

Figure 2 HYSPLIT 400 m a.g.l. back trajectories for NAM 12 km for 2015. For each day trajectories are shown ending at the two different TCCON receptor sites at 14:00 (UTC − 7). Magenta trajectories end at Caltech. Cyan trajectories end at AFRC. The black lines are coastlines and the geopolitical boundaries of the SoCAB. Blue lines are county borders.

Others interested in undertaking similar studies may also consider using the recently developed X-STILT (X-Stochastic Time-Inverted Lagrangian Transport) model to obtain footprints for column observations (Wu et al., 2018).

2.4 Inverse methods for comparing measured to model data

Different schemes can be applied to reduce the measured–model mismatch. One of the simplest is to find the ratio between the average enhancements in the observations and the forward model and then to scale the prior based on this ratio. Bayesian inversions are more complex, but can also improve information on the spatial distribution and intensity of fluxes (Lauvaux et al., 2016; Turner et al., 2016); they can be solved by analytical or adjoint methods (Kopacz et al., 2009; Rodgers, 2000). Different cost functions can be used, which might change the results. Here we test and compare three different methods. The first is a Kalman filter (described in Appendix D), which is computationally cheap but has only 1 degree of freedom. For scaling retrievals, using too few degrees of freedom can cause the results to be heavily weighted by the largest model results relative to the observations (Appendix 2). We also use Bayesian inversions based on the methods of Rodgers (2000) (described in Appendix E). One Bayesian inversion is based on a nonlinear forward model with 40 different scaling factors (Eq. E2 in Appendix), and the other is a linear forward model with up to nearly 35 000 scaling factors (Eq. E3), though only a fraction (<1000) of these are used. Because of potential bias in the first two methods, we focus on the linear forward model. Uncertainty estimates are stated for the linear forward model while disregarding the other methods.

2.5 Summary: data sources and methods

In summary, we have four sets of observations of X_gas differences: Caltech TCCON – AFRC TCCON (CO₂, CH₄, and CO) and OCO-2 – AFRC TCCON (CO₂). We use one gridded spatiotemporal inventory for both CO₂ and CO (ODIAC2016, with a weekly pattern for hourly emissions), one for sensitivity tests for CO₂ (Hestia-LA v2.5) and one gridded spatiotemporal inventory for CH₄ (Sect. 2.2). HYSPLIT is run with three dynamical models for the Caltech TCCON – AFRC TCCON differences (GDAS 0.5^∘, NAM 12 km, and HRRR 3 km for a subset) and is run with NAM 12 km for the OCO-2 – AFRC TCCON differences. Three different inversion techniques are used, including a Bayesian inversion with a linear forward model, a Bayesian inversion with a nonlinear forward model, and a Kalman filter. Unless specified, values reported are from the Caltech TCCON – AFRC TCCON difference with the NAM 12 km model and the Bayesian inversion with the linear forward model.

3 Typical X_gas enhancements

Several previous studies have discussed the SoCAB $X_{{\mathrm{CO}}_{\mathrm{2}}}$, $X_{{\mathrm{CH}}_{\mathrm{4}}}$, and X_CO enhancements from local anthropogenic activity (Hedelius et al., 2017a; Janardanan et al., 2016; Kort et al., 2012; Schwandner et al., 2017; Wunch et al., 2009, 2016). There have also been several studies that have discussed enhancements noted from CLARS (California Laboratory for Atmospheric Remote Sensing). CLARS has a viewing geometry that is more sensitive to the mixing layer than TCCON and nadir-viewing satellites, which leads to larger typical enhancements in CO₂ and CH₄ (Wong et al., 2016, 2015). For comparability we exclude enhancements from CLARS and in situ observations (Verhulst et al., 2017) in this section. Kort et al. (2012) noted that observing changes in typical X_gas enhancements from space-based instruments can provide a first-order estimate of how local emissions have changed year to year without the need for a full inversion. This requires similar year-to-year ventilation patterns and sufficiently large and representative sample sizes, which is becoming less of an issue as more space-based observations become available.

Figure 3 Histograms of X_gas enhancements observed in the SoCAB for all dates of this study (Sect. 2.1). Data are averaged for ± 30 min centered on the hour. (a) Enhancements are defined as Caltech TCCON observations minus AFRC TCCON observations (Appendix A). Colors represent the hour of day, and white lines with black dots in the top row are hourly medians. Enhancements peak in early afternoon from morning rush hour emissions being transported from downtown Los Angeles (southwest) to Caltech, and from mixed-layer dynamics. (b) Enhancements are OCO-2 observations minus AFRC TCCON observations. Colors represent the distance from the Caltech TCCON site.

Table 2 lists $X_{{\mathrm{CO}}_{\mathrm{2}}}$ enhancements observed over the SoCAB compared to an external background. An instrument with a smaller footprint (e.g., OCO-2, about 1.3 km × 2.25 km) could observe a wider range of $X_{{\mathrm{CO}}_{\mathrm{2}}}$ enhancements than an instrument with a larger footprint (e.g., GOSAT, about 10.5 km diameter). However, the footprint size should not affect the average enhancement over a domain much larger than an individual footprint. Most enhancements in Table 2 are of the order of 2–3 ppm, except for those reported by Schwandner et al. (2017), which are about double. Though their enhancements are within the range of $\mathrm{\Delta }{X}_{{\mathrm{CO}}_{\mathrm{2}}}$ enhancements in the v7r and v8r histograms in Fig. 3b, they are atypical. Their results are likely atypically large because of dynamics on the two particular dates analyzed and do not include enough data to determine typical enhancements, trends, and source and sink attribution. We disagree with their conclusions that these values are in agreement with Kort et al. (2012) and that TCCON validates this high of a typical SoCAB enhancement. Their conclusion that seasonal variations are 1.5–2 ppm does appear to be supported by previous work (Hedelius et al., 2017a). However, their full attribution of the seasonal cycle to biospheric processes within the basin is not supported by the findings of Newman et al. (2016), who found the excess CO₂ from the biosphere only varied from 8 % (summer) to 16 % (winter) of fossil fuel excess. More likely the changing enhancement reflects a small change in the biosphere and, most importantly, seasonal differences in the basin ventilation.

Kort et al. (2012)Janardanan et al. (2016)Hakkarainen et al. (2016)Hedelius et al. (2017a)Schwandner et al. (2017)

Table 2 SoCAB $X_{{\mathrm{CO}}_{\mathrm{2}}}$ enhancements.

^a Qualitative estimate based on Fig. 1 and Supplement Fig. 3 therein. ^b We modified the boundary condition compared to our previous work (see Appendix A); values are for 14:00 (UTC − 7).

Models that assimilate only global in situ (i.e., no total column) CO₂ data are biased by only about ±1 ppm (1σ ∼1 ppm) compared with TCCON observations (Kulawik et al., 2016). This highlights the need to understand bias and uncertainty in total column observations to the order of a few tenths of a part per million or better to provide new information. The TCCON-predicted bias uncertainty is 0.4 ppm or less (<0.1 %). A long-term CO₂ reduction goal is to reach 20 % of 1990 levels by 2050. This is about a 2 %–3 % decrease per year assuming a constant reduction. Thus a 0.4 ppm bias is of the order of 4–9 years worth of emission reductions.

4 SoCAB flux estimates

4.1 Carbon dioxide

Our flux estimate of CO₂ using the TCCON sites and linear model (Eq. E3) is 104±26 TgCO₂ yr⁻¹. An error assessment is described in Sect. 4.3. This estimate is shown, along with estimates from past studies, in Fig. 4. Our estimate is lower than those from Vulcan (Brioude et al., 2013; Fischer et al., 2017), Hestia-LA v1.0, Hestia v2.5, and the California Air Resources Board (CARB) 2017. Our result is also slightly lower than those of Ye et al. (2017), who estimated emissions by comparing OCO-2 observations with forward model results from a WRF-Chem model. Our result differs significantly from previous TD estimates from aircraft flights, EDGAR v4.0 (Wunch et al., 2009), and CARB 2011. Between 2011 and 2012 CARB changed how bunker fuels and aircraft emissions were reported for the state, which caused a significant decrease in reported emissions. Our posterior estimate is similar to that of EDGAR v4.2 and ODIAC2016, which is slightly less than ODIAC2015. The ODIAC2016 is based on disaggregation of CDIAC national total emissions. Thus, unlike locally developed emission inventories the interannual variations in subnational emissions are driven by the national emission trends. ODIAC could be low from incorrectly distributing too much of the emissions to rural areas due to blooming effects (Small et al., 2005). Blooming effects refer to the tendency for nightlights to exaggerate settlement areas compared with actual extent due to coarse gridded spatial resolution and indirect or nonelectrical light.

Most of the estimates from previous studies include only emissions from fossil fuel use. We have not separately accounted for biospheric uptake (emissions) in the model, and if it is significant, the anthropogenic flux would be larger (smaller) than our net estimate. In the GEOS-Chem model described by Liu et al. (2017) the nearby ocean is a neutral to weak sink, likely from biological activity.

Figure 4 Estimates of SoCAB CO₂ fluxes (annual estimates from TCCON are shown as black triangles and OCO-2 estimates are shown in pink) compared with previous studies. TD aircraft estimates are from Brioude et al. (2013). TD estimate from Ye et al. (2017) is based on OCO-2 observations and 5 % random uncertainty has been added. The Hestia-LA v1.0 estimate was inferred after a forward implementation into a WRF model (Hedelius et al., 2017a). EDGAR v4.2, ODIAC, and CARB emissions were calculated from databases. The EDGAR v4.0 value was reported by Wunch et al. (2009). All other values were found in a literature review. The large range of variability highlights the need for additional study of SoCAB CO₂ fluxes.

OCO-2 provides better spatial coverage than TCCON (Fig. 5), and the orbit tracks can change longitudinally with season or when the spacecraft moves for collision avoidance. However, observations only occur at the same local solar time, and are days to weeks apart. The estimate using OCO-2 data is slightly larger at 120±30 TgCO₂ yr⁻¹, which is in better agreement with the results of Ye et al. (2017). This value varies by up to 12 TgCO₂ yr⁻¹ depending on filtering methods (e.g., warn levels, Appendix 1).

Figure 5 A visualization of OCO-2 observations and the forward model used in the flux inversion on 20 June 2015. The nadir track is shown in red starting at the bottom and $\sim -\mathrm{117.6}$ and going towards the northwest. Observations are overlaid on the green ODIAC prior at 14:00 (UTC − 7). For every fifth sounding the set of back trajectories is shown in gray. Back trajectories originating from the AFRC site are shown in blue. Coastlines and the geopolitical boundaries of the SoCAB are shown in black. County borders are shown in blue.

4.2 CH₄ and CO

Using the same methodology we estimate a CH₄ flux of 360±90 Gg CH₄ yr⁻¹. This is less than the estimate by Wunch et al. (2009) but similar to estimates from Wong et al. (2016) and CARB (Fig. 6). CARB-based CH₄ fluxes for just the SoCAB were estimated by subtracting agriculture and forest emissions (53 %–61 % of total depending on version and year) and out-of-state electricity generation (0 %–0.1 %). The remaining flux was scaled by 42 % based on the population of the SoCAB, and 5 % of the agriculture and forestry emissions were added back in. Our estimate is slightly lower than previous estimates of CH₄ fluxes using in situ (tower and aircraft) data (Cui et al., 2015; Hsu et al., 2010; Peischl et al., 2013; Wecht et al., 2014; Wennberg et al., 2012).

Figure 6 SoCAB CH₄ and CO flux estimates. Annual estimates from this study are shown as light blue triangles. Brioude et al. (2013) and Wecht et al. (2014) (in situ estimates and GOSAT inv) used meteorological models to estimate fluxes, similar to the work presented here. Tracer–tracer relations were used for the TCCON (Wunch et al., 2009, 2016), CLARS (Wong et al., 2016, 2015), and in situ observations (Hsu et al., 2010; Wennberg et al., 2012). “W09 adj” shows the results of Wunch et al. (2009) when adjusted for our posterior CO₂ flux. Methane in situ results for CalNex (May–June 2010) are from Wennberg et al. (2012). CO in situ results in 2002 and 2010 are from Brioude et al. (2013).

We also estimate a CO flux of 487±122 Gg CO yr⁻¹. This is significantly less than the estimates by Wunch et al. (2009) of 1400±300 Gg CO yr⁻¹ from August 2007 to June 2008 and the estimate of 1440±110 Gg CO yr⁻¹ by Brioude et al. (2013) for summer 2010. Wunch et al. (2009) used a tracer–tracer relationship in which the assumed CO₂ was likely too large (191 TgCO₂ yr⁻¹). When their results are scaled down based on our posterior CO₂ fluxes (104–120 TgCO₂ yr⁻¹), the CO flux is 750–880 Gg CO yr⁻¹, which is in better agreement with the CARB inventory. The CARB CO inventories, specific to the SoCAB, have decreasing CO emissions; part of the difference could be from different observation periods. CARB2017 emissions are 581 Gg CO yr⁻¹ for 2015.

4.3 Sensitivity tests and error assessment

For a single estimate of the SoCAB flux, we have a sufficiently large sample that random uncertainty is small. This is supported by a bootstrap analysis in which we select a random subset of data equal in size to the original n=200 times (Efron and Gong, 1983). The random uncertainty estimate is 4 Tg CO₂ yr⁻¹ (2σ), or about 4 %. Persistent biases from a priori flux uncertainty, model errors, observation biases including boundary conditions, and poorly chosen initial values are more detrimental to our flux estimate.

Several variables (x_a, S_ϵ, S_a) need initial values (see Appendix 3), and how these are chosen can affect the final flux calculated. We evaluate four sensitivity tests (Fig. 7). For the first test, we filter out data for which the observations differ from the model above a threshold. We scale from the starting factor of 10× (Appendix 1). We also adjust values of x_a, S_ϵ, and S_a by factors of 2⁻¹⁰ to 2¹⁰. These results show the overall flux generally has low sensitivity to scaling S_ϵ and S_a but has some sensitivity when filtering more data and about a 30 % sensitivity to the scaling of x_a. The interannual variability, which we expect is less than about 25 %, increases for large S_a. Increasing S_a increases r and the degrees of freedom for the signal (dof_s) with only a small effect on the overall flux but also increases the interannual range. Decreasing S_ϵ increases r and dof_s, but it also increases χ² and the interannual range. We estimate an overall uncertainty of 10 % from these parameters.

Figure 7 Assessment of sensitivity to initial values for CO₂. (a) Reduced state vector with seven categories (overall, spatial, vertical, time of day, weekday–weekend, month, and year). λ values indicate how much the prior is scaled on average compared to other elements in its category. Gray lines are results from all tests, and colored lines are from the S_a test. As S_a gets larger, the variability in the retrieved γ factors increases. Missing elements represent lack of sensitivity. (b) Overall fit parameters, including the overall flux, dof_s, χ², Pearson's correlation coefficient between observed and post-inversion model values, and the interannual range. Note the log₂ axes, which indicate the magnitude of change in the sensitivity test compared with the base case (sf: scaling factor). Moving left, filtering (ftr) becomes more stringent, constraints on S_a or S_e are increased, or x_a is scaled down. For scalings less than about 8× the total flux change is small, except for scaling x_a, which increases the flux by about 30 % of the change in the prior. Here the goal was to simultaneously increase dof_s, decrease χ², and increase r while keeping the interannual variability below about 25 %.

We next test the sensitivity to different inversion and modeling schemes (Table 3). The Kalman filter (Appendix D) and the nonlinear inversion (Eq. E2) results are not unreasonable for CO₂. However, their CH₄ flux results are unreasonably low, likely from high model : measured values having unreasonably high weights in these particular schemes with few scaling factors (Appendix 2). GDAS and HRRR results are within uncertainty.

Table 3 Fluxes from various methods.

For a given gas, all the inversions use the same observed ΔX_gas (Caltech TCCON − AFRC TCCON) data. The top three rows are from using different meteorological models, with the same inversion scheme (Eq. E3). The last three rows are from using the same meteorological model (NAM 12 km) with different inversion schemes. Errors are 25 %. ^* For GDAS S_a is 20× smaller.

There is some uncertainty due to the accuracy and resolution of the emission inventories. Gately and Hutyra (2017) compared emission inventories over the northeastern US and noted inventory differences of 100 % for half of the 0.1^∘ grid cells in the domain. Lauvaux et al. (2016) and Oda et al. (2017) compared aggregate posteriori inversion results from different emission inventories and noted differences of only 5 %–8 % for the Indianapolis region despite large differences at the grid level. We make a similar comparison in which we use the more spatially and temporally accurate Hestia v2.5 fossil fuel inventory instead of ODIAC as the prior. We note that the correlation between the forward model data and TCCON is slightly higher with Hestia than ODIAC, and there are fewer outliers that differ by a factor of 10× or more. However, the flux estimate of 110±28 is similar to the posterior flux estimate using ODIAC.

Finally we consider the observation uncertainty. Hedelius et al. (2017b) reported a 2σ measurement bias of less than ∼0.2 ppm $X_{{\mathrm{CO}}_{\mathrm{2}}}$ (central estimate, maximum range <0.5 ppm) between the AFRC and Caltech TCCON sites, but even a bias of 0.2–0.3 ppm $X_{{\mathrm{CO}}_{\mathrm{2}}}$ will produce an error of ∼ 10 % in the flux. This bias could also arise from improper boundary conditions or application of averaging kernels.

In summary, we estimate 5 % random uncertainty from the bootstrap analysis, 10 % from our choice of initial values, 5 % from the prior flux, 10 % from observations and the boundary condition, and 20 % from model winds (Appendix B). The sum in quadrature is 25 %. By comparison, uncertainty estimates from other inversions were 11 % (inner 50 percentile range from an ensemble) for Indianapolis (Lauvaux et al., 2016) and 5 % for the Bay Area using pseudo-observations (Turner et al., 2016). Both of these studies benefited from additional sites (nine and 34, respectively) and custom WRF model runs. Ye et al. (2017) estimated an uncertainty of 5 % for the SoCAB flux by using data from 10 OCO-2 tracks; however this is not directly comparable with our result because it does not include uncertainty from biases in the forward model, observations, and inversion scheme.

5 Discussion

5.1 Emission ratios

Emission ratios can help us evaluate the inversion for the SoCAB. In previous studies it was noted that the Pasadena area is a good receptor site for the basin, so tracer–tracer ratios observed there should approximately correlate with emission ratios (Newman et al., 2016; Wunch et al., 2009, 2016). If the ratios are significantly different it could highlight an error in the inversion scheme or the a priori assumption of sources. However, errors in the model can be correlated for different tracers, which would obscure universal biases to all gases. For example, the CO : CO₂ flux ratio of 7.5 ppb : ppm is in good agreement with past literature, despite the absolute fluxes being lower on average.

We estimate emission ratios using the solar zenith angle (SZA) anomaly method described by Wunch et al. (2009, 2016) and using the average enhancement compared with AFRC or the Pasadena : Lancaster gradient ratio. Errors are assumed to equal the standard deviation of all the data, and a linear fit is made using the methods of York et al. (2004) on monthly timescales. We estimate the emission ratio from the work of Verhulst et al. (2017) using the weighted mean of the excess ratios from their five in-basin sites, with weights of $\frac{\mathrm{1}}{{\mathit{\sigma }}^{\mathrm{2}}}$. Emission ratios from the SZA anomaly method and the differenced enhancement are in agreement with ratios from previous studies (Fig. 8). The CH₄ : CO₂ from the inversion at 9.7 ppb : ppm is larger than past studies, suggesting either our CH₄ flux is too large or our CO₂ flux is too low (or some combination of both). If both were adjusted by 15 % the ratio would be 7.1 ppb : ppm. The CO : CO₂ ratio is in good agreement with the ratios using the SZA anomaly method, but is lower than past estimates. Based on the CARB inventories, a decrease is expected because CO emissions have decreased almost exponentially over the past decades, whereas CO₂ emissions have decreased only moderately (e.g., compare Figs. 4 and 6).

Figure 8 Emission ratios compared with previous studies. Numbers shown are central values from the different methods and studies. Overall fits are shown as dashed lines. The ratio from the gradients of Pasadena (Caltech) to Lancaster (AFRC) is based on the difference between the two TCCON sites. Values from Silva et al. (2013) and Silva and Arellano (2017) were part of global studies.

In November 2015, the large CH₄ : CO₂ ratio is from additional methane emissions from the Aliso Canyon gas leak (Conley et al., 2016). Though this leak persisted until February 2016, different wind patterns caused less of the highly methane-enriched air to be transported and observed in Pasadena after the first 2 months. The large CO : CO₂ ratios seen in summer 2016 are from wildfires. The San Gabriel Complex Fire was less than 25 km to the east and burned 22 km² over a month. It was close enough for ash to be transported to Pasadena. Eight other major fires within 150 km burned an additional 400 km² during June–August 2016 (http://cdfdata.fire.ca.gov/incidents/incidents_archived?archive_year=2016, last access: 12 November 2018).

Pollack et al. (2012)Brioude et al. (2013)

Table 4 Weekday : weekend emission ratios.

^a WD : WE CO ratios from Pollack et al. (2012) were calculated using the difference between the CalNex-Pasadena and c Flasks (Table 2 therein). For CO₂ we used the CO WD : WE ratios with the CO : CO₂ WD : WE ratios in Pasadena (Table 3 therein). ^b WD : WE ratios from Brioude et al. (2013) were calculated by assuming ratios between daytime and all-day emissions in the posterior were equal using Table 3 therein. ^c Temporal Improvements for Modeling Emissions by Scaling (TIMES) were reported for the contiguous United States by Nassar et al. (2013). ^d Hestia-LA is based on Fig. 2 from Hedelius et al. (2017a). This same ratio is used in the CO and CO₂ priors in this study.

5.2 Weekend effect

The weekday-to-weekend (WD : WE) flux ratios are listed in Table 4. The uncertainty is estimated to be ±0.10 based on changes in the ratio from the S_a scaling test up to 16× (Fig. 7). Weekday : weekend ratios are similar to those from previous studies for CO (Brioude et al., 2013; Pollack et al., 2012). The CO₂ and CO ratios are scaled down compared with the prior, which puts CO in better agreement with Hestia-LA v2.5. Methane has a ratio that is slightly larger than unity. Methane is not expected to vary as much as CO₂ or CO on weekdays compared to weekends because production from biogenic sources and fugitive losses from natural gas infrastructure are less time variant.

6 Conclusions

This study demonstrates a method to readily obtain estimates of net CO₂ fluxes over regions of the order of 10 000 km using remote-sensing observations. This work is a step towards estimating fluxes from a greater number of urban areas using space-based observations of CO₂. Our estimates of total annual CO₂ fluxes from the SoCAB using HYSPLIT with NAM 12 km as our dynamical model are on the low end of previous estimates (Fig. 4) and about 28 %–47 % less than inventory values reported in tracer–tracer flux estimate papers (Wong et al., 2015; Wunch et al., 2009). This has important implications for these studies, which would have overestimated CH₄ emissions if CO₂ emissions were also too large. Net CO fluxes are significantly less than previous studies, likely from an underestimate of about 20 % combined with a known decrease in emissions. Net CH₄ fluxes are in agreement with previous studies.

This study is one of only a few in which satellite observations were used to help infer the net flux of CO₂ from an urban area. Several lessons learned here will be important for future studies using space-based observations of CO₂ for flux estimates from other urban regions. We have shown a method for accounting for the sensitivity of the instrument to true changes in the atmospheric composition (i.e., accounting for averaging kernels). We have also shown a method to account for differences in column observations that could arise from different surface altitudes. In the Appendix we document how changing the number of elements in a retrieval state vector can, in some cases, bias the inferred flux result. This effect becomes increasingly important for inversions using only one scale factor with a large discrepancy between the forward model and observations. Finally, we describe the sensitivity of the results to filtering and parameters such as the a priori and the a priori covariance matrix S_a.

The overall uncertainty is 25 %, with the dynamical model contributing the most. X-STILT (Wu et al., 2018) and higher-resolution models (e.g., HRRR) may help reduce dynamical model uncertainty in future inversions. We consider an uncertainty of 25 % to be large and it shows additional work is needed to improve constraints. If errors are from persistent biases, then relative changes in time can be observed, though such changes might also be observed using just the observations without a model (Kort et al., 2012). Understanding contributions from the biosphere may also be important in future studies to diagnose how much carbon is from fossil fuels (Newman et al., 2016). The wide range of uncertainty suggests that CO₂ flux estimates from the SoCAB will benefit from additional measurements – such as the LA Megacity Carbon Project in situ tower network (Verhulst et al., 2017), the planned geostationary GeoCARB mission, and the OCO-3 mission, which has a raster mode that can scan throughout the basin. Further improvements in modeling and inversion techniques will also help, including assimilating all available observations (in situ network, TCCON, CLARS, OCO-2, and GOSAT). These additional surface and space-based observations can not only aid in improving the accuracy of the overall flux, but also may be incorporated into spatiotemporal inversions to map fluxes from subregions of the SoCAB with confidence.

Appendix A: Observation data filtering and boundary condition

GOSAT–ACOS v2.9 $X_{{\mathrm{CO}}_{\mathrm{2}}}$ levels are enhanced by only 3.2 ± 1.5 (1σ) ppm in the SoCAB (Kort et al., 2012). This means a bias of 0.3 ppm could lead to a 10 % bias in the flux. Thus it is critical to account for biases down to the tenths of a part per million level or better. This is a challenge given that the accuracy of OCO-2 (v7r) over land had been estimated as 0.65 ppm (Worden et al., 2017), and OCO-2 comparisons with TCCON range from −0.1 to 1.6 ppm (Wunch et al., 2017).

A1 Quality filters

Compared with the TCCON, OCO-2 spectra have a lower resolution. OCO-2 observations are also sensitive to surface albedo and are more sensitive to aerosol scattering than solar-viewing instruments. These sensitivities can cause spurious results, which need to be filtered out. Included in the OCO-2 data is a binary flag as well as warn levels (WLs) for quality filtering. WLs are a global metric of data quality, for which WLs less than or equal to (0, 1, 2, 3, 4, 5) correspond to about (50 %, 60 %, 70 %, 80 %, 90 %, 100 %) of data passing in v8r, and larger WLs generally correspond to less reliable data. WL definitions are different for v7 and v8, but here we use the binary $X_{{\mathrm{CO}}_{\mathrm{2}}}$ filter and only include v8 data with a WL ≤1. WL ≥4 data are already removed by the binary flag. We also exclude data that differ from the model by a factor of 10 or more. This factor of 10 is somewhat arbitrary and an argument could be made against using this criterion as a filter. However, a few large outliers can significantly affect inversion results (Appendix 2) so we opt to remove suspect values. A sensitivity test including different filter cutoffs for TCCON $X_{{\mathrm{CO}}_{\mathrm{2}}}$ is described in Sect. 4.3. After filtering, 2361 paired OCO-2–AFRC observations remain.

For TCCON observations we use the public data, which already have some static within-range filters applied. We also exclude data that differ from the model by a factor of 10 or greater, leaving 4872 observations.

A2 Background, boundary conditions, and averaging kernels

To eliminate the ambient X_gas levels that would be observed in the absence of local emissions, we subtract values measured by the AFRC TCCON site from both the Caltech TCCON and OCO-2 data obtained in the basin. We choose TCCON data as background for OCO-2 to reduce the likelihood of albedo-related bias from using OCO-2 observations over the Mojave Desert (Wunch et al., 2017) as well as the chance of inducing a bias from using different viewing modes by using ocean glint observations. In other studies of $X_{{\mathrm{CO}}_{\mathrm{2}}}$ enhancements, observations not directly influenced by the source were used as background. For example, Janardanan et al. (2016) categorized space-based observations of $X_{{\mathrm{CO}}_{\mathrm{2}}}$ by making a forward model estimate of $X_{{\mathrm{CO}}_{\mathrm{2}}}$ enhancements from fossil fuel combustion and setting a threshold to define as polluted or unpolluted. Such an approach could work globally but may have errors from the prior emissions or transport model.

Because we expect most of the difference in $X_{{\mathrm{CO}}_{\mathrm{2}}}$ to arise from polluted air near the surface, we divide the enhancements by the surface averaging kernels of the in-basin observations. OCO-2 surface averaging kernels in the basin are 0.986 ± 0.010 (1σ) with a 99 % confidence interval of 0.955 to 1.016. TCCON surface averaging kernels depend on surface pressure and solar zenith angle (SZA) and are 0.96 ± 0.14 (1σ) throughout the full range of observations.

Even in the absence of local anthropogenic emissions, the $X_{{\mathrm{CO}}_{\mathrm{2}}}$ measured within the SoCAB could be different from that measured at AFRC by a few tenths of a part per million because of different measurement heights and atmospheric CO₂ profiles (Hedelius et al., 2017a). We account for a boundary condition of the form

$$\begin{array}{}\text{(A1)}& b_{\mathrm{gas}}=\left({\displaystyle \frac{X_{\mathrm{gas},\mathrm{a},S}}{X_{\mathrm{gas},\mathrm{a},B}}}-\mathrm{1}\right){\widehat{X}}_{\mathrm{gas},B},\end{array}$$

where subscript “a” represents the a priori estimate, S represents a measurement within the SoCAB, B represents the background, and the circumflex represents a retrieved value. Equation (A1) can be interpreted as the difference that would be observed between two sites due to differences in the gas vertical profiles. The a priori profiles do not include local anthropogenic emissions. The result from Eq. (A1) is subtracted from the SoCAB–AFRC difference. We perform the same adjustment for CH₄ and CO.

Appendix B: Dynamical model error

Figure B1 Histogram of wind speed errors (HYSPLIT minus measured) compared to surface observations at the San Gabriel airport. The mean error is −0.9 m s⁻¹ (−30 %), the 95 % confidence interval is [−3.7, 2.3] m s⁻¹. The mean direction error is less than 5^∘, and 75 % of direction errors are within $\pm \mathrm{45}{}^{\circ }$.

Dynamical models can have errors in the PBL height estimation as well as in the wind speed and direction. In a case study for spring 2011 and 2012 primarily over the Midwestern US, a NAM temperature-derived PBL height had a mean bias of about −50 m, with an inner 50 percentile range of about ±250 m (Coniglio et al., 2013). For wind error we compare with 10 m winds from the San Gabriel (El Monte) Airport 10 km SE of Caltech (lat 34.083, long −118.033, 90 m a.s.l.). We assume the winds are the same at both locations. Airport meteorological data are obtained through the NOAA National Centers for Environmental Information (https://www.ncdc.noaa.gov/cdo-web/datatools/lcd, last access: 12 November 2018).

Trajectory speed and direction are estimated based on when and where trajectories ending at 50 m a.g.l. enter a 5 km radius circle around the receptor site. Results are shown in Fig. B1. The mean speed of HYSPLIT trajectories is less than what is expected by comparing with the surface winds. In contrast, previous studies have shown high model wind speed bias near the surface at the LAX airport, 34 km to the southwest, and near the coast (Angevine et al., 2012; Feng et al., 2016; Ye et al., 2017). The difference biases could in part be from the coast versus inland; however, Feng et al. (2016) also showed a high model bias closer to Caltech. Model differences, the 10 km horizontal and ∼150 m height difference between Caltech and the airport could also contribute to the discrepancy. We expect the average bias throughout the PBL to be lower than at the surface and assign an uncertainty of up to ∼20 % to the average wind.

Appendix C: Residence times from HYSPLIT

HYSPLIT mean trajectories are air parcel locations at different heights for select times (in our case, every 20 min). These are aggregated and normalized for each 0.01^∘ × 0.01^∘ cell and for each hour. First each trajectory is interpolated to 1 s positions. Then we determine the vertical fraction of the mixing layer the trajectory takes. This fraction is the ending vertical spacing between adjacent trajectories (hPa) divided by the local mixing-layer depth (hPa). The HYSPLIT mixing depth is based on the underlying Eulerian model. Parcels above the mixing layer get counted as zero. Then we count how long any parcel was in each cell (s) to obtain the residence time. Monthly average examples of this are shown in Fig. C1. The residence time is multiplied by the a priori flux to determine the column enhancement (g m⁻²). By dividing by a model estimate of the dry air (molecules m⁻²) based on model surface pressure, we obtain a forward model estimate of the X_gas enhancement from local sources (ppm or ppb).

Figure C1 Maps of monthly averaged residence times in the mixing layer per pixel for trajectories ending at 21:00 UTC, shown for all times leading up to the observation. Pixels are $\mathrm{0.01}{}^{\circ }\times \mathrm{0.01}{}^{\circ }$, or approximately 1.03 km². In July the origins were more predictable, but in January there was greater variation. Coastlines and the geopolitical boundaries of the SoCAB are shown in black. County borders are shown in blue.

Appendix D: Kalman filter

The Kalman filter used to estimate SoCAB CO₂ emissions is based on methods described by Kleiman and Prinn (2000) with modifications. This is an iterative approach using a single overall scaling factor. The difference in $X_{{\mathrm{CO}}_{\mathrm{2}}}$ between measurements in the SoCAB and AFRC is the observed measurement, y^obs. The error σ_y associated with each y^obs is estimated from the sum in quadrature of the error from each site, i.e.,

$$\begin{array}{}\text{(D1)}& {\mathit{\sigma }}_y=\sqrt{{\stackrel{\mathrm{‾}}{\widehat{y}}}_{\mathrm{err},\mathrm{C}}^{\mathrm{2}}+{\stackrel{\mathrm{‾}}{\widehat{y}}}_{\mathrm{err},\mathrm{A}}^{\mathrm{2}}},\end{array}$$

where the subscript C is for Caltech (or measurements in the SoCAB), and A is for AFRC (or “background”). The error of the averaged data for an individual site is estimated as

$$\begin{array}{}\text{(D2)}& {\stackrel{\mathrm{‾}}{\widehat{y}}}_{\mathrm{err}}=\sqrt{{\displaystyle \frac{\mathrm{1}}{n-\mathrm{1}}}{\displaystyle \frac{\mathrm{1}}{\sum {\widehat{z}}_{i,\mathrm{err}}^{-\mathrm{2}}}}\sum _i{\widehat{z}}_{i,\mathrm{err}}^{-\mathrm{2}}{\left({\widehat{z}}_i-\stackrel{\mathrm{‾}}{\widehat{y}}\right)}^{\mathrm{2}}},\end{array}$$

where n is the number of measurements, ${\widehat{z}}_i$ are the individual $X_{{\mathrm{CO}}_{\mathrm{2}}}$ measurements, ${\widehat{z}}_{i,\mathrm{err}}$ are the reported errors associated with the measurements, and $\stackrel{\mathrm{‾}}{\widehat{y}}$ is the weighted average using ${\widehat{z}}_{\mathrm{err}}^{-\mathrm{2}}$ as weights. Note that Eq. (D2) takes into account both the measurement errors and the spread of the measurements. However, we note that a similar equation underestimated the error compared with a bootstrap method (Gatz and Smith, 1995).

D1 Iterations

We initialize the iterations with an arbitrary scaling factor α₀=1 and an associated error of σ_α 0=0.7. These initial values have little influence on the final result.

We iterate over the k measurements by calculating the partial derivative:

$$\begin{array}{}\text{(D3)}& h_k={\displaystyle \frac{\partial y_k^{\mathrm{est}}}{\partial {\mathit{\alpha }}_k}}=\sum _j{s}_{j,k}{t}_{j,k},\end{array}$$

where subscript j is for a particular grid box, s is the a priori surface flux, and t is the residence time. Equation (D3) is identical to Eq. A2 in Kleiman and Prinn (2000). Because this is a scaling retrieval, h_k is the observation operator. We can multiply it by the state element (α) to obtain the estimated observation (Eq. D6). The gain scalar g_k and new state error are calculated by (Eqs. A4 and A5 in Kleiman and Prinn, 2000)

$$\begin{array}{}\text{(D4)}& {\displaystyle }& {\displaystyle }{g}_k={\mathit{\sigma }}_{\mathit{\alpha },k-\mathrm{1}}^{\mathrm{2}}{h}_k{\left(h_k^{\mathrm{2}}{\mathit{\sigma }}_{\mathit{\alpha },k-\mathrm{1}}^{\mathrm{2}}+{\mathit{\sigma }}_{z,k}^{\mathrm{2}}\right)}^{-\mathrm{1}},\text{(D5)}& {\displaystyle }& {\displaystyle }{\mathit{\sigma }}_{\mathit{\alpha },k}^{\mathrm{2}}={\mathit{\sigma }}_{\mathit{\alpha },k-\mathrm{1}}^{\mathrm{2}}\left(\mathrm{1}-h_k{g}_k\right).\end{array}$$

We make a modification to calculate the estimated measurement, omitting the term for the convergence of fluxes due to unresolved motions in the transport model. The estimated forward model is

$$\begin{array}{}\text{(D6)}& y_k^{\mathrm{est}}={\mathit{\alpha }}_{k-\mathrm{1}}{h}_k,\end{array}$$

and the state estimate is

$$\begin{array}{}\text{(D7)}& {\mathit{\alpha }}_k={\mathit{\alpha }}_{k-\mathrm{1}}+g_k\left(y_k^{\mathrm{obs}}-y_k^{\mathrm{est}}\right).\end{array}$$

D2 A note on single scale factor inversions with large outliers

Some single scale factor inversions can be written in the form

$$\begin{array}{}\text{(D8)}& \widehat{\mathit{y}}=\mathit{\lambda }{\mathit{y}}_{\mathrm{mod}},\end{array}$$

where “mod” represents the initial model values. We consider the case when cost function is of the form

$$\begin{array}{}\text{(D9)}& J_c=\sum _k{\displaystyle \frac{{\left(y_{\mathrm{obs},k}-\mathit{\lambda }{y}_{\mathrm{mod},k}\right)}^{\mathrm{2}}}{y_{\mathrm{err},k}^{\mathrm{2}}}},\end{array}$$

where the error term accounts for both the model and observation errors. If the error is not a function of λ then

$$\begin{array}{}\text{(D10)}& {\displaystyle \frac{\partial J_c}{\partial \mathit{\lambda }}}=\mathrm{2}\sum _k{\displaystyle \frac{\mathit{\lambda }{y}_{\mathrm{mod},k}^{\mathrm{2}}-y_{\mathrm{obs},k}{y}_{\mathrm{mod},k}}{y_{\mathrm{err},k}^{\mathrm{2}}}}.\end{array}$$

Setting Eq. (D10) equal to zero and solving for λ yields

$$\begin{array}{}\text{(D11)}& \mathit{\lambda }={\displaystyle \frac{{\mathrm{\Sigma }}_k{y}_{\mathrm{obs},k}{y}_{\mathrm{mod},k}}{{\mathrm{\Sigma }}_i{y}_{\mathrm{mod},k}^{\mathrm{2}}}}.\end{array}$$

Note the change from vector to summation notation. Equation (D11) is a first-order estimate of the overall scale factor λ. This indicates that λ can be low with high model : observation ratios, which heavily weight the result.

Figure D1 Effects of scaling a random subset of model data compared to no scaling. When fewer points are scaled, they are scaled by a larger amount. The Kalman filter and nonlinear inversion are more affected by a few strong outliers than the linear inversion.

This is demonstrated in a sensitivity test, in which we scale a subset of points (Fig. D1). We create pseudo-observed values by using the original model values. We create pseudo-model data by scaling a random subset of the original model data by $(\mathrm{1.1}{)}^{\frac{n}{s}}$, for which n is the total number of points, and s is the number in the subset. For example, when 100 % of the model points are adjusted we scale them all up by 10 %. The test is repeated multiple times, with fewer repeats for the nonlinear model because it takes the longest. These results show that having a few large outliers in the Kalman filter (one scale factor) and the nonlinear (40 scale factors, Sect. E) inversions can significantly pull the results compared with the linear (∼1000 scale factors, Sect. E) inversion.

Appendix E: Bayesian inversions

The Bayesian approach to solving atmospheric inverse problems has been described in more detail by Rodgers (2000). Turner et al. (2016) describes this approach for an urban region. Here we follow the notation of Rodgers (2000). For scaling retrievals, Bayesian inversions minimizes a cost function ($-\mathrm{2}\ln P\left(\mathit{y}\right|\mathit{x})$) of the form in Eq. (D9). This assumes error statistics are adequately known and are Gaussian for both the state vector x (length n) and the measurement vector y (length m).

Figure E1 Extent of the eight spatial subregions. C: center; Q1–Q4: SoCAB quadrants; O: ocean; CV: Central Valley, D: all other areas, mostly the Mojave Desert to the northeast.

E1 Forward model

The generalized forward model can be written as

$$\begin{array}{}\text{(E1)}& \mathit{y}=\mathbf{F}\left(\mathit{x}\right)+\mathit{ϵ},\end{array}$$

where ϵ is an error term. We test two similar forward models for the Bayesian inversions. The first is chosen to reduce the number of elements in the state vector. This choice was made based on having only two measurement locations. This model has 40 state vector scaling factors λ in seven different classes corresponding to year (6), month (12), weekday–weekend (2), time of day (6), vertical level (5), spatial bin (8), and overall (1). Time-of-day bins cover 4 h each with local ending times at 03:00, 07:00, 11:00, 15:00, 19:00, and 23:00. Aggregated vertical bins are each about 3.5 % of the atmosphere, split at 300, 612, 936, 1272, and 3200 m a.g.l. These are designed to help diagnose transport or footprint extent errors, and the upper two levels are weighted less when estimating the total SoCAB flux. Spatial bins (Fig. E1) were chosen with one over the ocean, one over Central Valley, one for the rest of the area outside the SoCAB, and five inside the SoCAB. Each SoCAB area has approximately the same influence on observations at the Caltech (abbreviated CIT) site based on residence times. This model is

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\mathbf{F}\left(\mathit{x}\right)={\mathit{\lambda }}_{all}\sum _{\mathrm{yr}=\mathrm{1}}^{\mathrm{6}}\sum _{\mathrm{mth}=\mathrm{1}}^{\mathrm{12}}\sum _{\mathrm{dow}=\mathrm{1}}^{\mathrm{2}}\sum _{\mathrm{tod}=\mathrm{1}}^{\mathrm{6}}\sum _{\mathrm{vbin}=\mathrm{1}}^{\mathrm{5}}\sum _{\mathrm{sbin}=\mathrm{1}}^{\mathrm{8}}\\ \text{(E2)}& {\displaystyle }& {\displaystyle }{\mathit{\lambda }}_{\mathrm{yr}}{\mathit{\lambda }}_{\mathrm{mth}}{\mathit{\lambda }}_{\mathrm{dow}}{\mathit{\lambda }}_{\mathrm{tod}}{\mathit{\lambda }}_{\mathrm{vbin}}{\mathit{\lambda }}_{\mathrm{sbin}}\left(m_{j,\mathrm{CIT}}-m_{j,\mathrm{AFRC}}\right).\end{array}$$

Here, $m=\mathrm{\Sigma }t\times s$ is the model amount determined by multiplying the residence time t by the a priori surface flux s and summing over all times and $\mathrm{0.01}{}^{\circ }\times \mathrm{0.01}{}^{\circ }$ grid boxes in the bin. We use j here as shorthand for the subscript “yr, mth, dow, tod, vbin, sbin”.

We also use a similar linear model of the form

$$\begin{array}{}\text{(E3)}& \mathbf{F}\left(\mathit{x}\right)=\sum _{j=\mathrm{1}}^{\mathrm{34}\mathrm{560}}{\mathit{\lambda }}_j\left(m_{j,\mathrm{CIT}}-m_{j,\mathrm{AFRC}}\right).\end{array}$$

In this form there are up to nearly 35 000 original elements in our state vector as opposed to the 40 elements in Eq. (E2). Most of the original elements are not ever sampled (e.g., during 2012 and 2017) and not used when reporting our total fluxes. We remove elements that are not linearly independent, which reduces the actual number used to fewer than (about one-eighth) the number of observations. We select the most important elements from matrix R found by performing a QR decomposition on the K matrix. Changing the cutoff (and S_a) affects the sensitivity to the prior (Sect. 3).

E2 Solutions

For the linear forward model (Eq. E3), the retrieved state vector ($\widehat{x}$) can be found in a single step,

$$\begin{array}{}\text{(E4)}& \widehat{\mathit{x}}={\mathit{x}}_{\mathrm{a}}+{\mathbf{S}}_{\mathrm{a}}{\mathbf{K}}^T{\left({\mathbf{KS}}_{\mathrm{a}}{\mathbf{K}}^T+{\mathbf{S}}_{\mathit{ϵ}}\right)}^{-\mathrm{1}}\left(\mathit{y}-\mathbf{K}{\mathit{x}}_{\mathrm{a}}\right).\end{array}$$

x_a denotes the a priori state vector. S_a is the a priori covariance matrix for the state vector (denoted B in some texts). K is the m×n Jacobian matrix (denoted H in some texts). S_ϵ is the m×m measurement error covariance matrix (denoted R in some texts), which includes errors from both the observations and the forward model. S_ϵ is often treated as a diagonal matrix, with ${\mathit{\sigma }}_k^{\mathrm{2}}$ values along the diagonal.

For a nonlinear forward model (e.g., Eq. E2), the inverse solution can be found using an iterative Levenberg–Marquardt method. This is described in more detail by Rodgers (2000). The iterative solution is

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{\mathit{x}}_{i+\mathrm{1}}={\mathit{x}}_i+{\left[\left(\mathrm{1}+\mathit{\gamma }\right){\mathbf{S}}_{\mathrm{a}}^{-\mathrm{1}}+{\mathbf{K}}_i^T{\mathbf{S}}_{\mathit{ϵ}}^{-\mathrm{1}}{\mathbf{K}}_i\right]}^{-\mathrm{1}}\\ \text{(E5)}& {\displaystyle }& {\displaystyle }\left\{{\mathbf{K}}_i^T{\mathbf{S}}_{\mathit{ϵ}}^{-\mathrm{1}}\left[\mathit{y}-\mathbf{F}\left({\mathit{x}}_i\right)\right]-{\mathbf{S}}_{\mathrm{a}}^{-\mathrm{1}}\left[{\mathit{x}}_i-{\mathit{x}}_{\mathrm{a}}\right]\right\}.\end{array}$$

The symbol γ is a factor chosen at each iteration to minimize the cost function based on how χ² changes, and i+1 denotes the current iteration.

E3 A priori values

We define values for x_a, S_ϵ, and S_a. First, our state vector is composed of scaling factors and all elements in x_a are unity for CO and CH₄. Because ODIAC2016 emissions for the SoCAB are low compared to other inventories, we multiply by 1.25 for CO₂. S_ϵ is a diagonal matrix. Along diagonal elements are the errors from the observations plus the errors from the transport model. Observation errors are ${\mathit{\sigma }}_y^{\mathrm{2}}$ determined from Eq. (D1). We assume transport errors are constant and equal to the overall median observation error.

For simplicity, S_a is chosen as a single scalar value for the linear model (Eq. E3). We tune two parameters, namely S_a, and the threshold for determining linear independence in the QR decomposition. This is a trade-off between maximizing the degrees of freedom and r, while avoiding unstable conditions, and minimizing χ². We scan over a variety of S_a and threshold values. We use interannual variability and dependence on the prior as noted by a sensitivity test (Fig. 7) to judge the quality. Generally as we increase the threshold fewer elements are allowed in the state vector, the dependence on the prior decreases, and the interannual range increases. As S_a increases, so does the interannual range, and the dependence on the prior decreases. We select values that keep the interannual variability under about 25 % and minimize dependence on the prior. We repeat this procedure for the three gases retrieved by TCCON and for OCO-2 observations. S_a is tuned to 0.01 for CO₂, 0.007 for CH₄, 0.0007 for CO, and 0.04 for CO₂ using OCO-2 observations. For the 40-factor inversion, S_a is a matrix and diagonal values are the same as the linear inversion. Off-diagonal values between adjacent elements (e.g., years, months) are one-third of those along the diagonal, which is a somewhat arbitrary choice based on our a priori guess of how strongly adjacent elements are related.