Top-down, data-driven models possess ample power to improve the accuracy of bottom-up carbon dioxide (CO_{2}) 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 CO_{2} fluxes. Here we present a Bayesian inverse modeling framework over Salt Lake City, Utah, which utilizes available CO_{2} 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 CO_{2} 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.
Salt Lake CityUrban CO_{2} emissionsBayesian inverse modelingOSSEError covariance parametersSynthetic data1. Introduction
Anthropogenic greenhouse gases (GHGs), particularly carbon dioxide (CO_{2}), 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 CO_{2} 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 CO_{2} 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 CO_{2} 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 CO_{2} 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 CO_{2} 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 CO_{2} enhancements with observed signals at 5 sites around the valley. However, no follow-up studies have taken advantage of the SLV CO_{2} 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 CO_{2} 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 CO_{2} 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 CO_{2} 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.
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 CO_{2} 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 CO_{2} 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 CO_{2} (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 CO_{2} 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 CO_{2} 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.
Prior (ODIAC), “True” (Hestia), Posterior, and residual emission maps. CO_{2} 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 CO_{2} 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 CO_{2} 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 3^{rd} and 30^{th} for a total of 4032 simulations. STILT has been widely used to interpret CO_{2} 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.
WRF-STILT footprint (influence) averages. (a) WRF-STILT footprints are shown on a log_{10} 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):
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), s_{p} 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 s_{p}. 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_{τ} * m_{s} where m_{τ} = number of timesteps (114) and m_{s} = number of grid cells (2386).
The posterior best estimate of emissions, ŝ, is given as the solution to the above cost-minimization function:
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. ŝ, s_{p}, s_{truth}) 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 s_{truth} 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:
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/):
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:
These matrices are computed using separation distances (X_{s}) and lag-times (X_{τ}) between cells and timesteps respectively, divided by their corresponding correlation range parameters l_{s} 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 l_{s}, 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:
Individual components of the model-data errors are summarized in Table 1, and are described as follows:
Description of components within the model-data mismatch matrix (R). All components are first given as hourly standard error (2^{nd} column) and, if applicable, correlation is applied between same-day observations using a decay time-scale parameter (3^{rd} column). Component errors are then aggregated to express observational error in enhancements that represent an average for a given afternoon (4^{th} column). Standard error values (in ppm) are squared to express variance (in ppm^{2}) 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
R_{part}
0.1 ppm (fixed value)
uncorrelated
0.041 ppm
R_{aggr}
16.7% and 25% of mean enhancement (spatial and temporal aggregation)
correlated with no decay in same afternoon
1.176 ppm
R_{transPBL}
7% of mean enhancement
correlated with no decay in same afternoon
0.197 ppm
R_{transWIND}
35% of mean enhancement
correlated with decay time-scale = 2.8 hours
0.881 ppm
R_{part} 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. R_{aggr} 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 (R_{transWIND}), vertical mixing-layer height error (R_{transPBL}), and unresolved eddy turbulence error (R_{eddy}). R_{transWIND} 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 CO_{2} with and without wind errors are then used to estimate the impacts of transport errors.
R_{transPBL} 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, R_{eddy} is accounted for within the above transport and aggregation errors and is neglected here.
Remaining errors in the equation (R_{instr}, R_{bkgd}, and R_{bio}) 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 R_{aggr}, R_{transPBL}, and R_{transWIND} for hourly observations within each afternoon (but not across days or sites). R_{aggr} and R_{transPBL} are considered to be fully correlated within each afternoon (as these error components are not expected to be random within a given afternoon), while R_{transWIND} 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 (l_{s} 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.
Sensitivity tests performed for this inversion analysis. Bolded values denoted with an asterisk (*) indicate baseline value from the inversion. R_{diagonal} 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
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 (r^{2}) 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):
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.
Afternoon timeseries of prior, true, and posterior (optimized) domain-average emissions. CO_{2} 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.
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, s_{p} = 3.54, and s_{truth} = 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, s_{p} = 3.45 and s_{truth} = 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.
Modeled observations (prior, posterior, and synthetic) at two sample sites. Derived afternoon CO_{2} 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, l_{s}, as shown in Figure 7, results in additional spread in the spatial extent of uncertainty reduction. With small length-scales (e.g. l_{s} = 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.
Uncertainty reduction maps with varying spatial correlation length scale (l_{s}). Maps are shown for l_{s} = 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
Space- and time-averaged CO_{2} emissions vs. spatial and temporal correlation length scales (l_{s} and l_{τ}). Afternoon domain-averaged CO_{2} 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 (l_{s} = 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.
Space- and time-averaged CO_{2} emissions and uncertainty reduction vs. model-data mismatch error. (a) Domain-average CO_{2} 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 χ_{r}^{2} 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 χ_{r}^{2} = 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.
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
Posterior minus True CO_{2} Emissions (Average) [μmol/(m^{2} 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 CO_{2} contribution maps for these two time periods at the RPK site (Figure 10), detailing the gridded contributions to true and prior CO_{2} 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.
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 CO_{2} 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}).
Afternoon timeseries of prior and true domain-average emissions, compared to baseline and four additional model scenarios. CO_{2} 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 CO_{2} 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 CO_{2} 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:
Text S1. Overview of computational methods for efficiency. DOI: https://doi.org/10.1525/elementa.375.s1
Text S2. Description of Monte-Carlo methods to determine grid-averaged emissions results. DOI: https://doi.org/10.1525/elementa.375.s1
Text S3. Additional details regarding error covariance and emissions calculations. DOI: https://doi.org/10.1525/elementa.375.s1
Text S4. Results using alternative methods for synthetic data generation. DOI: https://doi.org/10.1525/elementa.375.s1
Figure S1. Monte-Carlo simulation of domain-averaged emissions and reduced chi-square statistic. DOI: https://doi.org/10.1525/elementa.375.s1
Figure S3. Uncertainty reduction vs. spatial and temporal correlation length scales (l_{s} and l_{τ}). DOI: https://doi.org/10.1525/elementa.375.s1
Figure S4. Posterior minus true emissions averaged over September 2015. DOI: https://doi.org/10.1525/elementa.375.s1
Figure S5. Average footprint areas for September 8–16 and 19–22. DOI: https://doi.org/10.1525/elementa.375.s1
Table S1. Comparison of posterior modeled enha ncements vs. synthetic data at different sites across the Salt Lake Valley. DOI: https://doi.org/10.1525/elementa.375.s1
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.
Andres, RJ, Boden, TA and Higdon, DM. 2016. Gridded uncertainty in fossil fuel carbon dioxide emission maps, a CDIAC example. 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 CO_{2} emissions from atmospheric concentration measurements. 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, NO_{x} and CO_{2} and their impacts. 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. 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 CO_{2} emission from Houston, Texas. 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. C40 Cities. 2016. 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. Dhakal, S and Shrestha, RM. 2010. Bridging the research gaps for carbon emissions and their management in cities. Enting, IG. 2002. Enting, IG and Mansbridge, JV. 1989. Seasonal sources and sinks of atmospheric CO_{2} Direct inversion of filtered data. Enting, IG and Mansbridge, JV. 1991. Latitudinal distribution of sources and sinks of CO_{2}: results of an inversion study. 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). 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. Gately, CK and Hutyra, LR. 2017. Large uncertainties in urban-scale carbon emissions. Gerbig, C, Dolman, AJ and Heimann, M. 2009. On observational and modelling strategies targeted at regional carbon exchange over continents. Gerbig, C, Körner, S and Lin, JC. 2008. Vertical mixing in atmospheric tracer transport models: error characterization and propagation. 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?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. 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 CO_{2} Emissions at the Urban Scale in Four US Urban Areas. 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 FFCO_{2} emissions estimate in a large US urban area. Gurney, KR, Razlivanov, I, Song, Y, Zhou, Y, Benes, B and Abdul-Massih, M. 2012. Quantification of fossil fuel CO_{2} emissions on the building/street scale for a large U.S. city. 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 CO_{2}, CH_{4}, and CO flux estimates using ground- and space-based remote sensing and a Lagrangian model. Hogue, S, Marland, E, Andres, RJ, Marland, G and Woodard, D. 2016. Uncertainty in gridded CO_{2} emissions estimates. 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. International Energy Agency (IEA). 2008. 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.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.Kort, EA, Angevine, WM, Duren, R and Miller, CE. 2013. Surface observations for monitoring urban fossil fuel CO_{2} emissions: Minimum site location requirements for the Los Angeles megacity. 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 CO_{2} emissions during the dormant season of the Indianapolis flux experiment (INFLUX). 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 CO_{2} from Davos, Switzerland: The first real-time monitoring system using an atmospheric inversion technique. 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 CO_{2} budget of the corn belt: Exploring uncertainties from the assumptions in a mesoscale inverse system. Lin, JC and Gerbig, C. 2005. Accounting for the effect of transport errors on tracer inversions. 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. 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. CO_{2} and carbon emissions from cities: Linkages to air quality, socioeconomic activity, and stakeholders in the Salt Lake City urban area. 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.). Mallia, DV, Kochanski, A, Urbanski, S and Lin, JC. 2018. Optimizing Smoke and Plume Rise Modeling Approaches at Local Scales. Mallia, DV, Kochanski, A and Wu, D. 2017. Wind-Blown Dust Modeling Using a Backward-Lagrangian Particle Dispersion Model. Mallia, DV, Lin, JC, Urbanski, S, Ehleringer, J and Nehrkorn, T. 2015. Impacts of upwind wildfire emissions on CO, CO_{2}, and PM2.5 concentrations in Salt Lake City, Utah. 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. 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. 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. Michalak, AM, Randazzo, NA and Chevallier, F. 2017. Diagnostic methods for atmospheric inversions of long-lived greenhouse gases. 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). 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. 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. Nassar, R, Napier-Linton, L, Gurney, KR, Andres, RJ, Oda, T, Vogel, FR and Deng, F. 2013. Improving the temporal and spatial distribution of CO_{2} emissions from global fossil fuel emission data sets. National Research Council (NRC) Committee on Methods for Estimating Greenhouse Gas Emissions. 2010. 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. 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 CO_{2} emissions and top-down atmospheric measurements in the Los Angeles megacity. Nickless, A, Rayner, PJ, Engelbrecht, F, Brunke, E-G, Erni, B and Scholes, RJ. 2018. Estimates of CO_{2} fluxes over the City of Cape Town, South Africa, through Bayesian inverse modelling. 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 CO_{2} emissions in urban atmospheric inversions: A case study for Indianapolis, IN. Oda, T and Maksyutov, S. 2011. A very high-resolution (1 km ×1 km) global fossil fuel CO_{2} emission inventory derived using a point source database and satellite observations of nighttime lights. Oda, T and Maksyutov, S. 2015. ODIAC Fossil Fuel CO_{2} Emissions Dataset (Version name: ODIAC2016), Center for Global Environmental Research, National Institute for Environmental Studies, Japan. DOI: 10.17595/20170411.001Oda, T, Maksyutov, S and Andres, RJ. 2018. The Open-source Data Inventory for Anthropogenic CO_{2}, version 2016 (ODIAC2016): a global monthly fossil fuel CO_{2} gridded emissions data product for tracer transport simulations and surface flux inversions. 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 CO_{2} emissions quantification and exploration of emission drivers for potential policy applications. 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 CO_{2} modeling: Model intercomparison. Plumer, B. 2018. Climate Negotiators Reach an Overtime Deal to Keep Paris Pact Alive. 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. 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 CO_{2} fluxes in the Boston urban region. Tabuchi, H and Fountain, H. 2017. Bucking Trump, these cities, states and companies commit to Paris accord. Tarantola, A. 1987. 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 CO_{2} Emission Estimates from Multiple Methods from the Indianapolis Flux Project (INFLUX). Turner, AJ, Shusterman, AA, McDonald, BC, Teige, V, Harley, RA and Cohen, RC. 2016. Network design for quantifying urban CO_{2} emissions: Assessing trade-offs between precision and network density. United Nations – Department of Economic and Social Affairs (UN-DESA) Population Division. 2014. US Conference of Mayors. 2017. Wang, R. 2012. Leaders, Followers, and Laggards: Adoption of the US Conference of Mayors Climate Protection Agreement in California. 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 CO_{2} fluxes in an urban environment: An observing system simulation experiment to assess the impact of multiple uncertainties. 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. Yadav, V and Michalak, AM. 2013. Improving computational efficiency in large linear inverse problems: an example from carbon dioxide flux estimation.