A- A+
Alt. Display

# Bayesian inverse estimation of urban CO2 emissions: Results from a synthetic data simulation over Salt Lake City, UT

## Abstract

Top-down, data-driven models possess ample power to improve the accuracy of bottom-up carbon dioxide (CO2) emission inventories, and more work is needed to explore the merger of top-down and bottom-up estimates to better inform the metrics used to monitor global CO2 fluxes. Here we present a Bayesian inverse modeling framework over Salt Lake City, Utah, which utilizes available CO2 emission inventories to establish a synthetic data simulation aimed at exploring model uncertainties. Prescribing a high-resolution, urban-scale data product (Hestia) as the “true” emissions in the model, we combine prior emissions with an atmospheric transport model to derive modeled afternoon CO2 enhancements at six monitoring sites within the Salt Lake Valley during the month of September 2015. A global high-resolution gridded emissions data product (ODIAC) is used as the prior, and objective uncertainty structures are defined for both the a priori estimates and the transport model-data relationship which consider non-negligible spatial and temporal covariances. Optimized (posterior) emissions over the Salt Lake Valley agree closely with the assumed “true” emissions during afternoon times, while results including unconstrained times (e.g. night-time) lack such agreement. Both spatial and temporal correlations of prior errors were found to be necessary for obtaining a robust posterior estimate. Model sensitivity analyses are performed, which examine correlation length and time scales, model-data mismatch error, and measurement site network variability. Through these analyses, one measurement site is identified as being particularly prone to introducing bias into posterior emissions due to influences from a nearby point source. Increasing model-data mismatch error at this site is shown to reduce bias in the posterior without significantly compromising agreement with monthly averaged true emissions.

##### Knowledge Domain: Atmospheric Science
How to Cite: Kunik, L., Mallia, D.V., Gurney, K.R., Mendoza, D.L., Oda, T. and Lin, J.C., 2019. Bayesian inverse estimation of urban CO2 emissions: Results from a synthetic data simulation over Salt Lake City, UT. Elem Sci Anth, 7(1), p.36. DOI: http://doi.org/10.1525/elementa.375
Published on 05 Sep 2019
Accepted on 03 Aug 2019            Submitted on 29 Jan 2019
Domain Editor-in-Chief: Detlev Helmig; Institute of Alpine and Arctic Research, University of Colorado Boulder, US
Associate Editor: Lori Bruhwiler; Global Monitoring Division, NOAA Earth System Research Laboratory, US

## 1. Introduction

Anthropogenic greenhouse gases (GHGs), particularly carbon dioxide (CO2), are the driver of rising global mean surface temperature, a trend which climate models predict may surpass 3.5°C higher than pre-industrial levels by 2100 without additional mitigation efforts (IPCC, 2014). Limiting future global mean surface temperature rise to 1.5°C could mitigate significant risks of climate change such as sea level rise, ecosystem stress and food security; however, achieving this goal will require reduction of CO2 emissions to near net-zero within the next 30 years (IPCC, 2018). Cities account for ~70% of global emissions (IEA, 2008), and as urban populations grow by a projected 2.5 billion people by 2050 (UN-DESA, 2014), cities are expected to play a major role in reducing carbon emissions in coming decades.

Recent climate negotiations as part of the United Nations COP-24 meeting in Katowice, Poland have demonstrated global motivation to enact ambitious strategies to reduce GHG emissions (Plumer, 2018), leaving cities with significant responsibilities to meet these goals. Numerous cities across the world have already pledged to reduce GHG emissions, with programs such as the C40 Climate Leadership Group (C40 Cities, 2016), the Covenant of Mayors (Lombardi et al., 2014), and the U.S. Conference of Mayors (Wang, 2012) exemplifying cooperation across the globe. This is especially significant in the wake of the U.S. decision to withdraw from the Paris Climate Accord, in which U.S. sub-national entities such as cities and states are bearing extra responsibility in upholding emission reduction strategies (Tabuchi and Fountain, 2017; US Conf. Mayors, 2017).

Improved understanding of the urban carbon cycle is needed to independently evaluate emissions over time and to provide detailed guidance on potential urban emission mitigation opportunities. Mitigation and management techniques that focus on specific sectors or sub-city domains can help uncover opportunities for energy efficiency and emission reduction actions (Salon et al., 2010; Dhakal and Shrestha, 2010). “Bottom-up” emission estimates seek to achieve these goals via a collection of approaches such as direct flux monitoring, energy consumption statistics and activity modeling; however, only a limited number of data products exist for urban domains which comprehensively distinguish activity sectors at very fine resolutions (Gurney et al., 2012, 2017; Newman et al., 2016; Patarasuk et al., 2016; Gately and Hutyra, 2017; Bun et al., 2018). When bottom-up inventories are not available, researchers often need to rely on data products that downscale global/national estimates using spatial proxies (e.g. CDIAC, EDGAR, ODIAC). These global products are widely available, however, oftentimes the products are less resolved representations of local emissions behavior, which can misallocate emissions in time and space due to the spatial proxies used (Gurney et al., 2018; Gately and Hutyra, 2017).

Uncertainties unique to both bottom-up and global downscaling methods can be addressed and constrained with “top-down” (atmospheric) estimates, which infer emission patterns utilizing real-time measurements of concentrations and atmospheric transport models. Yet, top-down methods are prone to uncertainties that exist within atmospheric transport models (Peylin et al., 2011; Díaz Isaac et al., 2014), along with measurement representation errors and inconsistent data availability (Gerbig et al., 2009; Lauvaux et al., 2012; Turner et al., 2016). To merge bottom-up and top-down methods, Bayesian inverse frameworks can provide a quantitative and rigorous approach for flux estimation by pairing atmospheric observational constraints with spatial and temporal structures describing prior estimates and uncertainties in both fluxes and model-data relationships (Enting et al., 2002).

Flux inversion methodologies have long been used by the atmospheric science community (Tarantola, 1987; Enting and Mansbridge, 1989, 1991), aiding in flux estimation for trace gases using a combination of data and transport models. Recent studies have applied these frameworks to examine GHG behavior over cities, providing top-down optimized emissions estimates of CO2 and methane (Sargent et al., 2018; Miles et al., 2017; McKain et al., 2012, 2015). In slight contrast to these frameworks, Bayesian-style approaches incorporate uncertainties of observations and fluxes as constraining inputs, and have been successfully applied over urban ground-based measurement networks to optimize existing bottom-up estimates of CO2 fluxes over Cape Town (Nickless et al., 2018), Indianapolis (Lauvaux et al., 2016; Oda et al., 2017; Gurney et al., 2017; Turnbull et al., 2019), Paris (Bréon et al., 2015), and Davos (Lauvaux et al., 2013). Studies have also used aircraft campaign data as the atmospheric measurements for urban inverse analyses over Los Angeles (Gourdji et al., 2018; Cui et al., 2015; Brioude et al., 2013) and Houston (Brioude et al., 2011, 2012), while satellite and ground-based data were used to analyze carbon fluxes over California (Fisher et al., 2017, Hedelius et al., 2018).

Given that the accuracy of inversion estimates is often not adequate to verify reported GHG emissions (National Research Council (NRC), 2010), recent studies have emphasized the exploration of underlying uncertainty covariance structures, which drive Bayesian inverse estimates. Within these studies, specific attention has been paid to the correlation of prior emissions uncertainties in either the spatial (Wu et al., 2018; Lauvaux et al., 2016; Oda et al., 2017) or temporal dimension (Nickless et al., 2018; Breon et al., 2015). However, the combination of space and time has yet to be formally examined within an urban anthropogenic CO2 inversion context. Similarly, correlation of errors within the observation-model transport relationship (i.e. “model-data mismatch”) has seen limited experimentation using non-negligible off-diagonal correlation terms (Lauvaux et al., 2016; Oda et al., 2017). In this study, we aim to further explore these structures by incorporating variations in both prior error and model-data mismatch covariance matrices within an atmospheric inversion over Salt Lake City, Utah.

Salt Lake City lies at the northern end of the Salt Lake Valley, geographically positioned between the Wasatch and Oquirrh mountain ranges in northern Utah. Encompassing Salt Lake City and its surrounding metropolitan area, the Salt Lake Valley (SLV) region is a rapidly developing urban area with a population of over 1 million people. Within the SLV, a number of stationary atmospheric monitoring sites and mobile platforms constitute an extensive, long-standing network of CO2 measurements (Lin et al., 2018; Mitchell et al., 2018a, 2018b). To date, the sole inversion study using this network was by McKain et al. (2012), which matched modeled anthropogenic CO2 enhancements with observed signals at 5 sites around the valley. However, no follow-up studies have taken advantage of the SLV CO2 network for inverse analysis of carbon emissions. Because the Salt Lake Valley is one of a handful of urban regions (e.g. Indianapolis and Los Angeles) where extensive measurement networks have aligned with research interests to create a high-resolution emissions data product (Patarasuk et al., 2016; Gurney et al., 2012, 2017; Newman et al., 2016), a SLV Bayesian inversion analysis can contribute towards a growing pool of urban studies that seek to better understand the urban carbon cycle.

As a preliminary assessment of the SLV inversion framework, we implement a “synthetic data” experiment in which known urban emissions are used to generate synthetic CO2 observations, which are then used in an atmospheric inversion to solve for emissions. Because the “true” emissions are known, synthetic data experiments allow researchers to evaluate the impact of different data/flux parameters on possible deviations from the known answer. These studies, often referred to as observing-system simulation experiments (OSSEs), have been used by numerous researchers to explore uncertainties in various inversion frameworks (Wu et al., 2018; Turner et al., 2016; Kort et al., 2013; Wu et al., 2011; Gerbig et al., 2006). Within the OSSE framework, no boundary in-flow or biogenic enhancements are considered, and transport is consistent within the model and synthetic data; thus, biases from these model components can be ignored, allowing us to better examine the underlying uncertainty structures in the inversions and their effects on model performance. In this OSSE study, we define a baseline configuration for model parameters and incorporate spatially- and temporally-resolved prior uncertainty, along with spatial and temporal correlations of prior errors and component-level derivation of model-data mismatch error structures. Thereafter, we perform sensitivity tests to assess this baseline setup and examine covariance within prior uncertainty and model data mismatch error, as well as the relative influence of network measurement sites used in our study.

## 2. Methods

Our domain of interest follows the spatial bounds of Salt Lake County, which occupies the Salt Lake Valley and parts of surrounding mountains to the East and West (Figure 1). Salt Lake County encompasses six measurement sites within the SLV CO2 observational network which are used in this study. Sites are shown in the domain map in Figure 1 and are labelled with their respective site codes – UoU (University of Utah), DBK (Daybreak), RPK (Rose Park), SUG (Sugarhouse), MUR (Murray), and SUN (Suncrest) (Lin et al., 2018; Mitchell et al., 2018a). As the focus of our inverse method, we solve for CO2 emissions at a fine-scale spatial resolution of 0.01° × 0.01°, resulting in a total of 2386 grid cells within our Salt Lake County domain. We adopt a temporal resolution of 6-hourly time-steps corresponding to 0–6; 6–12; 12–18; and 18–24 UTC (local time is Mountain Daylight Time (MDT), equivalent to UTC-6 hours). Emissions are optimized starting from September 1, 18:00 UTC through September 30, 6:00 UTC, with a total of 114 time-steps which span over a roughly 4-week window.

Figure 1

Map of the Salt Lake Valley domain. The study domain in northern Utah is shown as the red box, with Salt Lake County represented by the unshaded region on the right. Yellow circles represent CO2 monitoring locations across Salt Lake County used in this study. DOI: https://doi.org/10.1525/elementa.375.f1

### 2.1. Emission Inventories

Two independent CO2 emission inventories are used as state vectors in the inversion (Section 2.3), one acting as our prior estimate (primarily a global downscaling method), and the other as the set of “true” emissions (a bottom-up approach) which we aim to recover through optimization of the prior.

We use the 2016 version of the Open-source Data Inventory for Anthropogenic CO2 (ODIAC2016, Oda and Maksyutov, 2015) as the a priori input estimate for our inversion. ODIAC is a global high-resolution (1 km) monthly fossil fuel CO2 emission data product which is designed and developed primarily for global and regional tracer transport and inversion applications, but also has been used for several urban studies (e.g. Lauvaux et al., 2016; Oda et al., 2017). ODIAC is based on downscaling of bottom-up CO2 emissions estimates using power plant emissions and geolocation information taken from the Carbon Monitoring and Action (CARMA) database (www.carma.org), as well as satellite-observed nightlight data (Oda and Maksyutov, 2011; Oda et al., 2018). The native temporal and spatial resolutions of ODIAC are monthly and Arc 30 second respectively. Here, we use a temporal-downscaling and re-gridding method described by Nassar et al. (2013) to scale the emissions to a 6-hourly temporal resolution and 0.01° grid resolution. CARMA point source geolocation used in ODIAC places various point sources at locations which are misaligned with their counterparts in Hestia after re-gridding ODIAC to 0.01°. These differences are described in further detail in Gurney et al., (2018). Because inversion frameworks generally cannot correct these gridding discrepancies on their own, misaligned point sources are manually re-located within our prior state vector to match locations in Hestia, and the re-located grid cells are replaced with a mean inventory value solely for the purpose of this study. A map of afternoon-averaged (18–24 UTC) ODIAC emissions with re-aligned power plants is shown in Figure 2a for the month of September 2015.

Figure 2

Prior (ODIAC), “True” (Hestia), Posterior, and residual emission maps. CO2 emissions over Salt Lake Valley are averaged for afternoon time-steps over the month of September 2015 (measurement sites denoted with black points). The 4 panels show: (a) Prior emissions (ODIAC), (b) True emissions (Hestia), (c) posterior (optimized) emissions and (d) emission corrections (posterior minus prior). Emissions are capped at ±15 μmol m–2 s–1 for visualization. Posterior emissions retain prior spatial structure but adapt point sources and faint on-road emission structure from truth (Hestia) visible in panels (c) and (d). DOI: https://doi.org/10.1525/elementa.375.f2

The Hestia-SLC v2018.01.24 fossil fuel CO2 emissions data product was defined as the “true” emissions within our synthetic data framework, due to its street-level granularity. The Hestia project consists of fine-scale mappings of CO2 emissions estimates using a bottom-up approach for select cities across the U.S. (including Salt Lake City). Hestia-SLC is comprised of eight different source sectors in Salt Lake County – Commercial (buildings and point), Industrial, On-road, Non-road, Railroad, Airport, Residential, and Electricity-production (Patarasuk et al., 2016; Gurney et al., 2012). Native grid spacing in Hestia is comprised of point, line and polygon sources, and are gridded for atmospheric modeling to 0.002° (Figure 2b). For the purposes of this study, emissions were re-gridded to 0.01° resolution, and the native hourly temporal resolution is used. Afternoon averaged Hestia emissions can be seen in Figure 2b over the same time frame as Figure 2a.

### 2.2. Model Transport

This study uses the Weather Research and Forecasting (WRF) model coupled with the Stochastic Time-Inverted Linear Transport (STILT) model to simulate atmospheric transport for this inversion (Lin et al., 2003; Nehrkorn et al., 2010). The domain setup and parameterization configuration used here is similar to that of Mallia et al., (2015), with the exception of the simulation time, which is set for September of 2015. WRF winds are used to drive an ensemble of 1000 STILT backward trajectories, which are released and traced backwards in time for 24 hours from each of the six site locations for all hours between September 3rd and 30th for a total of 4032 simulations. STILT has been widely used to interpret CO2 and other trace gases, including studies focusing on inverse modeling methods (Lin et al., 2003; Lin and Gerbig, 2005; Mallia et al., 2015; Lin et al., 2017; Mallia et al., 2018; Fasoli et al., 2018). Turbulent dispersion within STILT is parameterized as stochastic motions within the backward trajectories. These trajectories are then used to map the surface influence or “footprint” for each measurement time and location, quantifying the sensitivity of the observation to upwind source regions using units of concentration per unit flux (Lin et al., 2003).

Model transport acts as the linkage between gridded emissions and individual observations; therefore, footprint receptors are selected to match afternoon hourly-averaged observations from 12–17 local time (18–23 UTC), when the planetary boundary layer (PBL) is well-mixed within the valley at each of the six measurement sites used. Within “real-data” contexts, analysis is generally limited by large uncertainties in model transport outside of the afternoon, particularly from erroneous PBL height estimation. Here we adopt “real-data” conventions by excluding non-afternoon times from our observations and footprints. To encompass the first 24-hour period of emissions in the state vector, receptor times begin 24 hours after the beginning of the emissions vector. A total of 1008 footprints/observations are initially composed of 6 sites each, with 6 afternoon hours over 28 days (September 2, 18:00 UTC through September 30, 23:00 UTC). Data and footprints are further aggregated to represent daily afternoon averages in order to minimize the influence of model biases in PBL height on the variability of data between hours, following “real-data” conventions used in other urban studies (e.g. McKain et al., 2015; Miles et al., 2017). This method condenses data down to a total of n = 168 observations, using footprints averaged over all 6 hours. Figure 3 shows monthly average footprint mapping for each of the measurement sites over the month of September 2015.

Figure 3

WRF-STILT footprint (influence) averages. (a) WRF-STILT footprints are shown on a log10 scale across SLV domain (shaded region) including all observation locations. Footprint average includes afternoon observation times only (18:00–23:00 UTC) for September 2015. (b) Time series of percent influence within Salt Lake County domain vs. number of hours back from start of backward trajectories, averaged over all footprints. DOI: https://doi.org/10.1525/elementa.375.f3

### 2.3. Bayesian Inverse Methods

Posterior emission estimates are optimized through the minimization of the cost function described by (Tarantola, 1987) and (Enting, 2002):

(1)

Where z (n × 1) is a vector of observed enhancements, H (n × m) is a Jacobian matrix of footprint values which relate the measurements to the state vector of unknowns (s (m × 1)), R (n × n) is a square and symmetric matrix describing the covariance of model-data mismatch errors (also referred to as observational errors), sp is the state vector of prior emissions, and Q (m × m) is a square and symmetric matrix describing the covariance of deviations between the true field s and prior field sp. In the context of matrix dimensions, n is equal to the number of constraining observations, and m is equal to the number of unknown grid cells in the state vector. For this study, m is equal to mτ * ms where mτ = number of timesteps (114) and ms = number of grid cells (2386).

The posterior best estimate of emissions, ŝ, is given as the solution to the above cost-minimization function:

(2)

In addition, we construct the posterior uncertainty covariance matrix Vŝ, given as:

(3)

The solution to the cost function and resulting posterior uncertainty given above are solved following the computational methods described by Yadav & Michalak (2013), which express Q as a Kronecker product (explained in section 2.5) to solve for posterior fluxes at a significant cost reduction when the size of the state vectors (e.g. ŝ, sp, struth) become large. A minor variation of these methods is introduced when treating prior uncertainty as a spatiotemporally heterogeneous vector, for which we utilize methods described by NOAA’s CarbonTracker-Lagrange software (www.esrl.noaa.gov/gmd/ccgg/carbontracker-lagrange/). Further details on these computational methods are described in Text S-1 of the Supplemental Material.

### 2.4. Synthetic data framework

A key feature of this study is the use of synthetic data to drive corrections to the a priori estimate. The goal of this framework is to recover true emissions (Hestia) starting from the prior inventory (ODIAC). In order to achieve this, synthetic data are created which align with enhancements from the set of “true” emissions at all receptor locations and times, using footprints from the transport model to convert emissions to observed signals. Daily afternoon synthetic observations from 12–17 MDT are expressed as a single afternoon-aggregated observation for each site, where modeled signals from struth are obtained via convolution with daily afternoon-averaged footprints (described in the H matrix). Random noise is then added to the synthetic data, resembling possible “true” perturbations equal to the observational error described in the model-data mismatch (R) matrix in section 2.5. The equation for synthetic data generation is given as:

$\mathbf{\text{z}}=\mathbf{\text{H}}{\mathbf{\text{s}}}_{truth}+\epsilon$
(4)

where errors (ε) follow a Gaussian distribution with a mean of 0 and standard deviation equal to the standard error within the diagonal R matrix.

Because synthetic data are generated with an additional random component (ε), posterior estimates of emissions can vary due to random perturbations. To analyze inversion results, we obtain expected posterior values by running a Monte-Carlo-style ensemble of 10,000 inversion iterations and taking the mean posterior best estimate. Details of this averaging method are given in Text S-2 of the Supplemental Material.

### 2.5. Error Covariance Parameters

As mentioned earlier, the prior error covariance (Q) matrix describes both the variance in prior emissions uncertainty and the spatial and temporal correlation of these uncertainties. Here we derive Q from three distinct components: spatial covariance (E), temporal covariance (D), and prior error variance (σ2), where the spatial and temporal error covariance matrices are combined via a Kronecker product, described as in CarbonTracker-Lagrange documentation (www.esrl.noaa.gov/gmd/ccgg/carbontracker-lagrange/):

$\mathbf{\text{Q}}={\mathbf{\text{I}}}_{\sigma }\left(\mathbf{\text{D}}\otimes \mathbf{\text{E}}\right){\mathbf{\text{I}}}_{\sigma }$
(5)

Under this method, Iσ is a diagonal matrix whose elements describe the uncertainty of prior emissions, defined here as the magnitude of difference between the prior and true emission inventories, with a minimum value of 1 μmol m–2 s–1 assigned within the vector. The square of this uncertainty term (Iσ) in equation 5 comprises the prior error variance (σ2) component of Q. By this definition of Q, distinct uncertainty values are given to individual grid-cells for each time step in order to describe the heterogeneous nature of deviations between the prior and true emission inventories. Accurate representation of grid-level uncertainties is especially important in the real-data context, as these parameters are largely responsible for the spatially-explicit nature of flux inversion corrections. However, despite the challenging nature of assigning grid-level uncertainties in both global downscaling and bottom-up methods (Andres et al., 2016; Hogue et al., 2016; Hutchins et al., 2017; Gurney et al., 2018), within the context of our synthetic-data analysis we can assume Hestia gridded values to be “true” and therefore free from errors. Using this approach, we assume a prior uncertainty state vector for which each grid cell describes the deviation between its prior and “true” values. Thus, uncertainty is here equivalent to the 1σ range, reflecting the absolute value of gridded differences between ODIAC and Hestia.

For the remaining components of Q, spatial and temporal covariance matrices are defined using exponential decay equations as shown below:

$\mathbf{\text{E}}=\left[\mathit{\text{exp}}\left(-\frac{{\mathbf{\text{X}}}_{s}}{{l}_{s}}\right)\right]$
(6)
$\mathbf{\text{D}}=\left[\mathit{\text{exp}}\left(-\frac{{\mathbf{\text{X}}}_{\tau }}{{l}_{\tau }}\right)\right]$
(7)

These matrices are computed using separation distances (Xs) and lag-times (Xτ) between cells and timesteps respectively, divided by their corresponding correlation range parameters ls and lτ. Within the temporal correlation matrix, temporal correlation is assumed in fluxes across days only for time periods at equivalent hours of the day (non-matching times of day are considered uncorrelated and given a value of zero on the off-diagonal across the flux time domain). The correlation range parameters defined here describe the distance and time at which errors in the prior emissions are considered uncorrelated.

In order to determine the extent of spatial and temporal variability within prior errors, we employ the use of two objective methods for defining these range parameters. To find the spatial length-scale ls, we fit a variogram using a mapping of the average afternoon difference in the ODIAC and Hestia emission inventories over Salt Lake County, determining the spatial length-scale to be approximately 6 km. To find the temporal length-scale lτ, we analyzed the autocorrelation function of daily afternoon-averaged inventory differences over the month of September 2015, finding the temporal length-scale to be equal to ~2 days.

These covariance parameters, while subtle, control the extent to which corrections to prior emissions propagate or spread to neighboring grid cells and time steps. Despite particular attention given to these terms, known uncertainties exist due to downscaling methods used in the prior inventory (e.g., mis-specified point source geolocation in ODIAC). Because of this, we assess the validity of each correlation parameter value by exploring the effects of length scale variability on posterior emission results in section 3.2. Figures and further explanation of length scale determination can be found in Text S-3 of the Supplemental Material.

The model-data mismatch matrix (R) describes errors relating the transport model to the observations, and can be expressed as the sum of uncertainties given as:

$\begin{array}{l}\mathbf{\text{R}}={\mathbf{\text{R}}}_{\mathit{\text{part}}}+{\mathbf{\text{R}}}_{\mathit{\text{aggr}}}+{\mathbf{\text{R}}}_{\mathit{\text{transPBL}}}+{\mathbf{\text{R}}}_{\mathit{\text{transWIND}}}\\ \text{ }+{\mathbf{\text{R}}}_{\mathit{\text{eddy}}}+{\mathbf{\text{R}}}_{\mathit{\text{instr}}}+{\mathbf{\text{R}}}_{\mathit{\text{bkgd}}}+{\mathbf{\text{R}}}_{\mathit{\text{bio}}}\end{array}$
(8)

Individual components of the model-data errors are summarized in Table 1, and are described as follows:

Table 1

Description of components within the model-data mismatch matrix (R). All components are first given as hourly standard error (2nd column) and, if applicable, correlation is applied between same-day observations using a decay time-scale parameter (3rd column). Component errors are then aggregated to express observational error in enhancements that represent an average for a given afternoon (4th column). Standard error values (in ppm) are squared to express variance (in ppm2) in the R matrix. DOI: https://doi.org/10.1525/elementa.375.t1

Model-data mismatch component Hourly Standard Error Correlation Method Standard Error after afternoon aggregation

Rpart 0.1 ppm (fixed value) uncorrelated 0.041 ppm
Raggr 16.7% and 25% of mean enhancement (spatial and temporal aggregation) correlated with no decay in same afternoon 1.176 ppm
RtransPBL 7% of mean enhancement correlated with no decay in same afternoon 0.197 ppm
RtransWIND 35% of mean enhancement correlated with decay time-scale = 2.8 hours 0.881 ppm

Rpart is the error which stems from the release of a finite number of particles within trajectory ensembles. This value was determined following methods used in Mallia et al., (2015) but applied to the given emissions domain and found to be small (<0.1 ppm) given that a sufficiently large number of particles (1000) are released within STILT trajectory ensembles. Raggr is the error introduced in the model from aggregating spatially and temporally heterogeneous fluxes into single homogeneous cells and timesteps. Here we define this as the root mean square error (RMSE) of the difference between our chosen resolution (0.01 deg/6-hourly) and available finer resolutions (0.002 deg/1-hourly).

Transport model error in WRF-STILT is broken down into components describing horizontal wind error (RtransWIND), vertical mixing-layer height error (RtransPBL), and unresolved eddy turbulence error (Reddy). RtransWIND is estimated by comparing sets of WRF-STILT runs that included wind error, which are assumed to follow a Gaussian distribution and are thus unbiased. These wind errors are then incorporated as additional stochastic motions (Lin and Gerbig, 2005). The difference in variance between WRF-STILT simulated CO2 with and without wind errors are then used to estimate the impacts of transport errors.

RtransPBL is ideally constrained by comparing radiosonde measurements to WRF-modeled PBL heights; however, due to the lack of available radiosonde measurements in the afternoon over Salt Lake City, an approximated error of 7% mean enhancement was adopted from Gerbig et al., (2008) to represent errors using a high-resolution model for afternoon-only time periods. Finally, Reddy is accounted for within the above transport and aggregation errors and is neglected here.

Remaining errors in the equation (Rinstr, Rbkgd, and Rbio) do not pertain to our synthetic data approach, as physical instruments, boundary in-flow, and biogenic fluxes can be disregarded here, and are thus neglected in the baseline case.

Model-data mismatch errors in inversion contexts are often considered statistically independent, which is problematic in that these errors are often temporally correlated with other errors within a time window. Here, we attempt to characterize these correlations by assigning error correlations to the off-diagonal elements in Raggr, RtransPBL, and RtransWIND for hourly observations within each afternoon (but not across days or sites). Raggr and RtransPBL are considered to be fully correlated within each afternoon (as these error components are not expected to be random within a given afternoon), while RtransWIND is given decaying correlation with increasing time between observations, assuming a correlation timescale of 2.8 hours as in Mallia et al. (2017). In this study we assume error to be constant across sites, and do not account for spatial correlation of measurements between sites; however, it should be noted that in a real-data application, errors are likely to be correlated between towers, especially if separation distances between sites are small (as they are in this study area). Model-data mismatch errors are also likely to vary between sites based on a variety of factors such as instrument reliability, inlet height, local sources, and surrounding topography.

After hourly error correlations are applied, model-data mismatch errors are aggregated to the daily scale for each measurement site. Details covering this aggregation method are described further in Text S-3 of the Supplemental Material. The final standard error used here for observations across sites is 1.48 ppm, expressed on the diagonal of the R matrix as a variance of (1.48 ppm)2.

### 2.6. Sensitivity Analyses and Method Validation

In addition to our baseline case, we compare the results to a variety of sensitivity tests which evaluate parameters’ effects on model performance. In this analysis, we test the influence of spatial and temporal correlation range parameters (ls and lτ), model-data mismatch variance, and the site array. Values and configurations for each model parameter are varied to encompass a wide range of model scenarios, and a summary of the tested parameters is given in Table 2. Results for these sensitivity analyses are discussed in section 3.

Table 2

Sensitivity tests performed for this inversion analysis. Bolded values denoted with an asterisk (*) indicate baseline value from the inversion. Rdiagonal values here express standard error but are squared to express variance and are applied uniformly to the diagonal of R. DOI: https://doi.org/10.1525/elementa.375.t2

Parameter Values tested

ls 0, 1, 3, 5, 6*, 7, 8, 10, 20, 30 kilometers
lτ 0, 1, 2*, 3, 5, 14, 30 days
Rdiagonal 0, 0.5, 0.75, 1, 1.25, 1.5, 1.75, 2, 2.5, 3, 4, 5, 6, 7, 10 ppm
site network array S1* = UoU, RPK, DBK, MUR, SUG, SUN (all sites, baseline case)
S2 = RPK, DBK, MUR, SUG, SUN (UoU excluded)
S3 = UoU, DBK, MUR, SUG, SUN (RPK excluded)
S4 = UoU, RPK, MUR, SUG, SUN (DBK excluded)
S5 = UoU, RPK, DBK, SUG, SUN (MUR excluded)
S6 = UoU, RPK, DBK, MUR, SUN (SUG excluded)
S7 = UoU, RPK, DBK, MUR, SUG (SUN excluded)
S8 = RPK, DBK, MUR, SUN (UoU and SUG excluded)
S9 = RPK (all others excluded)

In order to evaluate the efficacy of the inversion model under certain parameters, we assess multiple measures of performance. Many possible measures exist which serve to validate atmospheric inverse models (see Michalak et al., 2017), and within this study we focus on a number of statistical validation methods that seek to verify the consistency of our assumptions within estimates of prior and posterior emissions and their respective uncertainties.

Our primary validation method compares posterior emission estimates with the “true” emissions on an aggregated space- and time-domain scale, averaged over a sufficiently large number of synthetic data iterations. To obtain these values, grid-scale emissions are first averaged in time, and then weighted by grid cell area to average in space. Similar methods are used to determine error reduction by calculating prior and posterior uncertainty aggregated in space and time, with prior and posterior covariances included within this aggregation technique. Through this, we obtain a measure of the amount of uncertainty that is reduced as a result of our constraints on emissions. Further description of these aggregation methods can be found in Text S-3 the Supplemental Material.

Modeled posterior observations, which are computed by multiplying footprints with posterior emissions, can be compared to synthetic data to evaluate model efficacy. Here we compute standard error (RMSE) and coefficient of determination (r2) between posterior and observed enhancements at each measurement site.

To assess the validity of posterior emissions given our prescribed error terms, we verify our results by calculating the reduced chi-squared value from our model residuals. Squared data and emissions residuals from our inversion are normalized by their respective variances in R and Q, and are expected to follow a chi-squared distribution with n degrees of freedom (among n + m residuals from data and emissions, respectively). A single reduced chi-squared value can be obtained by following the equation described by (Tarantola, 1987):

${\chi }_{r}^{2}=\frac{1}{\nu }\left[{\left(\mathbf{\text{z}}-\mathbf{\text{H}}\stackrel{^}{\mathbf{\text{s}}}\right)}^{T}{\mathbf{\text{R}}}^{-1}\left(\mathbf{\text{z}}-\mathbf{\text{H}}\stackrel{^}{\mathbf{\text{s}}}\right)+{\left(\stackrel{^}{s}-{\mathbf{\text{s}}}_{p}\right)}^{T}{\mathbf{\text{Q}}}^{-1}\left(\stackrel{^}{s}-{\mathbf{\text{s}}}_{p}\right)\right]$
(9)

where ν is the number of degrees of freedom, in this case equal to the number of observations in the inverse problem. A value of reduced chi-squared = 1 indicates alignment of residuals with their prescribed errors. Given that residuals are based on randomly-generated synthetic data, an ensemble of inversions must be generated in order to obtain an expected chi-squared value. Details of this calculation and computational limitations are discussed in Text S-1 of the Supplemental Material.

## 3. Results

We will first assess the agreement of posterior emissions, obtained using our baseline configuration, with Hestia, the reference set of true emissions. Further exploration of inverse model performance under varying conditions of spatial and temporal uncertainty, model-data mismatch, and observational network configuration are described in sections 3.2–3.4.

### 3.1. Model performance under baseline configuration

Figure 2c shows the posterior emissions (afternoon-averaged) over the Salt Lake Valley. Spatial patterns of optimized emissions still largely resemble those of the prior, with some additional point sources and on-road activity recovered from the true emissions grid (Figure 2b). Posterior corrections (posterior minus prior emissions, Figure 2d) are mostly positive in the northern downtown region of the domain, and largely control corrections to the prior at the domain-averaged scale. Differences between posterior and true emission maps can be seen in Figure S4 in the Supplemental Material. Figure 4 shows the time series for temporally resolved emissions averaged over the total SLV spatial domain, with only afternoon emissions shown (corresponding to times of day which are constrained with synthetic data). Posterior emissions are generally closer to truth than prior; however, emissions are overestimated from September 8–16 and even correct in the wrong direction, away from truth, from September 19–22 and 27–29.

Figure 4

Afternoon timeseries of prior, true, and posterior (optimized) domain-average emissions. CO2 emissions are aggregated over SLV domain for all afternoon timesteps in September 2015 (timestep intervals are 6-hourly). Prior emissions averages are shown in blue, true emissions in black, and posterior (optimized) in red. Prior and posterior emission uncertainties are displayed as light blue and orange shaded regions, respectively. September 19–22 and 27–29 are shown in dashed gray boxes to highlight days when posterior emissions correct downwards (in the opposite direction of the truth). DOI: https://doi.org/10.1525/elementa.375.f4

Due to inherent loss of information from atmospheric mixing, along with properties of both the prior baseline emissions values and uncertainty structures prescribed within the Q and R matrices, we do not expect the resulting corrections to exactly reproduce the truth at precise time and space resolutions. We therefore analyze the optimized results at the domain-aggregated monthly average scale, inferring a single value of net emissions over the Salt Lake County domain for September 2015. Using Monte Carlo methods described in section 2.4 and text S-2 of the Supplemental Material, the expected value of posterior emissions aggregated over the afternoon domain is 4.63 ± 0.03 μmol m–2 s–1, which compares favorably to the true emissions (4.63 μmol m–2 s–1) relative to the prior emissions (3.97 μmol m–2 s–1). The standard error of the average estimate, 0.03 μmol m–2 s–1, is less than 1 percent of true emissions. The corresponding afternoon domain-averaged posterior uncertainty reduction for the baseline case is 39.32%. Maps of prior uncertainty and uncertainty reduction (both in afternoon) are shown in Figure 5. It should be noted that within this synthetic data context, true emissions are known; thus, uncertainty reduction should be used solely as a comparison metric between scenarios of this inversion and is a statistical measure of the information content gained from the constraining inputs and uncertainties.

Figure 5

Prior uncertainty and uncertainty reduction maps. (a) Prior emissions uncertainty over SLV (afternoon average, September 2015), capped at 30 μmol m–2 s–1 for visualization. Uncertainty is defined as the absolute difference between ODIAC and Hestia at the 0.01° grid-cell level, with a minimum uncertainty of 1 μmol m–2 s–1 per cell. (b) Percent uncertainty reduction from prior after baseline inversion run (afternoon average, September 2015). Most uncertainty reduction occurs in the north of SLV where observation sites are concentrated. Observation sites are denoted by black diamonds in both figures. DOI: https://doi.org/10.1525/elementa.375.f5

Despite strong agreement in the afternoon at the spatially and temporally averaged scale, constraining influence from afternoon footprints diminishes when moving backwards in time throughout the day (Figure 3b). Consequently, morning hours which are intentionally unconstrained (i.e. 6–12 local time) retain some residual influence from afternoon footprints and are corrected to some degree. However, these corrections are based solely on the prior-true difference in the afternoon, giving overestimated emissions for the time period (ŝ = 3.63, sp = 3.54, and struth = 3.52 μmol m–2 s–1). Further back from the afternoon (i.e. 0–6 and 18–24 local time bins), no more than 10% of corrective influence remains within the domain, and posterior emissions revert to the prior. As a result, emissions averaged over all hours of day (ŝ = 3.64, sp = 3.45 and struth = 3.39 μmol m–2 s–1) also reflect differences in prior and true emissions from the constrained afternoon timeframe and are not as accurate a depiction of inversion performance as afternoon emissions.

Prior, observed and posterior enhancements at each site are shown in Figure 6. Similar to afternoon emissions, the prior enhancements (blue) often underestimate the observed “true” signal (black). Posterior corrections to the prior enhancements (red) align closely with the observed enhancements in general, and misalignments between observed and posterior signals are generally limited to instances at the RPK site when prior signals switch from underestimation to overestimation, and vice versa (i.e. when blue and black lines cross). This is likely due to the lag effect of the emissions corrections due to temporal correlation of posterior corrections.

Figure 6

Modeled observations (prior, posterior, and synthetic) at two sample sites. Derived afternoon CO2 enhancements are calculated using WRF-STILT footprints multiplied by prior, posterior, and (for synthetic data) “true” emissions. Results are shown for each measurement site with prior enhancements in blue, synthetic enhancements in black, and posterior in red. Synthetic and posterior observations are expressed as an average of an ensemble of 200 inversion runs, with the standard deviation of their observations shaded in gray and orange, respectively. September 19–2 and 27–29 are shown in dashed gray boxes to showcase abnormal behavior at the RPK site, when prior enhancements are higher than synthetic data values. DOI: https://doi.org/10.1525/elementa.375.f6

The expected value of reduced chi-squared value of fit for the baseline case is equal to 0.92. While this is slightly lower than the ideal value of 1.0, this represents relatively good fit of residuals to error matrices. As noted by Nickless et al. (2018), a reduced chi-squared value < 1 indicates that prior flux covariance and model-data mismatch errors may be overestimated. This could potentially be explained by the assignment of small modifications to the prior uncertainty structure (e.g. assigning a floor value of 1 μmol m–2 s–1).

### 3.2. Effects of spatial and temporal correlation on additive corrections

Results of sensitivity analyses on temporal and spatial correlation parameters are shown in Figures 7 and 8. Increasing correlation length, ls, as shown in Figure 7, results in additional spread in the spatial extent of uncertainty reduction. With small length-scales (e.g. ls = 0 km), isolated areas of correction (and therefore high uncertainty reduction) are present immediately surrounding the observational sites, despite lower domain-scale uncertainty reduction at these length-scales (Figure 8a). In contrast, larger length-scales show widespread reduction without such concentrated spikes around observational sites (Figure 7). Overall reduction in uncertainty increases when applying larger length scales, reflecting increased spread of corrective power to distant neighboring grid cells. The same effect is seen in the correlation time scale, where increasing time scale values show temporal spreading of corrections to neighboring days, strengthening the corrective trend (Figure 8b). As can be inferred here, mis-specification of spatial and temporal correlation scales may result in biased posterior emissions, where optimized emissions are overly influenced by strong trends in the difference between observations and prior-modeled enhancements.

Figure 7

Uncertainty reduction maps with varying spatial correlation length scale (ls). Maps are shown for ls = 0, 3, 10, and 30 km. Uncertainty reductions are averaged over September 2015. At small length scales (top panels) uncertainty reduction is confined to small areas of high reduction surrounding observation sites, whereas at larger length scales (bottom panels) uncertainty reduction is more widespread across the county domain. DOI: https://doi.org/10.1525/elementa.375.f7

Figure 8

Space- and time-averaged CO2 emissions vs. spatial and temporal correlation length scales (ls and lτ). Afternoon domain-averaged CO2 emissions (with shaded posterior uncertainty) are shown at (a) increasing spatial length scales and (b) increasing temporal length scales. Results from baseline configuration are marked by vertical lines. At small length scales for both parameters, posterior emissions underestimate the “truth”, while larger length scales show overestimation. DOI: https://doi.org/10.1525/elementa.375.f8

An additional inversion run was generated using a setup that neglects both spatial and temporal correlation in prior errors (ls = 0 km, lτ = 0 days). This configuration yields afternoon-averaged posterior emissions equal to 4.16 μmol m–2 s–1, which underestimates the true afternoon grid-averaged emissions by nearly 0.5 μmol m–2 s–1. It is thus apparent that given the prior (ODIAC) and true (Hestia) emissions used here, neglecting either spatial or temporal correlations in prior errors results in significant underestimation of optimized emissions.

### 3.3. Model-data mismatch error influence on inversion performance

Mis-specified model-data mismatch errors also have potential to introduce bias in posterior emissions estimates. Figure 9 shows comparisons of posterior emissions and uncertainty reduction to model-data mismatch error on the diagonal of the R matrix. Results show the expected behavior that with increasing uncertainty in the model-data relationship, correction to the prior diminishes. While this analysis does not consider any sort of systematic bias in the observations themselves, it asserts that underestimations of the model-data mismatch error result in misleadingly high reduction in uncertainty and overly-confident corrections to the prior (which are not necessarily correct) given observed enhancements with normally-distributed errors. Following this, overestimation of model-data mismatch error results in a posterior estimate that remains closer to the prior, limiting the model’s capability to provide full correction to match the true emissions.

Figure 9

Space- and time-averaged CO2 emissions and uncertainty reduction vs. model-data mismatch error. (a) Domain-average CO2 afternoon emissions and (b) uncertainty reduction percentages are plotted for scenarios using different values of model-data mismatch standard error (square root of the R matrix diagonal). Posterior emissions in (a) are shown in blue (with posterior uncertainty shaded) and are compared against prior (ODIAC) and “true” (Hestia) averages, dotted in black and violet, respectively. Posterior emissions begin to underestimate the truth as errors increase, approaching the prior average value at high error values. Uncertainty reduction in (b) is shown in green for afternoon-only timesteps (18–23 UTC) and in yellow for all timesteps combined (0–23 UTC). Uncertainty reduction is highest at ~60% in afternoon with no prescribed model-data mismatch error but approaches 0% as prescribed error increases. Results from baseline configuration are marked by vertical lines in both figures. DOI: https://doi.org/10.1525/elementa.375.f9

It should be noted that within this sensitivity analysis, model-data mismatch errors drive the magnitude of random perturbations (which follow a normal distribution) in the “true” signal, and thus analysis of these observational errors reflect increasingly varying synthetic observations. The calculated χr2 value of the baseline case (0.92) is less than 1, suggesting that prescribed uncertainties in R and Q may be slightly overestimated. However, as noted by Michalak et al. (2005), a χr2 = 1 metric is not in itself a comprehensive indicator of properly estimated covariance parameters. Additionally, “true” anthropogenic enhancements are no longer known in a real-data case (due to uncertainty in boundary conditions), leaving it up to the modeler to estimate any errors or systematic biases intrinsic to the model-data relationship. While not performed in this study, a potential additional method to quantify error covariance parameters in R, as well as Q, is via Restricted Maximum Likelihood estimation (Michalak et al., 2005). This method uses a top-down approach to optimize parameters in R and Q given a degree of prior knowledge about their underlying structures. Because the defining parameters and underlying structures of these covariance matrices are generally unknown and can have significant uncertainties, this method may act as a strong verification tool in future work.

### 3.4. Relative leverage of specific sites on emissions estimates

Within our synthetic data framework, we are able to assume 100% data availability at each observational site within our network. However, in the real world, routine maintenance of supporting instruments and hardware often prevents this from being the case. Thus, the final portion of this analysis focuses on the leverage that individual monitoring sites within our observational network have on the performance of the inversion. Using these results as guidance, future inversion studies over the SLV network may be better equipped to interpret the behavior of different sites’ data, and to identify time periods with maximally-informative data availability. Table 3 describes results of various iterations of the inversion performed with specific sites removed from the observational network. Sites whose removal results in a large change in domain-averaged emissions (e.g., UoU, SUG) contribute significantly to the corrective magnitude of the inversion, whereas sites whose removal results in significantly less uncertainty reduction (e.g. DBK) have more independent and unique upwind influence regions. It should be noted that if each site contributed equally and independently to the total uncertainty reduction (~39%), each site’s uncertainty reduction would be ~6.5%. This is, however, not the case, indicating that a degree of overlap exists among sites’ footprints. DBK comes closest at 6.17% and thus provides the largest degree of independent information of all sites, due to its isolated location and unique footprint over the western domain (Figure 3). Small changes to overall uncertainty reduction (<1%) indicate that a given site’s contribution to the information content of the inversion is less unique relative to other sites’ contributions. The UoU and SUG sites are examples of sites that exhibit low changes in uncertainty reduction; however, their combined removal results in a much larger loss of information to the inversion because their combined footprint covers a unique region within the geographical domain. Thus, neither site is individually imperative for the information network of the study, but their combined information content is needed for maximum performance of the inversion. While it is clear that all six measurement sites are needed to maximize the information content from an inversion over Salt Lake City, the relative contributions found here can serve as a guide to the utility of each site in a real data inversion context where data availability is less consistent.

Table 3

Inversion results with measurement sites (far-left column) excluded sequentially. Averages are computed for afternoon times only and are averaged over the entire Salt Lake County domain for the month of September 2015. Posterior uncertainty is expressed as Vŝ aggregated over afternoon times and full spatial domain. The change in uncertainty reduction for each configuration is relative to the baseline simulation (row 1), which has an original uncertainty reduction of 39.32%. DOI: https://doi.org/10.1525/elementa.375.t3

Excluded sites Posterior CO2 Emissions (Average) [μmol/(m2 s)] Posterior minus True CO2 Emissions (Average) [μmol/(m2 s)] ΔUncertainty Reduction [%]

None 4.64 0.00
UoU 4.55 –0.09 –0.90
DBK 4.50 –0.13 –6.17
RPK 4.69 0.05 –3.98
SUG 4.52 –0.12 –0.43
MUR 4.70 0.06 –3.51
SUN 4.62 –0.01 –2.27
UoU & SUG 4.24 –0.39 –3.88

#### 3.4.1. Bias introduced from proximity to point sources

A useful outcome of the network analysis is the ability to monitor unexpected behaviors in the optimized posterior with changing site network configurations. Here we examine more closely the time periods of September 19–22 and 27–29, when valley-averaged emissions appear to overshoot and are corrected downwards despite true emissions being higher than the prior (see Figure 4). In contrast to these periods, we compare with September 8–16, where the posterior corrections notably overestimate the truth on the valley-average.

To investigate this anomalous behavior specifically from the underestimation period, we note that exclusion of the RPK site in the network sensitivity analysis is the only configuration that does not see this period’s under-estimation biases. RPK site is unique among the six sites in being located within <5 km of four high-emitting point sources (according to Hestia), suggesting that a significant degree of hyper-local behavior may influence this site. As shown in Figure 6, large positive differences in observed vs. prior enhancements at RPK align with systematic afternoon over-correction in emissions from September 8–16. While the true signal is generally larger than prior at other northern sites, it is smaller than prior at RPK between September 19–22 and 27–29, aligning with underestimation of true emissions over the valley average (boxed portions of Figures 4 and 6).

To explore the linkage between these phenomena and uncover their driving behaviors, we examine the afternoon mean CO2 contribution maps for these two time periods at the RPK site (Figure 10), detailing the gridded contributions to true and prior CO2 enhancements (footprints multiplied by emissions) in upwind source regions. During September 8–16, the RPK site’s highest contributing grid cell from the true emissions is a point source located just south of the site. This same cell, however, exhibits almost no contributions from the true emissions during the September 19–22 and 27–29 periods. Prior and true enhancements are driven by the same transport, which infers that temporal differences in the proportional contributions must therefore be driven by temporal variability in footprints or the emissions themselves. However, because footprints averaged over each period do not display large differences in upwind source regions (see Figure S5 in the Supplemental Material), variability in emissions are the most likely cause. Indeed, the lower panel of Figure 10 shows that this particular source, identified as the Gadsby power plant, experiences large variations within the month. These variations align with the temporal patterns in RPK’s prior and true enhancements, as well as over- and under- estimations of posterior emissions from truth on the valley-average. The coincidence of these enhancements and domain-average biases enforces the idea that point sources can be a driving cause of biased behavior in the posterior emissions.

Figure 10

Signal contribution to the RPK site during two distinct periods, shown with emission time series from an influential nearby point source. Maps of average gridded contribution to prior and true CO2 signals at the RPK site are shown for September 8–16 (left column) and September 19–22, 27–29 (right column). A prominent point source is identified in the circled grid cell and the time series of emissions at the grid cell are shown below for the month of September 2015. Respective time periods are shown in the emission maps by dotted boxes in red (September 8–16) and violet (September 19–22, 27–29). DOI: https://doi.org/10.1525/elementa.375.f10

This naturally leads to the question: what can be done to mitigate the model bias at smaller temporal sub-scales introduced from a site like RPK? The answer hinges on whether the cause is due to misspecification of prior error covariance (Q), or of the model-data mismatch error (R). To explore this question, we examine results of additional model configurations to test if these biases can be addressed while preserving the baseline’s agreement with the truth. Two model configurations are run in addition to baseline and RPK-excluded setups: one with prior uncertainty of four nearby point sources reduced to 10% of original in the Q matrix (variation A), and one with 4 ppm of additional error introduced to RPK observations in the R matrix (variation B). Time series of spatially-averaged emissions from these variations are compared in Figure 11, showing that during the September 19–22 and 27–29 periods, the underestimations are actually amplified from baseline in the reduced sigma case. In contrast, these underestimations are minimized with additional error applied to RPK observations; afternoon domain-averaged emissions from variation B are equal to 4.64 μmol m–2 s–1, maintaining close alignment with the true afternoon average (4.63 μmol m–2 s–1).

Figure 11

Afternoon timeseries of prior and true domain-average emissions, compared to baseline and four additional model scenarios. CO2 emissions are in the same format as Figure 4, comparing baseline posterior (red) to additional scenarios. Scenario A is the same as baseline, but with uncertainty at four point sources surrounding RPK reduced to 10%. Scenarios B and C are the same as baseline, but with 2 ppm and 4 ppm RMSE, respectively, added to RPK observations on the R diagonal. The RPK-excluded scenario is same as baseline but with RPK removed from observations (only 5 sites used). RPK-excluded emissions do not exhibit negative corrective bias within the grey shaded columns; Scenarios B and C minimize this bias from baseline, and Scenario A shows exaggerated bias. DOI: https://doi.org/10.1525/elementa.375.f11

The reduced bias from increasing RPK observational error supports the notion that the biases introduced by the RPK site originate in the large enhancements from nearby point sources, rather than these sources’ prescribed uncertainty. Thus, despite true emissions being consistently higher than prior, the posterior is lower than prior during September 19–22 and 27–29 because of the corrections driven by these periods’ low signals from truth relative to prior at the RPK site.

Interestingly, RMSE between posterior and “true” enhancements is lowest at RPK of all sites (Table S-1 in the Supplemental Material), showing that this is not an indicator of bias introduced to the flux space. In this case, adjustment of the site’s error to “de-weight” observations proportional to other sites (basically attempting to increase RMSE between data and posterior signals) is one way to address concerns of hyperlocal enhancements disproportionally correcting emissions across the domain. The quantity of or extent to which error should be applied to this site (or similar sites within other domains) depends on the degree to which the bias introduced skews posterior emissions over the space and time domain in question. Thus, we are limited here by our chosen domain, as a single month is likely too limited a time scale to make accurate assertions about both the level of bias present in domain-averaged emissions and the appropriate error to assign to observations at this site.

Following the conclusions of Turner et al. (2016), a denser array of observational sites is likely to simultaneously reduce the bias in corrections introduced from the RPK site as well as minimize the information content lost from de-weighting these observations. This suggests that the SLV observational network could significantly benefit from additional receptors such as mobile measurements described in Mitchell et al. (2018b). Given the current functional network in this study with only six sites, we still seek to take optimal advantage of all available data, including those from the RPK site which ultimately contribute valuable information content to the inversion. Seeing as the RPK-excluded set-up is essentially a scenario which assigns infinite error to these observations, the appropriate additional error to minimize bias depends on the nature and magnitude of influence from surrounding point sources, as well as the density of remaining measurement sites.

### 3.5. Limitations

While the synthetic data inversion setup is useful in that it allows us to examine more closely the errors within our model, some limitations exist in the extent to which these findings can be fully applied to a real-data context. As mentioned earlier, the measures which are used here to justify the accuracy of optimized posterior emissions are limited to some degree by our selection of temporal domain. The time frame used here represents a month-long snapshot of emissions at 6-hourly resolution, over which variations in emissions could differ from other time periods at different scales. This month-long “snapshot” aspect of the chosen time domain limits our ability to quantify the magnitude of bias introduced by the interaction between the chosen site network and the true emissions field. Ideally, the observational network of an urban inverse study should be dense enough to capture all hyperlocal emissions activity that occurs within the domain; however, emission optimization using this Bayesian inverse framework is only as accurate as site network density and model transport allow. While baseline results presented here exactly match “true” emissions during the constrained period, this value is sensitive to the flux timespan and specific site configuration and could likely change if viewing over a different month or site network.

Another limitation of this study is the simplification of synthetic data, which are generated using model transport that is consistent with that used in the inversion. Realistically, data are prone to potentially large additional errors in biogenic CO2 signals (from flux estimates and transport) and other boundary conditions. These errors (along with horizontal and vertical transport errors) may be prone to systematic biases that are unknown and likely difficult to quantify. Correlation of observational errors both in time and space (i.e. between measurement sites) may be unknown or ambiguous as well, and within the fluxes themselves, a level of uncertainty exists in the methods used to determine spatial and temporal length scales. For these reasons, applying this framework to real-data contexts (especially over other urban domains) would necessitate additional considerations regarding unresolved uncertainties.

## 4. Conclusion and Discussion

In this study, we describe methods used to produce optimized emissions within a synthetic data framework, with emphasis on quantifying prior error covariance length scales and observational uncertainties. Expected afternoon averages from domain-averaged output closely match those of prescribed “true” emissions, with a standard error of about 0.03 μmol m–2 s–1 (less than 1% of the aggregated emissions). Domain-averaged uncertainties of afternoon emissions were reduced by around 39%. Emissions were constrained only during the afternoon (18–23:00 UTC), excluding other times due to concerns of large observational/transport uncertainties outside of these times. Our ability to correct for non-afternoon emissions given these constraints is limited based on these errors and results in a loss of around 50% information content when not isolating the constrained afternoon times from results. As previous synthesis reports and studies (e.g. NRC, 2010; McKain et al., 2012) have suggested, column-based constraints on emissions may be one way to strengthen the information content of the inversion, due to the decreased sensitivity of column-based signals to errors in the PBL height, which remains a significant modeling challenge outside the afternoon.

Proper estimation of correlation range parameters is proven to be important in this study in order to produce results in agreement with true emissions. We show here that neglecting spatial and/or temporal correlation, in the context of the given domain and prior/true emissions vectors, greatly limits the corrective power of the inversion model and results in poor agreement with true emissions. Spatial and temporal correlation of prior errors result in a degree of corrective spreading from areas with large influence/uncertainty into neighboring cells; however, when estimates of correlation range parameters are set too high, overcorrection can take place. Thus, objective determination of these parameters is highly recommended for optimizing model efficacy. One approach mentioned but not implemented in this study is Restricted Maximum Likelihood estimation, which could be a valuable tool to aid in determining covariance parameters for future studies.

Varying levels of model-data mismatch error are found to be a factor in introducing biases into the optimized emissions based on this synthetic data approach with normally-distributed random errors in observed CO2 enhancements. This study provides a comprehensive breakdown of total observational uncertainty which thoroughly assesses most components of error; deviating from this estimate of error was shown to compromise agreement with true emissions given the normally-distributed random perturbation method used in this study.

For future inverse analyses involving the Salt Lake Valley measurement network used here, we have shown that the DBK site contributes the most unique and independent information for constraining emissions. We have also noted the overlapping influence of the UoU and SUG sites. Overall, our network of sites displays sufficient spread in influence from upwind areas for the purpose of this analysis, but further work is needed to quantify gaps in information content resulting from unused locations within the urban domain where potential future sites could contribute to more effective inversion analyses.

Within this site network analysis, point sources are shown to have powerful influence on nearby sites, as is exemplified at the RPK site where large differences in prior and true signals propagated strong biased corrections across the emissions domain at certain times of the month. Increasing model-data mismatch error for RPK-specific observations is shown to reduce this bias without significantly compromising agreement with “true” emissions. While the spatial recoverability of missing point sources is shown to be adequate based on the setup used here, the location and magnitude of missing point source emissions are rarely known a priori in the real world. Thus, while recoverability is inherently limited by the resolution of emissions and transport models, as well as prior estimates of gridded uncertainty, further exploration of site network density and error covariance structures (e.g. using Maximum Likelihood Estimation or non-exponential spatial correlation to estimate Q) is needed to identify best practices for recovering missing point sources. Additionally, a multi-pronged approach to reducing bias from nearby point sources would also include the improvement of point source definitions within urban emission bottom-up data products. While global emission inventories do not specifically cater to high-resolution urban inversion analyses, adoption of more comprehensive point source allocation (which are sometimes available for individual cities or nations) in global inventories such as ODIAC will likely better suit them for future urban inversion applications.

## Data Accessibility Statement

Source code in R, sample input files, and select sample footprints are uploaded as online supporting information at https://github.com/lkunik/bayesian-osse-R-sample/. Updates and modifications to the code are ongoing, and contributions to the code framework are welcome and can be initiated by pull request. The version of software referenced for this study at the time of publication can be found here: https://doi.org/10.5281/zenodo.2655990. Questions about the code and files can be directed to the corresponding author.

## Supplemental files

The supplemental files for this article can be found as follows:

## Acknowledgements

We thank the Center of High Performance Computing at the University of Utah for providing the computational resources needed to carry out the modeling analysis. We would also like to acknowledge Arlyn Andrews, Logan Mitchell, Ryan Bares, and Ben Fasoli for their valuable input on modeling and observational aspects of this study.

## Funding information

This study was supported by NOAA Climate Program Office’s Atmospheric Chemistry, Carbon Cycle, and Climate program, award # NA17OAR4310084 (Utah).

## Competing interests

The authors have no competing interests to declare.

## Author contributions

LK and JCL designed the study. LK, DVM and JCL implemented the analysis. LK and DVM drafted the manuscript. DLM and KRG provided the Hestia emissions dataset, and TO provided the ODIAC emissions dataset used in this study. DVM provided the WRF-STILT footprints. All authors contributed valuable input to the final paper.

## References

1. Andres, RJ, Boden, TA and Higdon, DM. 2016. Gridded uncertainty in fossil fuel carbon dioxide emission maps, a CDIAC example. Atmos Chem Phys 16(23): 14979–14995. DOI: 10.5194/acp-16-14979-2016

2. Breón, FM, Broquet, G, Puygrenier, V, Chevallier, F, Xueref-Remy, I, Ramonet, M, Dieudonné, E, Lopez, M, Schmidt, M, Perrussel, O and Ciais, P. 2015. An attempt at estimating Paris area CO2 emissions from atmospheric concentration measurements. Atmos Chem Phys 15(4): 1707–1724. DOI: 10.5194/acp-15-1707-2015

3. Brioude, J, Angevine, WM, Ahmadov, R, Kim, SW, Evan, S, McKeen, SA, Hsie, EY, Frost, GJ, Neuman, JA, Pollack, IB, Peischl, J, Ryerson, TB, Holloway, J, Brown, SS, Nowak, JB, Roberts, JM, Wofsy, SC, Santoni, GW, Oda, T and Trainer, M. 2013. Top-down estimate of surface flux in the Los Angeles Basin using a mesoscale inverse modeling technique: Assessing anthropogenic emissions of CO, NOx and CO2 and their impacts. Atmos Chem Phys 13(7): 3661–3677. DOI: 10.5194/acp-13-3661-2013

4. Brioude, J, Kim, SW, Angevine, WM, Frost, GJ, Lee, SH, McKeen, SA, Trainer, M, Fehsenfeld, FC, Holloway, JS, Ryerson, TB, Williams, EJ, Petron, G and Fast, JD. 2011. Top-down estimate of anthropogenic emission inventories and their interannual variability in Houston using a mesoscale inverse modeling technique. J Geophys Res 116(20): 1–19. DOI: 10.1029/2011JD016215

5. Brioude, J, Petron, G, Frost, GJ, Ahmadov, R, Angevine, WM, Hsie, EY, Kim, SW, Lee, SH, McKeen, SA, Trainer, M, Fehsenfeld, FC, Holloway, JS, Peischl, J, Ryerson, TB and Gurney, KR. 2012. A new inversion method to calculate emission inventories without a prior at mesoscale: Application to the anthropogenic CO2 emission from Houston, Texas. J Geophys Res 117(5): 1–15. DOI: 10.1029/2011JD016918

6. Bun, R, Nahorski, Z, Horabik-Pyzel, J, Danylo, O, See, L, Charkovska, N, Topylko, P, Halushchak, M, Lesiv, M, Valakh, M and Kinakh, V. 2018. Development of a high-resolution spatial inventory of greenhouse gas emissions for Poland from stationary and mobile sources. Mitig Adapt Strat Gl , 1–28. DOI: 10.1007/s11027-018-9791-2

7. C40 Cities. 2016. C40 Cities Climate Leadership Group . [online] Available at: https://www.c40.org/about. [Accessed 30 August 2018].

8. Cui, YY, Brioude, J, McKeen, SA, Angevine, WM, Kim, SW, Gregory, J, Frost, GJ, Ahmadov, R, Peischl, J, Bousserez, N, Liu, Z, Ryerson, TB, Steve, C, Wofsy, SC, Gregory, W, Santoni, GW, Eric, A, Kort, EA, Marc, L, Fischer, ML and Trainer, M. 2015. Top-down estimate of methane emissions in California using a mesoscale inverse modeling technique: The South Coast Air Basin. J Geophys Res 120(13): 6698–6711. DOI: 10.1002/2014JD023002

9. Dhakal, S and Shrestha, RM. 2010. Bridging the research gaps for carbon emissions and their management in cities. Energy Policy 38(9): 4753–4755. DOI: 10.1016/j.enpol.2009.12.001

10. Enting, IG. 2002. Inverse Problems in Atmospheric Constituent Transport . Cambridge, UK: Cambridge University Press. DOI: 10.1017/CBO9780511535741

11. Enting, IG and Mansbridge, JV. 1989. Seasonal sources and sinks of atmospheric CO2 Direct inversion of filtered data. Tellus 41(2): 111–126. DOI: 10.3402/tellusb.v41i2.15056

12. Enting, IG and Mansbridge, JV. 1991. Latitudinal distribution of sources and sinks of CO2: results of an inversion study. Tellus B: Chemical and Physical Meteorology 43(2): 156–170. DOI: 10.1034/j.1600-0889.1991.00010.x

13. Fasoli, B, Lin, JC, Bowling, DR, Mitchell, L and Mendoza, D. 2018. Simulating atmospheric tracer concentrations for spatially distributed receptors: updates to the Stochastic Time-Inverted Lagrangian Transport model’s R interface (STILT-R version 2). Geosci Model Dev 11(7): 2813–2824. DOI: 10.5194/gmd-11-2813-2018

14. Fischer, ML, Parazoo, N, Brophy, K, Cui, X, Jeong, S, Liu, J, Keeling, R, Taylor, TE, Gurney, K, Oda, T and Graven, H. 2017. Simulating estimation of California fossil fuel and biosphere carbon dioxide exchanges combining in situ tower and satellite column observations. J Geophys Res Atmos 122(6): 3653–3671. DOI: 10.1002/2016JD025617

15. Gately, CK and Hutyra, LR. 2017. Large uncertainties in urban-scale carbon emissions. J Geophys Res Atmos 122(20): 11242–11260. DOI: 10.1002/2017JD027359

16. Gerbig, C, Dolman, AJ and Heimann, M. 2009. On observational and modelling strategies targeted at regional carbon exchange over continents. Biogeosciences 6(10): 1949–1959. DOI: 10.5194/bg-6-1949-2009

17. Gerbig, C, Körner, S and Lin, JC. 2008. Vertical mixing in atmospheric tracer transport models: error characterization and propagation. Atmos Chem Phys 8(3): 591–602. DOI: 10.5194/acp-8-591-2008

18. Gerbig, C, Lin, JC, Munger, JW and Wofsy, SC. 2006. What can tracer observations in the continental boundary layer tell us about surface-atmosphere fluxes? Atmos Chem Phys 6(2): 539–554. DOI: 10.5194/acp-6-539-2006

19. Gourdji, S, Yadav, V, Karion, A, Mueller, KL, Conley, S, Ryerson, T, Nehrkorn, T and Kort, EA. 2018. Reducing errors in aircraft atmospheric inversion estimates of point-source emissions: the Aliso Canyon natural gas leak as a natural tracer experiment. Environ Res Lett 13: 045003. DOI: 10.1088/1748-9326/aab049

20. Gurney, KR, Liang, J, O’Keeffe, D, Patarasuk, R, Hutchins, M, Huang, J, Rao, P and Song, Y. 2018. Comparison of Global Downscaled Versus Bottom-Up Fossil Fuel CO2 Emissions at the Urban Scale in Four US Urban Areas. J Geophys Res Atmos 124(5): 2823–2840. DOI: 10.1029/2018JD028859

21. Gurney, KR, Liang, J, Patarasuk, R, O’Keeffe, D, Huang, J, Hutchins, M, Lauvaux, T, Turnbull, JC and Shepson, PB. 2017. Reconciling the differences between a bottom-up and inverse-estimated FFCO2 emissions estimate in a large US urban area. Elem Sci Anth 5: 44. DOI: 10.1525/elementa.137

22. Gurney, KR, Razlivanov, I, Song, Y, Zhou, Y, Benes, B and Abdul-Massih, M. 2012. Quantification of fossil fuel CO2 emissions on the building/street scale for a large U.S. city. Environ Sci Technol 46(21): 12194–12202. DOI: 10.1021/es3011282

23. Hedelius, JK, Liu, J, Tomohiro Oda, T, Maksyutov, S, Roehl, CM, Iraci, LT, Podolske, JR, Hillyard, PW, Liang, J, Gurney, KR, Wunch, D and Wennberg, PO. 2018. Southern California megacity CO2, CH4, and CO flux estimates using ground- and space-based remote sensing and a Lagrangian model. Atmos Chem Phys 18(22): 16271–16291. DOI: 10.5194/acp-18-16271-2018

24. Hogue, S, Marland, E, Andres, RJ, Marland, G and Woodard, D. 2016. Uncertainty in gridded CO2 emissions estimates. Earth’s Future 4(5): 225–239. DOI: 10.1002/2015EF000343

25. Hutchins, MG, Colby, JD, Marland, G and Marland, E. 2017. A comparison of five high resolution spatially-explicit, fossil-fuel, carbon dioxide emission inventories for the United States. Mitigation and Adaptation Strategies for Global Change 22(6): 947–972. DOI: 10.1007/s11027-016-9709-9

26. International Energy Agency (IEA). 2008. World Energy Outlook 2008 . Paris: Organisation for Economic Co-operation and Development/International Energy Agency.

27. IPCC. 2014. Climate Change 2014: Synthesis Report. Contribution of Working Groups I, II and III to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change. In: Core Writing Team, Pachauri, RK and Meyer, LA (eds.), 151. Geneva, Switzerland: IPCC.

28. IPCC. 2018. Global warming of 1.5°C: An IPCC Special Report on the impacts of global warming of 1.5°C above pre-industrial levels and related global greenhouse gas emission pathways, in the context of strengthening the global response to the threat of climate change, sustainable development, and efforts to eradicate poverty. In: Masson-Delmotte, V, Zhai, P, Pörtner, HO, Roberts, D, Skea, J, Shukla, PR, Pirani, A, Moufouma-Okia, W, Péan, C, Pidcock, R, Connors, S, Matthews, JBR, Chen, Y, Zhou, X, Gomis, MI, Lonnoy, E, Maycock, T, Tignor, M and Waterfield, T (eds.). In Press.

29. Kort, EA, Angevine, WM, Duren, R and Miller, CE. 2013. Surface observations for monitoring urban fossil fuel CO2 emissions: Minimum site location requirements for the Los Angeles megacity. J Geophys Res Atmos 118(3): 1–8. DOI: 10.1002/jgrd.50135

30. Lauvaux, T, Miles, NL, Deng, A, Richardson, SJ, Cambaliza, MO, Davis, KJ, Gaudet, B, Gurney, KR, Huang, J, O’Keefe, D, Song, Y, Karion, A, Oda, T, Patarasuk, R, Razlivanov, I, Sarmiento, D, Shepson, P, Sweeney, C, Turnbull, J and Wu, K. 2016. High-resolution atmospheric inversion of urban CO2 emissions during the dormant season of the Indianapolis flux experiment (INFLUX). J Geophys Res Atmos 121(10): 5213–5236. DOI: 10.1002/2015JD024473

31. Lauvaux, T, Miles, NL, Richardson, SJ, Deng, A, Stauffer, DR, Davis, KJ, Jacobson, G, Rella, C, Calonder, GP and DeCola, PL. 2013. Urban emissions of CO2 from Davos, Switzerland: The first real-time monitoring system using an atmospheric inversion technique. J Applied Meteorol Clim 52(12): 2654–2668. DOI: 10.1175/JAMC-D-13-038.1

32. Lauvaux, T, Schuh, AE, Uliasz, M, Richardson, S, Miles, N, Andrews, AE, Sweeney, C, Diaz, LI, Martins, D, Shepson, PB and Davis, KJ. 2012. Constraining the CO2 budget of the corn belt: Exploring uncertainties from the assumptions in a mesoscale inverse system. Atmos Chem Phys 12(1): 337–354. DOI: 10.5194/acp-12-337-2012

33. Lin, JC and Gerbig, C. 2005. Accounting for the effect of transport errors on tracer inversions. Geophys Res Lett 32(1): 1–5. DOI: 10.1029/2004GL021127

34. Lin, JC, Gerbig, C, Wofsy, SC, Andrews, AE, Daube, BC, Davis, KJ and Grainger, CA. 2003. A near-field tool for simulating the upstream influence of atmospheric observations: The Stochastic Time-Inverted Lagrangian Transport (STILT) model. J Geophys Res 108(D16): 4493. DOI: 10.1029/2002JD003161

35. Lin, JC, Mitchell, L, Crosman, E, Mendoza, DL, Buchert, M, Bares, R, Fasoli, B, Bowling, DR, Pataki, D, Catharine, D, Strong, C, Gurney, KR, Patarasuk, R, Baasandorj, M, Jacques, A, Hoch, S, Horel, J and Ehleringer, J. 2018. CO2 and carbon emissions from cities: Linkages to air quality, socioeconomic activity, and stakeholders in the Salt Lake City urban area. BAMS 99(11): 2325–2339. DOI: 10.1175/BAMS-D-17-0037.1

36. Lombardi, M, Rana, R, Pazienza, P and Tricase, K. 2014. The European Policy for the Sustainability of Urban Areas and the “Covenant of Mayors” Initiative: A Case Study. In: Salomone, R and Saija, G (eds.). Pathways to Environmental Sustainability . Cham: Springer. DOI: 10.1007/978-3-319-03826-1_18

37. Mallia, DV, Kochanski, A, Urbanski, S and Lin, JC. 2018. Optimizing Smoke and Plume Rise Modeling Approaches at Local Scales. Atmosphere 9(5): 166. DOI: 10.3390/atmos9050166

38. Mallia, DV, Kochanski, A and Wu, D. 2017. Wind-Blown Dust Modeling Using a Backward-Lagrangian Particle Dispersion Model. J Appl Meteorol Climatol 56: 2845–2867. DOI: 10.1175/JAMC-D-16-0351.1

39. Mallia, DV, Lin, JC, Urbanski, S, Ehleringer, J and Nehrkorn, T. 2015. Impacts of upwind wildfire emissions on CO, CO2, and PM2.5 concentrations in Salt Lake City, Utah. J Geophys Res Atmos 120(1): 147–166. DOI: 10.1002/2014JD022472

40. McKain, K, Down, A, Raciti, SM, Budney, J, Hutyra, LR, Floerchinger, C, Herndon, SC, Nehrkorn, T, Zahniser, MS, Jackson, RB, Phillips, N and Wofsy, SC. 2015. Methane emissions from natural gas infrastructure and use in the urban region of Boston, Massachusetts. P Natl Acad Sci 112(7): 1941–1946. DOI: 10.1073/pnas.1416261112

41. McKain, K, Wofsy, SC, Nehrkorn, T, Eluszkiewicz, J, Ehleringer, JR and Stephens, BB. 2012. Assessment of ground-based atmospheric observations for verification of greenhouse gas emissions from an urban region. P Natl Acad Sci 109(22): 8423–8428. DOI: 10.1073/pnas.1116645109

42. Michalak, AM, Hirsch, A, Bruhwiler, L, Gurney, KR, Peters, W and Tans, PP. 2005. Maximum likelihood estimation of covariance parameters for Bayesian atmospheric trace gas surface flux inversions. J Geophys Res 110(D24107). DOI: 10.1029/2005JD005970

43. Michalak, AM, Randazzo, NA and Chevallier, F. 2017. Diagnostic methods for atmospheric inversions of long-lived greenhouse gases. Atmos Chem Phys 17(12): 7405–7421. DOI: 10.5194/acp-17-7405-2017

44. Miles, NL, Richardson, SJ, Lauvaux, T, Davis, KJ, Balashov, NV, Deng, A, Turnbull, JC, Sweeney, C, Gurney, KR, Patarasuk, R, Razlivanov, I, Cambaliza, MOL and Shepson, PB. 2017. Quantification of urban atmospheric boundary layer greenhouse gas dry mole fraction enhancements in the dormant season: Results from the Indianapolis Flux Experiment (INFLUX). Elem Sci Anth 5: 27. DOI: 10.1525/elementa.127

45. Mitchell, LE, Crosman, ET, Jacques, AA, Fasoli, B, Leclair-Marzolf, L, Horel, J, Bowling, DR, Ehleringer, JR and Lin, JC. 2018b. Monitoring of greenhouse gases and pollutants across an urban area using a light-rail public transit platform. Atmos Environ 187: 9–23. DOI: 10.1016/j.atmosenv.2018.05.044

46. Mitchell, LE, Lin, JC, Bowling, DR, Pataki, DE, Strong, C, Schauer, AJ, Bares, R, Bush, SE, Stephens, BB, Mendoza, D, Mallia, D, Holland, L, Gurney, KR and Ehleringer, JR. 2018a. Long-term urban carbon dioxide observations reveal spatial and temporal dynamics related to urban characteristics and growth. P Natl Acad Sci 115(12): 2912–2917. DOI: 10.1073/pnas.1702393115

47. Nassar, R, Napier-Linton, L, Gurney, KR, Andres, RJ, Oda, T, Vogel, FR and Deng, F. 2013. Improving the temporal and spatial distribution of CO2 emissions from global fossil fuel emission data sets. J Geophys Res 118(2): 917–933. DOI: 10.1029/2012JD018196

48. National Research Council (NRC) Committee on Methods for Estimating Greenhouse Gas Emissions. 2010. Verifying greenhouse gas emissions: Methods to support international climate agreements . Washington, DC: National Academies Press.

49. Nehrkorn, T, Eluszkiewicz, J, Wofsy, SC, Lin, JC, Gerbig, C, Longo, M and Freitas, S. 2010. Coupled Weather Research and Forecasting – Stochastic Time-Inverted Lagrangian Transport (WRF – STILT) model. Meteor Atmos Phys 107(1–2): 51–64. DOI: 10.1007/s00703-010-0068-x

50. Newman, S, Xu, X, Gurney, KR, Hsu, YK, Li, KF, Jiang, X, Keeling, R, Feng, S, O’Keefe, D, Patarasuk, R, Wong, KW, Rao, P, Fischer, ML and Yung, YL. 2016. Toward consistency between trends in bottom-up CO2 emissions and top-down atmospheric measurements in the Los Angeles megacity. Atmos Chem Phys 16(6): 3843–3863. DOI: 10.5194/acp-16-3843-2016

51. Nickless, A, Rayner, PJ, Engelbrecht, F, Brunke, E-G, Erni, B and Scholes, RJ. 2018. Estimates of CO2 fluxes over the City of Cape Town, South Africa, through Bayesian inverse modelling. Atmos Chem Phys 18(7): 4765–4801. DOI: 10.5194/acp-18-4765-2018

52. Oda, T, Lauvaux, T, Lu, D, Rao, P, Miles, NL, Richardson, SJ and Gurney, KR. 2017. On the impact of granularity of space-based urban CO2 emissions in urban atmospheric inversions: A case study for Indianapolis, IN. Elem Sci Anth 5: 28. DOI: 10.1525/elementa.146

53. Oda, T and Maksyutov, S. 2011. A very high-resolution (1 km ×1 km) global fossil fuel CO2 emission inventory derived using a point source database and satellite observations of nighttime lights. Atmos Chem Phys , 11(2): 543–556. DOI: 10.5194/acp-11-543-2011

54. Oda, T and Maksyutov, S. 2015. ODIAC Fossil Fuel CO2 Emissions Dataset (Version name: ODIAC2016), Center for Global Environmental Research, National Institute for Environmental Studies, Japan. DOI: 10.17595/20170411.001

55. Oda, T, Maksyutov, S and Andres, RJ. 2018. The Open-source Data Inventory for Anthropogenic CO2, version 2016 (ODIAC2016): a global monthly fossil fuel CO2 gridded emissions data product for tracer transport simulations and surface flux inversions. Earth Syst Sci Data 10: 87–107. DOI: 10.5194/essd-10-87-2018

56. Patarasuk, R, Gurney, KR, O’Keeffe, D, Song, Y, Huang, J, Rao, P, Buchert, M, Lin, JC, Mendoza, D and Ehleringer, JR. 2016. Urban high-resolution fossil fuel CO2 emissions quantification and exploration of emission drivers for potential policy applications. Urban Ecosyst 19(3): 1013–1039. DOI: 10.1007/s11252-016-0553-1

57. Peylin, P, Houweling, S, Krol, MC, Karstens, U, Rodenbeck, C, Geels, C, Vermeulen, A, Badawy, B, Aulagnier, C, Pregger, T, Delage, F, Pieterse, G, Ciais, P and Heimann, M. 2011. Importance of fossil fuel emission uncertainties over Europe for CO2 modeling: Model intercomparison. Atmos Chem Phys 11(13): 6607–6622. DOI: 10.5194/acp-11-6607-2011

58. Plumer, B. 2018. Climate Negotiators Reach an Overtime Deal to Keep Paris Pact Alive. New York Times [online] Day: 15 Month: December 2018. Available at: https://www.nytimes.com/2018/12/15/climate/cop24-katowice-climate-summit.html. [Accessed 20 December 2018].

59. Salon, D, Sperling, D, Meier, A, Murphy, S, Gorham, R and Barrett, J. 2010. City carbon budgets: A proposal to align incentives for climate-friendly communities. Energy Policy 38(4): 2032–2041. DOI: 10.1016/j.enpol.2009.12.005

60. Sargent, M, Barrera, Y, Nehrkorn, T, Hutyra, LR, Gately, CK, Jones, T, McKain, K, Sweeney, C, Hegarty, J, Hardiman, B, Wang, JA and Wofsy, SC. 2018. Anthropogenic and biogenic CO2 fluxes in the Boston urban region. P Natl Acad Sci 115(29): 7491–7496. DOI: 10.1073/pnas.1803715115

61. Tabuchi, H and Fountain, H. 2017. Bucking Trump, these cities, states and companies commit to Paris accord. New York Times . [online] Day: 1 Month: June 2017. Available at: https://www.nytimes.com/2017/06/01/climate/american-cities-climate-standards.html [Accessed 21 August 2018].

62. Tarantola, A. 1987. Inverse Problem Theory Methods for Data Fitting and Model Parameter Estimation . Amsterdam, NL: Elsevier Science B.V.

63. Turnbull, JC, Karion, A, Davis, KJ, Lauvaux, T, Miles, NL, Richardson, SJ, Sweeney, C, McKain, K, Lehman, SJ, Gurney, KR, Patarasuk, R, Liang, J, Shepson, PB, Heimburger, A, Harvey, R and Whetstone, J. 2019. Synthesis of Urban CO2 Emission Estimates from Multiple Methods from the Indianapolis Flux Project (INFLUX). Environ Sci Tech 53(1): 287–295. DOI: 10.1021/acs.est.8b05552

64. Turner, AJ, Shusterman, AA, McDonald, BC, Teige, V, Harley, RA and Cohen, RC. 2016. Network design for quantifying urban CO2 emissions: Assessing trade-offs between precision and network density. Atmos Chem Phys 16: 13465–13475. DOI: 10.5194/acp-16-13465-2016

65. United Nations – Department of Economic and Social Affairs (UN-DESA) Population Division. 2014. World Urbanization Prospects: The 2014 Revision, Highlights . New York: United Nations.

66. US Conference of Mayors. 2017. Mayors Undeterred by Paris Climate Accord Withdrawal . [press release] Day: 2 Month: June 2017. Available at https://www.usmayors.org/2017/06/02/mayors-undeterred-by-paris-climate-accord-withdrawal/ [Accessed 30 August 2018].

67. Wang, R. 2012. Leaders, Followers, and Laggards: Adoption of the US Conference of Mayors Climate Protection Agreement in California. Environ Plann C 30(6): 1116–1128. DOI: 10.1068/c1122

68. Wu, K, Lauvaux, T, Davis, KJ, Deng, A, Lopez Coto, I, Gurney, KR and Patarasuk, R. 2018. Joint inverse estimation of fossil fuel and biogenic CO2 fluxes in an urban environment: An observing system simulation experiment to assess the impact of multiple uncertainties. Elem Sci Anth 6(1): 17. DOI: 10.1525/elementa.138

69. Wu, L, Bocquet, M, Lauvaux, T, Chevallier, F, Rayner, P and Davis, K. 2011. Optimal representation of source-sink fluxes for mesoscale carbon dioxide inversion with synthetic data. J Geophys Res 116(D21): 304. DOI: 10.1029/2011JD016198

70. Yadav, V and Michalak, AM. 2013. Improving computational efficiency in large linear inverse problems: an example from carbon dioxide flux estimation. Geosci Model Dev 6: 583–590. DOI: 10.5194/gmd-6-583-2013