Southern California Megacity CO 2 , CH 4 , and CO flux estimates using remote sensing and a Lagrangian model

We estimate the overall CO2, CH4, 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 CO2 (ODIAC). Using TCCON data we estimate the direct net CO2 flux from the SoCAB to be 139± 35 Tg CO2 yr−1 for the study period of July 2013–August 2016. We obtain a slightly lower estimate of 118± 29 Tg CO2 yr−1 using OCO-2 data. These 5 CO2 emission estimates are in general agreement with previous work. Our net CH4 (325± 81 Gg CH4 yr−1) flux estimate is slightly lower than central values from previous top-down studies going back to 2010 (342–440 Gg CH4 yr−1). CO emissions are estimated at 555± 136 Gg CO yr−1, much lower than previous top-down estimates (1440 Gg CO yr−1). 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 inversions schemes or a coarser dynamical model influence the emission estimates. Overall, 10 the uncertainty is estimated to be 25 %, with the largest contribution from the dynamical model. The methods described are scalable and can be used to estimate direct net CO2 fluxes from other urban regions.


Introduction
About 43 % of global anthropogenic carbon dioxide (CO 2 ) emissions come directly from urban areas, and urban final energy use accounts for about 76 % of CO 2 emissions (Seto and Dhakal, 2014).Associations of cities that recognize their significant emissions of CO 2 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.servations on 29 different overpass days when the AFRC TCCON site also collected background observations before filtering.
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.

A priori flux estimates
Our flux estimate involves scaling the a priori spatial inventory, or sub-regions of the prior up or down to reduce the measured−model mismatch.More important than the total prior absolute flux is the distribution of sources.Hestia-LA v2.0 is likely the most accurate spatiotemporal inventory for the SoCAB, however it is not available globally.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 2 (ODIAC2016) which is available globally at a resolution of 30 arcseconds from 2000-2015 (Oda andMaksyutov, 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.For the Indianapolis region, Lauvaux et al. (2016) noted little difference in the aggregate inversion flux when using ODIAC compared with Hestia.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 0.25 • × 0.25 • can be applied to spa- tial inventories such as ODIAC to improve temporal emissions (Nassar et al., 2013).However, TIMES has a single peak for mid-day 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 downscale the ODIAC to a 0.01 • ×0.01 • grid over the domain 121.5 • W-114.5 • W and 30.5 • N-37.5 • N.This same prior is used for CO, but total emissions are 1 % of CO 2 emissions on a molar basis (0.6 % of mass).Figure 1 shows the ODIAC2016 prior for one month.
We make our own 0.01 • ×0.01 • methane prior using landfills, nightlights, expected total emissions, and the Harvard-Environmental Protection Agency (EPA) United States (U.S.) inventory (Maasakkers et al., 2016) shown in Fig. 1.A more detailed CH 4 inventory is also available for the SoCAB, which we do not use because it would be difficult to scale globally (Carranza et al., 2018).First, we distribute emissions from landfills as point sources (available 2010-2015, https://ghgdata.epa.gov/ghgp/main.do) and use 2015 emissions for 2016.Emissions from the Puente Hills landfill were doubled because the EPA estimate (average 13.6 Gg CH 4 yr −1 ) is low compared to previous estimates of 34 Gg CH 4 yr −1 (Peischl et al., 2013).After doubling Puente Hills emissions, EPA total SoCAB (144 Gg CH 4 yr −1 ) and Olinda Alpha (13.5 Gg CH 4 yr −1 ) landfill emissions are similar enough to other studies (164 Gg CH 4 yr −1 and 12.5 Gg CH 4 yr −1 respectively, Peischl et al., 2013) that we do not double emissions from other landfills in the SoCAB.Chino dairy emissions were added in as a ∼ 0.1 • × 0.1 • source (Chen et al., 2016;Viatte et al., 2017).Outside of the SoCAB CH 4 manure and enteric fermentation were added from the 0.1 • × 0.1 • Harvard-EPA inventory (Maasakkers et al., 2016).SoCAB emissions are assumed to sum to 400 Gg CH 4 yr − 1 based on the  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 arcseconds).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 4 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 2 emissions from the SoCAB.Literature estimates range from 99 Tg CO 2 yr −1 (Vulcan, Fischer et al., 2017) to 211 Tg CO 2 yr −1 (EDGAR v4.0, as reported by Wunch et al., 2009).Table 1 lists statistics for the SoCAB.We assume the non-residential natural gas (NG) use is for industry or power 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 U.S. and Central Valley (CV), there is a CO 2 return flux to the croplands from both human respiration and food waste.In the U.S. 60 million metric tonnes (MMT) of food are lost annaully 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 2 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 4 :CO 2 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 2 d −1 person −1 (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.Based on these various statistics we estimate a bottom-up net flux on the order of 110 Tg CO 2 yr −1 from the SoCAB.

Dynamical models
A dynamical model is needed in conjunction with the a priori flux estimates to generate forward model X gas enhancements.
This may be as complex as a custom high-resolution Weather Research and Forecasting (WRF) model (e.g., Lauvaux et al., 2016) or as simple as an average mixed layer wind velocity (e.g., Chen et al., 2016).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 an 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 measured−model mismatch for the regions of interest.
We use the North American Mesoscale Forecast System (NAM) at 12 km resolution (3 hr 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 began publicly releasing 3 km, 1 hr archived data from the High Resolution Rapid Refresh (HRRR) model that covers the U.S. (Benjamin et al., 2016).This product holds the potential to improve flux estimates at smaller scales.
We use HYSPLIT-4 (Hybrid Single Particle Lagrangian Integrated Trajectory-4;Stein et al., 2015) with the 3 archived NOAA data products described above.Our base method is to use mean 48 hr 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, October 2015, and January, April 2016.Other studies (e.g., Janardanan et al., 2016;Fischer et al., 2017) 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 10 without the additional turbulence term.
Figure 2 shows back trajectories for one layer and 2 different times that end at the observation sites.Trajectories from multiple vertical levels are combined to determine residence times or footprints as described in Appendix C. HYSPLIT shows 3 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 Central Valley (northwest, see Fig. E1).

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 compared with 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 (e.g., Turner et al., 2016;Lauvaux et al., 2016); they can be solved by analytical or adjoint methods (Rodgers, 2000;Kopacz et al., 2009).Different cost functions can be used, which might change the results.Here we test 3 different methods.The first is a Kalman filter (described in Appendix D) which is computationally cheap, but has only one 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 3 Typical X gas enhancements Several previous studies have discussed the SoCAB X CO2 , X CH4 , and X CO enhancements from local anthropogenic activity (Wunch et al., 2009;Kort et al., 2012;Janardanan et al., 2016;Wunch et al., 2016;Hedelius et al., 2017a;Schwandner et al., 2017).There have also been several studies which have discussed enhancements noted from the 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 2 and CH 4 (Wong et al., 2015(Wong et al., , 2016)).For comparability we exclude enhancements from CLARS and in situ observations (e.g., Verhulst et al., 2017) in this section.Kort et al. (2012) noted that observing changes in typical X gas enhancements from space-borne instruments can provide a first order estimate of how local emissions have changed year-to-year.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.Changes in X gas enhancements can provide a first-order estimate of how much local emissions have decreased without the need for a full inversion.
Table 2 lists X CO2 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 CO2 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.In Fig. 3 are histograms of enhancements for all dates of this study.Most enhancements are on order of 2-3 ppm except for those from the recently published paper by Schwandner et al. (2017), which are about double.Though their enhancements are within the range of ∆X CO2 enhancements in the v7r and v8r histograms in Fig. 3 (bottom row), 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 2 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.
Models that assimilate only global in situ (i.e., no total column) CO 2 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 ppm or better to provide new information.The TCCON-predicted bias uncertainty is 0.4 ppm or less (<0.1 %).A long-term CO 2 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 on order of 4-9 years worth of emission reductions.

Carbon dioxide
Our flux estimate of CO 2 using the TCCON sites and linear model (Eq.E3) is 139±35 TgCO 2 yr −1 .An error assessment is described in Sec.4.3.This estimate is shown, along with estimates from past studies, in Fig. 4. Our estimate is comparable to Vulcan (Brioude et al., 2013;Fischer et al., 2017)  Histograms are for all dates of this study (Section 2.1).
v4.0 (as reported by 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 larger than 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 to rural areas.Blooming effects  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.EDGAR v4.0 value was reported by Wunch et al. (2009).All other values were found in a literature review.
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.
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 lower at 118±29 TgCO 2 yr −1 .This value varies by up to 43 TgCO 2 yr −1 depending on filtering methods (e.g., warn levels, WL).Warn levels are a global metric of data quality, where 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.After excluding just the 2 largest WLs, the total net flux varies by only 23 TgCO 2 yr −1 when applying additional filters.

CH 4 and CO
Using the same methodology we estimate a CH 4 flux of 325±81 Gg CH 4 yr −1 .This is less than the estimate by Wunch et al. (2009), but similar to estimates from CARB (Fig. 6).CARB-based CH 4 fluxes for just the SoCAB were estimated For every 5th sounding the set of backtrajectories is shown.
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 lower than previous estimates of CH 4 fluxes using in situ (tower and aircraft) data (Hsu et al., 2010;Wennberg et al., 2012;Peischl et al., 2013;Wecht et al., 2014;Cui et al., 2015).
We also estimate a CO flux of 555±136 Gg CO yr −1 .This is significantly less than the estimates by Wunch et al. (2009)  When their results are scaled down based on our posterior CO 2 fluxes (118-139 TgCO 2 yr −1 ), the CO flux is 830-970 Gg CO yr −1 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 −1 for 2015.

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 where 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 8 Tg CO 2 yr −1 (2σ), or about 6 %.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 E3), and how these are chosen can affect the final flux calculated.We evaluate 4 sensitivity tests (Fig. 7).For the first test, we filter out data where the observations differ from the Tracer-tracer relations were used for the TCCON (Wunch et al., 2009(Wunch et al., , 2016)), CLARS (Wong et al., 2015(Wong et al., , 2016)), and in situ observations (Hsu et al., 2010;Wennberg et al., 2012).The "W09 adj" are the results of Wunch et al. (2009) when adjusted for our posterior CO2 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).
model above a threshold.We start with a threshold that is a factor of 64 and adjust from there.We also adjust values of x a , S , and S a by factors of 2 −10 to 2 10 .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 10 % 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 χ 2 and the interannual range.We estimate an overall uncertainty of 10 % from these parameters.Hedelius et al. (2017b) reported 2σ measurement bias of less than ∼0.2 ppm X CO2 (central estimate, maximum range <0.5 ppm) between the AFRC and Caltech TCCON sites, but even a bias of 0.2-0.3ppm X CO2 will produce an error of ∼ 10 % in the flux.This bias could also arise from improper boundary conditions or application of averaging kernels.
We test the sensitivity to different inversion and modeling schemes (  from high model:measured values having unreasonably high weights in these particular schemes with few scaling factors (Appendix D2).GDAS and HRRR results are within uncertainty.
In summary, we estimate a 20 % uncertainty from model winds, 10 % from our choice of initial values, 10 % from observations and the boundary condition, 5 % from the prior flux (based on results of Lauvaux et al., 2016), and 5 % from additional For a given gas, all the inversions use the same observed ∆Xgas (Caltech TCCON − AFRC TCCON).The top 3 rows are from using different meteorological models, with the same inversion scheme (Eq.E3).The last 3 rows are from using the same meteorological model (NAM 12 km) with different inversion schemes.Errors are 25 %. a For GDAS Sa is 30× smaller random uncertainty.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 pseudoobservations (Turner et al., 2016).Both of these studies benefited from additional sites (9 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.

Emission Ratios
Emission ratios can help us evaluate the inversion for the SoCAB.Previous studies (Newman et al., 2016;Wunch et al., 2009Wunch et al., , 2016) ) have 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.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.Interannual ranges for CO, CH 4 , and CO 2 are 19 %, 13 %, and 11 % compared to their central estimates.The interannual range of ratios for CO:CH 4 , and CO:CO 2 are similar at 21 % and 22 % respectively, but the CH 4 :CO 2 range is smaller at 2 %.
We estimate emission ratios using the solar zenith angle (SZA) anomaly method described by Wunch et al. (2009Wunch et al. ( , 2016)), as well as from 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).We estimate the emission ratio from the work of Verhulst et al. (2017) using the weighted mean of the excess ratios from their 5 in-basin sites, with weights 1 σ 2 .Emission ratios from the SZA anomaly method and the differenced enhancement are in agreement with ratios from previous studies (Fig. 8).The CH 4 :CO 2 from the inversion is slightly lower than but similar to other studies.The CO:CO 2 ratio is also lower.Based on the CARB inventories, a decrease is expected because CO emissions have decreased more than CO 2 emissions over the past decade.
In November 2015, the large CH 4 :CO 2 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 2 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 2 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 2 during June-August 2016 (https://firetracker.scpr.org/wildfires/archives/).

Weekend effect
The weekday to weekend flux ratios (WD:WE) are listed in Table 4.The uncertainty is estimated to be ±0.17 based on changes 10 in the ratio from the S a scaling test up to 8× (Fig. 7).Weekday:weekend ratios are larger than those from previous studies for CO (Pollack et al., 2012;Brioude et al., 2013).Compared with the prior, the CO 2 ratio is scaled down, and the CO ratio remains equal to the prior.Methane has a ratio that is larger than unity.In contrast to our inversion results, methane is not expected This study 1.11 ± 0.17 1.23 ± 0.17 1.17 ± 0.17 a WD:WE CO ratios from Pollack et al. (2012) were calculated using the difference between the CalNex-Pasadena and Mt.Wilson Flasks (Table 2 therein).For CO2 we used the CO WD:WE ratios with the CO/CO2 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 to vary significantly on weekdays compared to weekends because production from biogenic sources and fugitive losses from natural gas infrastructure are less time variant than CO and CO 2 emissions.

Conclusions
This study demonstrates a method to readily obtain estimates of net CO 2 fluxes over regions on order of 10,000 km using only remote sensing observations.This method could be applied almost anywhere globally using only OCO-2 or other space-based observations of CO 2 (e.g., GOSAT) without the need for ground observations, or a specialized model.Our estimates of total annual CO 2 fluxes from the SoCAB using HYSPLIT with NAM 12 km as our dynamical model are similar to some previous estimates (Fig. 4), but less than inventory values reported in tracer-tracer flux estimate papers (Wunch et al., 2009;Wong et al., 2015).This has important implications for these studies, which would have overestimated CH 4 emissions if CO 2 emissions were also too large.Net CO and CH 4 fluxes are slightly less than previous studies.
The overall uncertainty is 25 %, with the dynamical model contributing the most.We consider an uncertainty of 25 % to be large and 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 (e.g., Kort et al., 2012).The wide range of uncertainty suggests that CO 2 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 aid in not only improving the accuracy of the overall flux, but also may be incorporated into spatiotemporal inversions to map fluxes from sub- in quadrature of the error from each site, i.e., 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 calculated by where n is the number of measurements, ẑi are the individual X CO2 measurements, ẑi,err are the reported errors associated with the measurements, and ȳ is the weighted average using ẑ−2 err as weights.Note that Eq.D2 takes into account both the measurement errors as well as the spread of the measurements.

D1 Iterations
We initialize the iterations with an arbitrary scaling factor α 0 = 1 and an associated error of σ α0 = 0.7.These initial values We iterate over the k measurements by calculating the partial derivative: 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 (Eq.A4 and A5 in Kleiman and Prinn (2000)): 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 and the state estimate is D2 A note on single scale factor inversions with large outliers Some single scale factor inversions can be written in the form where "mod" represents the initial model values.We consider the case when cost function is of the form: This is demonstrated in a sensitivity test, where 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 (1.1) n s where 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 non-linear model because it takes the longest.These results show that having a few large outliers in the Kalman filter (1 scale factor) and the non-linear (40 scale factor, Sect.E) inversions can significantly pull the results compared with the linear (∼1,000 scale factor, Sect.E) inversion.

Appendix E: Bayesian inversions
The Bayesian approach to solving atmospheric inverse problems has been described in more detail by Rodgers (2000) (see Section 2.3.2).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 (−2 ln P (y|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).

E1 Forward model
The generalized forward model can be written as where is an error term.We test 2 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 2 measurement locations.This model has 40 state vector scaling factors λ in 7 different classes corresponding to year (6), month (12), weekday/end (2), time of day (6), vertical level (5), spatial bin (8), and overall (1).Time of day bins cover 4 hours each with local ending times at 3:00, 7: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 agl.These are designed to help diagnose transport or footprint extent errors, and the upper 2 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 (E2) Here, m = Σt×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 0.01 • ×0.01 • 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 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., from 2012 and 2017) and not used when reporting our total fluxes.We remove elements which are not linearly independent which reduces the actual number used to less than (about one-fifth) the number of

Figure 1 .
Figure 1.A priori flux maps for CO2 (left) and CH4 (right) for select months.The same spatiotemporal prior for CO2 (ODIAC2016) was used for CO, but scaled to 1 % on a per mole basis.The ODIAC product was downscaled to 0.01 • resolution.The methane prior was created based on point sources, total emissions, and the population distribution.

Figure 2 .
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 2 different TCCON receptor sites at 14:00 (UTC-7).Magenta trajectories end at Caltech.Cyan trajectories end at AFRC.

10
to the observations (Appendix D2).We also use Bayesian inversions based on the methods ofRodgers (2000) (described in Appendix E).One Bayesian inversion is based on a non-linear forward model with 40 different scaling factors (Eq.E2), 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 Atmos.Chem.Phys.Discuss., https://doi.org/10.5194/acp-2018-517Manuscript under review for journal Atmos.Chem.Phys.Discussion started: 31 May 2018 c Author(s) 2018.CC BY 4.0 License.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 methodsIn summary, we have 4 sets of observations of X gas differences: Caltech TCCON − AFRC TCCON (CO 2 , CH 4 , and CO), and OCO-2 − AFRC TCCON (CO 2 ).We use one gridded spatiotemporal inventory for CO 2 and CO (ODIAC2016, with a weekly pattern for hourly emissions), and one gridded spatiotemporal inventory for CH 4 (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 3km 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 non-linear 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.

Figure 3 .
Figure 3. Histograms of Xgas enhancements observed in the SoCAB.Data are averaged for ± 30 minutes centered on the hour.Top row: 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 getting transported from downtown Los Angeles (southwest) to Caltech, and from mixed layer dynamics.Bottom row: Enhancements are OCO-2 observations minus AFRC TCCON observations.Colors represent the distance from the Caltech TCCON site.

Figure 4 .
Figure 4. Estimates of SoCAB CO2 fluxes (annual estimates from TCCON are shown as black triangles and OCO-2 are shown in pink) compared with previous studies.Top-down (TD) flight estimates are from Brioude et al. (2013).TD estimate from Ye et al. (2017) is based on

Figure 5 .
Figure 5.A visualization of OCO-2 observations and the forward model used in the flux inversion on June 20, 2015.The nadir track is shown starting at the bottom and ∼ 117.6 • W and going towards the northwest.Observations are overlaid on the ODIAC prior at 14:00 (UTC-7).

Figure 6 .
Figure 6.SoCAB CH4 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.

Figure 7 .
Figure 7. Assessment of sensitivity to initial values.On the left is the reduced state vector with 7 categories (overall, spatial, vertical, time of day, weekday/end, 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 Sa test.As Sa gets larger, the variability in the retrieved γ factors increases.The right shows overall fit parameters, including the overall flux, the dofs, χ 2 , Pearson's correlation coefficient between observed and postinversion model values, and the interannual range.Note the log2 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 Sa or Se are increased, or xa is scaled down.Generally the total flux is unchanged except for scaling xa which increases the flux by about 10 % of the change in the prior.Here the goal was to simultaneously increase dofs, decrease χ 2 , and increase r while keeping the interannual variability below about 25 %.

Figure 8 .
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 Pasadena (Caltech) to Lancaster (AFRC) gradients is based on the difference between the TCCON observations.Values from Silva et al. (2013); Silva and Arellano (2017) were part of global studies.

Figure C1 .
Figure C1.Maps of monthly averaged residence times in the ML per pixel for trajectories ending at 21 UTC, shown for all times leading up to the observation.Pixels are 0.01 • × 0.01 • , or approximately 1.03 km 2 .In July the origins were more predictable, but in January there was greater variation.

Figure D1 .
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 non-linear inversion are more affect by a few strong outliers than the linear inversion.

Figure E1 .
Figure E1.Extent of the 8 spatial sub-regions.C=center, Q1-Q4 are SoCAB quadrants, O=ocean, CV=Central Valley, D=all other areas, mostly the Mojave Desert to the northeast.
m j,CIT − m j,AFRC ) .(E3) tions.The 0.1 • methane inventory was produced by Harvard University in collaboration with the EPA.The authors gratefully acknowledge the NOAA Air Resources Laboratory (ARL) for the provision of the HYSPLIT transport model (http://www.ready.noaa.gov)as well as the gridded archived meteorological data used in this publication.Atmos.Chem.Phys.Discuss., https://doi.org/10.5194/acp-2018-517Manuscript under review for journal Atmos.Chem.Phys.Discussion started: 31 May 2018 c Author(s) 2018.CC BY 4.0 License.
a Qualitative estimate based on Fig.1and Supplemental Fig.3therein.b We modified the boundary condition compared to our previous work (see Appendix A), values are for 14:00 (UTC-7).

Table 3 .
Fluxes from various methods.

Table 3
Hedelius et al. (2017a)ovements for Modeling Emissions by Scaling (TIMES) were reported for the contiguous United States byNassar et al. (2013).dHestia-LA is based on Fig.2fromHedelius et al. (2017a).This same ratio is used in the CO and CO2 priors in this study.