Evaluating consistency between total column CO₂ retrievals from OCO-2 and the in situ network over North America: implications for carbon flux estimation

Feedbacks between the climate system and the carbon cycle represent a key source of uncertainty in model projections of Earth's climate, in part due to our inability to directly measure large-scale biosphere–atmosphere carbon fluxes. In situ measurements of the CO₂ mole fraction from surface flasks, towers, and aircraft are used in inverse models to infer fluxes, but measurement networks remain sparse, with limited or no coverage over large parts of the planet. Satellite retrievals of total column CO₂ ($X_{{\mathrm{CO}}_{\mathrm{2}}}$), such as those from NASA's Orbiting Carbon Observatory-2 (OCO-2), can potentially provide unprecedented global information about CO₂ spatiotemporal variability. However, for use in inverse modeling, data need to be extremely stable, highly precise, and unbiased to distinguish abundance changes emanating from surface fluxes from those associated with variability in weather. Systematic errors in $X_{{\mathrm{CO}}_{\mathrm{2}}}$ have been identified and, while bias correction algorithms are applied globally, inconsistencies persist at regional and smaller scales that may complicate or confound flux estimation. To evaluate $X_{{\mathrm{CO}}_{\mathrm{2}}}$ retrievals and assess potential biases, we compare OCO-2 v10 retrievals with in situ data-constrained $X_{{\mathrm{CO}}_{\mathrm{2}}}$ simulations over North America estimated using surface fluxes and boundary conditions optimized with observations that are rigorously calibrated relative to the World Meteorological Organization X2007 CO₂ scale. Systematic errors in simulated atmospheric transport are independently evaluated using unassimilated aircraft and AirCore profiles. We find that the global OCO-2 v10 bias correction shifts the distribution of retrievals closer to the simulated $X_{{\mathrm{CO}}_{\mathrm{2}}}$, as intended. Comparisons between bias-corrected and simulated $X_{{\mathrm{CO}}_{\mathrm{2}}}$ reveal differences that vary seasonally. Importantly, the difference between simulations and retrievals is of the same magnitude as the imprint of recent surface flux in the total column. This work demonstrates that systematic errors in OCO-2 v10 retrievals of $X_{{\mathrm{CO}}_{\mathrm{2}}}$ over land can be large enough to confound reliable surface flux estimation and that further improvements in retrieval and bias correction techniques are essential. Finally, we show that independent observations, especially vertical profile data, such as those from the National Oceanic and Atmospheric Administration aircraft and AirCore programs are critical for evaluating errors in both satellite retrievals and carbon cycle models.

1 Introduction

Interannual variability in the growth rate of atmospheric CO₂ is largely driven by variability in uptake and release by terrestrial ecosystems (Heimann and Reichstein, 2008; Piao et al., 2020). Oceanic fluxes also respond to variability in climate (e.g., DeVries et al., 2019; Riebesell et al., 2009), but the amplitude of oceanic flux variability is thought to be considerably less than for terrestrial fluxes. While individual component fluxes (e.g., photosynthesis) are currently not directly measurable at scales larger than leaf or soil chambers, well-calibrated and precise measurements of CO₂ have allowed us to track the accumulation of this greenhouse gas in the atmosphere and associated radiative feedbacks on the global climate, as well as its spatiotemporal variability (e.g., Tans et al., 1989). These measurements continue to provide valuable insights into surface flux processes and feedbacks (e.g., Tans et al., 1990; Ballantyne et al., 2012; Keeling et al., 2017; Arora et al., 2020). Observed spatial and temporal gradients in the CO₂ mole fraction (relative to dry air) can be combined with a numerical model of atmospheric transport to infer surface fluxes (i.e., exchange of CO₂ between the atmosphere and the underlying ocean or land surface), in a top-down or inverse modeling framework (e.g., Peters et al., 2007; Gurney et al., 2002). An ever-increasing global greenhouse gas measurement network and progress in modeling techniques have tremendously improved our understanding of surface processes. However, measurement networks remain sparse and continue to undersample large parts of the world, including large parts of North America, which can limit our understanding of surface flux processes in those regions. Furthermore, incompatibility across data sets that arises from inconsistent calibrations or systematic errors can significantly corrupt surface flux estimates, leading to inaccurate models of carbon–climate interactions and subsequent errors in climate forecasts.

Satellite retrievals of total column CO₂ mole fraction ($X_{{\mathrm{CO}}_{\mathrm{2}}}$), such as those from NASA's Orbiting Carbon Observatory-2 (OCO-2), have the potential to provide unprecedented information about spatiotemporal patterns and variability in the Earth's atmosphere. However, space-based observations of $X_{{\mathrm{CO}}_{\mathrm{2}}}$ must be extremely stable, highly precise, and free from bias to detect and quantify abundance changes caused by a change in surface fluxes (Rayner and O'Brien, 2001; Olsen, 2004; Miller et al., 2007; Houweling et al., 2004). Regional flux of terrestrial net non-fire ecosystem exchange (NEE) of CO₂ can be small, as it is composed of two opposing fluxes (photosynthesis and respiration). Furthermore, NEE is ubiquitous on the terrestrial surface (unlike, for example, spatially discrete point sources from industrial emissions). Lastly, CO₂ has a long lifetime in the atmosphere (and is therefore well mixed). Together, these imply that the CO₂ mole fraction changes due to NEE over large regions (e.g., temperate North America) can be hard to distinguish from variability in the CO₂ mole fraction resulting from flux processes and transport upwind. CO₂ mole fraction changes in $X_{{\mathrm{CO}}_{\mathrm{2}}}$ from NEE at the surface are diluted over the path length of the atmosphere and are largely obscured by meteorological variability (Basu et al., 2018; Feng et al., 2019).

In situ measurements that comprise global networks, such as the National Oceanic and Atmospheric Administration (NOAA)'s Global Greenhouse Reference Network (https://www.esrl.noaa.gov/gmd/ccgg/ggrn.php, last access: 20 September 2021), are rigorously evaluated and carefully calibrated relative to the World Meteorological Organization (WMO) calibration scale (data used here are reported on the X2007 scale), thus ensuring the fidelity of these measurements over timescales of seasons to decades (Andrews et al., 2014; Hall et al., 2021). The open-path nature of space-based $X_{{\mathrm{CO}}_{\mathrm{2}}}$ measurements does not allow for direct calibration. Satellite retrievals require a forward model of radiative transfer that is run through an inversion system along with satellite-obtained absorption spectra of atmospheric O₂ and CO₂ to infer $X_{{\mathrm{CO}}_{\mathrm{2}}}$. While a great amount of progress has been made to identify and eliminate sources of uncertainty emanating from this chain of processes, e.g., in the molecular absorption model and spectroscopy (Thompson et al., 2012; Payne et al., 2020; Hobbs et al., 2020), considerable sources of uncertainty remain. These are attributed to the presence of aerosols in the column (Connor et al., 2016), clouds and cloud shadows (Massie et al., 2021), interference of jointly retrieved parameters (Kulawik et al., 2019), surface properties, and details of the instrumentation. Connor et al. (2016) estimate that aerosol-dependant biases for retrievals over land may be as large as ∼ 2 ppm (parts per million dry air mole fraction). Moreover, sensors typically degrade over time, and limited information is available to characterize resulting time-dependent systematic errors. Post-launch data corrections are routinely performed and have generally reduced $X_{{\mathrm{CO}}_{\mathrm{2}}}$ bias. For example, mean bias in land nadir $X_{{\mathrm{CO}}_{\mathrm{2}}}$ relative to the Total Carbon Column Observing Network (TCCON) in the v8 product was reduced from 0.72±1.22 to 0.30±1.04 ppm (O'Dell et al., 2018), and a correction in a geolocation error resulted in a decrease in across-scene standard deviation from 1.35 ppm in v8 to 0.74 ppm in the v9 data product (Kiel et al., 2019).

Currently, satellite-derived $X_{{\mathrm{CO}}_{\mathrm{2}}}$ retrievals are linked to the WMO scale most directly through a set of in situ profiles obtained over a network consisting of 26 ground-based Fourier transform infrared spectrometers that comprise the TCCON (Wunch et al., 2017). However, the TCCON itself provides remotely sensed information about $X_{{\mathrm{CO}}_{\mathrm{2}}}$, and comparison with aircraft profiles have revealed errors ∼ 1 ppm (Wunch et al., 2011). Seasonal and site-dependent biases associated with the validation of OCO-2 retrievals via TCCON have also been reported (Wunch et al., 2017). OCO-2 retrievals are additionally corrected for bias by comparing with 4D CO₂ mole fraction fields from global inverse models and a small area approximation, but both methods are prone to smoothing across fine-scale variability in $X_{{\mathrm{CO}}_{\mathrm{2}}}$ (O'Dell et al., 2018; Corbin et al., 2008). While bias correction generally reduces inferred surface flux uncertainty when retrievals are assimilated in atmospheric inversions (Basu et al., 2013), even small retrieval errors can lead to large errors in inferred flux (Takagi et al., 2014; Chevallier et al., 2014; Villalobos et al., 2020).

Thus, a method to routinely evaluate satellite retrievals is necessary. In this study, we propose such a method that takes advantage of the relatively dense in situ network of surface, tower-based, and vertical profile CO₂ mole fraction measurements over North America and leverage an ensemble of optimized flux estimates derived using the high-resolution CarbonTracker–Lagrange CO₂ inverse modeling framework (Hu et al., 2019). We demonstrate an approach for constructing $X_{{\mathrm{CO}}_{\mathrm{2}}}$ that offers optimal consistency with in situ measurements of CO₂ dry air mole fraction calibrated relative to the WMO X2007 scale (Hall et al., 2021). We compare our simulated $X_{{\mathrm{CO}}_{\mathrm{2}}}$ retrievals with the v10 OCO-2 product, hereafter $X_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{ret}}$, evaluate the extent to which observed differences are consistent with rigorous uncertainties on the simulated fields, and potentially correct biases in satellite retrievals. We essentially use the CarbonTracker–Lagrange modeling framework to interpolate the existing in situ measurements to the time and location of $X_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{ret}}$. The in situ measurement uncertainty is generally ∼ 0.15 ppm (Andrews et al., 2014), and we have carefully accounted for uncertainty in regional boundary conditions and uncertainties in the optimized fluxes. While a single realization of the simulated atmospheric transport is used here (a limitation of the approach that we aim to address in future work by using transport ensembles), we evaluate transport uncertainty using independent vertical in situ profiles of CO₂. CarbonTracker–Lagrange simulations are driven by meteorological simulations from the Weather Research Forecast model system optimized for particle dispersion modeling (Nehrkorn et al., 2010), with a resolution of 10 km over continental USA and 30 km over the rest of North America, considerably higher than global in situ informed simulations that have so far been used for OCO-2 evaluation (Kiel et al., 2019; Crowell et al., 2019; Miller and Michalak, 2020). Simulated $X_{{\mathrm{CO}}_{\mathrm{2}}}$ are compared with OCO-2 v10 retrievals, both before and after global bias correction, thus providing an independent evaluation of the global bias correction over North America. In principle, differences can result from errors in the retrievals or from errors in the CarbonTracker–Lagrange modeling framework, but we show that for certain seasons these differences are unlikely to result from the latter. North America is a useful test bed for evaluating consistencies and for developing improved model simulations and retrieval bias correction strategies, given the relatively dense sampling network (compared to other regions) and that the best surface flux estimates are likely to come from approaches that combine in situ measurements and satellite retrievals (Basu et al., 2013; Fischer et al., 2017; Byrne et al., 2020). The interagency North American Carbon Program (Wofsy et al., 2002) has supported intensely focused research for nearly 2 decades and has resulted in a wealth of data sets and model-to-data fusion activities that have informed model development.

2 Methods

Simulated $X_{{\mathrm{CO}}_{\mathrm{2}}}$, hereafter $X_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{sim}}$, is constructed by estimating impact of different surface fluxes (${\mathrm{\Delta }}_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{flux}}$) on the total column. This involves imposing a time, latitude, longitude, and altitude-dependent lateral boundary condition or background which accounts for changes in $X_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{sim}}$ originating outside our model domain. The chain of events that link the in situ data with the simulations is shown in Fig. 1. We use multiple ensembles of surface flux and boundary conditions to assess uncertainty in each. Comparisons with independent unassimilated aircraft and AirCore profiles are used to assess combined random and systematic errors in surface flux, background, and atmospheric transport.

Figure 1 Flowchart linking in situ measurements of CO₂ to simulated columns.

2.1 Convolution method

We follow the recommended protocol for comparing satellite retrievals with modeled CO₂ columns from the Atmospheric CO₂ Observations from Space (ACOS) retrieval algorithm v10 (O'Dell et al., 2012).

$$\begin{array}{}\text{(1)}& X_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{sim}}=\sum _{i=\mathrm{1}}^N{w}_i\left[a_i\cdot {\mathit{\chi }}_{{\mathrm{CO}}_{\mathrm{2}},i}^{\mathrm{model}}+(\mathrm{1}-a_i)\cdot {\mathit{\chi }}_{{\mathrm{CO}}_{\mathrm{2}},i}^{\mathrm{pri}}\right].\end{array}$$

Here, $X_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{sim}}$ (parts per million; hereafter ppm), the total column CO₂, is computed as a pressure-weighted sum of the modeled column ($X_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{model}}={\sum }_{i=\mathrm{1}}^N{w}_i$. ${\mathit{\chi }}_{{\mathrm{CO}}_{\mathrm{2}},i}^{\mathrm{model}}$), comprising N model (i.e., not OCO-2) levels from the surface to the top of the atmosphere (0.01 hPa). ${\mathit{\chi }}_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{model}}$ is convolved with the OCO-2 averaging kernel profile (a_i) and the OCO-2 prior profile (${\mathit{\chi }}_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{pri}}$) and summed according to a pressure weighting function (w; identical to h in Connor et al., 2008). w is calculated as follows:

$$\begin{array}{}\text{(2)}& \begin{array}{rl}w& =\sum _{i=\mathrm{1}}^N{w}_i=\sum _{i=\mathrm{1}}^N\left|\left(-p_i+{\displaystyle \frac{p_{i+\mathrm{1}}-p_i}{\ln \left(\frac{p_{i+\mathrm{1}}}{p_i}\right)}}\right)\right.\\ & \left.+\left(-p_i+{\displaystyle \frac{p_i-p_{i-\mathrm{1}}}{\ln \left(\frac{p_i}{p_{i-\mathrm{1}}}\right)}}\right)\right|{\displaystyle \frac{\mathrm{1}}{p_{\mathrm{surf}}}},\end{array}\end{array}$$

where p_i and p_surf are WRF-modeled pressure at level i and at the surface, respectively. The profile sum of w is always unity. Additionally, since $X_{{\mathrm{CO}}_{\mathrm{2}}}$ is reported as dry air mole fraction, WRF total pressure is converted to dry air pressure at all receptor levels.

${\mathit{\chi }}_{{\mathrm{CO}}_{\mathrm{2}},i}^{\mathrm{model}}$ (ppm) is constructed as follows:

$$\begin{array}{}\text{(3)}& X_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{model}}=\sum _{i=\mathrm{1}}^N{\mathit{\chi }}_{{\mathrm{CO}}_{\mathrm{2}},i}^{\mathrm{bkg}}+\sum _{i=\mathrm{1}}^{N-\mathrm{3}}{\mathrm{\Delta }}_{{\mathrm{CO}}_{\mathrm{2}},i}^{\mathrm{flux}}.\end{array}$$

Here, ${\mathit{\chi }}_{{\mathrm{CO}}_{\mathrm{2}},i}^{\mathrm{bkg}}$ (ppm) represents background (i.e., lateral boundary condition; described in Sect. 2.4) and ${\mathrm{\Delta }}_{{\mathrm{CO}}_{\mathrm{2}},\mathrm{i}}^{\mathrm{flux}}$ (ppm) denotes the impact of surface flux at level i of the column. ${\mathrm{\Delta }}_{{\mathrm{CO}}_{\mathrm{2}},\mathrm{i}}^{\mathrm{flux}}$ is computed at discrete levels from the surface to 14 km, whereas ${\mathit{\chi }}_{{\mathrm{CO}}_{\mathrm{2}},i}^{\mathrm{bkg}}$ is computed at three additional levels. These additional levels represent the upper troposphere and the stratosphere, where the influence of recent surface flux is assumed to be zero. If there are cases where recent surface fluxes influence the upper tropospheric and stratospheric air, those are accounted for as part of the background estimation. This is because models used to estimate background are also constrained by in situ measurements. Note that there may still be rare cases, e.g., in the proximity of large fires (Hooghiem et al., 2020), where surface flux influence in the upper troposphere and lower stratosphere may not be captured by this approach. ${\mathrm{\Delta }}_{{\mathrm{CO}}_{\mathrm{2}},\mathrm{i}}^{\mathrm{flux}}$ is estimated as follows:

$$\begin{array}{}\text{(4)}& {\mathrm{\Delta }}_{{\mathrm{CO}}_{\mathrm{2}},\mathrm{i}}^{\mathrm{flux}}=H_i\left(s_{\mathrm{bio}}+s_{\mathrm{ff}}+s_{\mathrm{bmb}}+s_{\mathrm{ocn}}\right),\end{array}$$

where H_i ($\frac{\mathrm{ppm}{\mathrm{CO}}_{\mathrm{2}}}{\mathrm{\mu }\mathrm{mol}{\mathrm{CO}}_{\mathrm{2}}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}}$) represents the sensitivity at pressure level i of simulated $X_{{\mathrm{CO}}_{\mathrm{2}}}$ to upwind surface fluxes (detailed in Sect. 2.3), and s indicates surface flux (µmol m⁻² s⁻¹). s_bio denotes net ecosystem exchange (i.e., the sum of ecosystem photosynthesis and respiration), s_bmb denotes biomass burning, s_ff corresponds to fossil fuel emissions, and s_ocn is the net ocean–atmosphere flux. Fluxes are described in detail in Sect. 2.5.

2.2 OCO-2 retrievals and receptor selection criterion

We construct $X_{{\mathrm{CO}}_{\mathrm{2}}}$ for valid land nadir and land glint $X_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{ret}}$ from the v10 data product (Osterman et al., 2020) over North America between September 2014 and August 2015. Globally, OCO-2 retrievals are primarily obtained in three operation modes, i.e., land glint, land nadir, and ocean glint. Additionally, retrievals are obtained in target mode, for evaluation against TCCON, and a transition mode, where the sensor switches between modes. While the nadir mode is mostly used over land, over darker ocean the satellite sensor is able to receive a higher fraction of directly reflected sunlight in a separate glint mode (Eldering et al., 2017). Here, we evaluate soundings over land obtained in both the nadir and glint modes. Additionally, only retrievals that passed the aerosol and cloud-screening filers are considered (i.e., quality flagged as 0 or good).

Over North America, there can be a few thousand to tens of thousands of valid retrievals on any given day (Fig. 2a). Each retrieval covers an area of approximately 1.29 km × 2.25 km on the surface. Individual satellite retrievals are known as footprints. The satellite collects eight simultaneous footprints, and the next row of footprints are spaced 300 ms apart. For each satellite overpass (along track), we select locations every 2 s over the continental USA (i.e.,∼12 km) and 4 s (i.e.,∼24 km) over the rest of North America (Fig. 2b). These locations are usually called receptors in a Lagrangian particle dispersion model (LPDM) as they represent locations and times from which a set of particles are released and tracked back in time. In an LPDM, an ensemble of particles is released from each receptor, and the residence time of particles in the planetary boundary layer is used to calculate sensitivities describing the relationship between upwind surface fluxes and mole fraction at the receptor location. Ultimately, a library of sensitivity arrays is generated corresponding to $X_{{\mathrm{CO}}_{\mathrm{2}}}$ retrieval locations. Note that these sensitivity arrays are sometimes called footprints or influence functions. Here we avoid this use of footprints so as not to create confusion with OCO-2 footprints (i.e., scenes). This method enables improved simulation of near-field transport compared to Eulerian gridded models, as particle locations are not restricted to grid boxes and meteorological fields can be interpolated to subgrid-scale locations (Lin et al., 2003), and has been used extensively in estimating regional trace gas fluxes in inverse models using in situ measurements (e.g., Schuh et al., 2009; Gourdji et al., 2012; Lauvaux et al., 2013; Alden et al., 2016; Hu et al., 2019).

Figure 2 All valid satellite retrievals for June 2015 (left). A magnified map (right) shows all retrievals for a section of an individual sounding and locations of along-track receptors (stars). Black lines in the magnified map indicate 1^∘ latitude bands.

In this study, a vertical profile of receptors corresponding to a range of altitudes is created corresponding to the location of a single valid satellite retrieval. We preferentially select receptors that are in the middle of the OCO-2 track to provide the most spatially representative sample and minimize footprint-dependant biases (O'Dell et al., 2012). For this analysis, we assume that $X_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{sim}}$ from a given receptor location is representative of all $X_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{ret}}$ within 1 s. In total, ∼ 32 000 unique receptor profiles associated with valid retrievals are created between September 2014 and August 2015, representing ∼ 1.61 million retrievals that passed quality flags between September 2014 and August 2015. At each receptor, 24 unique $X_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{sim}}$ are created from combinations of six flux and four background ensemble members (Fig. 1). Using an ensemble of 24 flux–background combinations provides an estimate of unresolved variability in the simulations. In this analysis, we report $X_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{sim}}$ as the mean and standard deviation of these simulations. Similarly, $X_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{ret}}$ represents the mean of all OCO-2 footprints within ±1 s for a selected retrieval.

A profile of receptors for each selected satellite scene consists of discrete altitude levels approximately representing the lowermost 850 hPa of the atmosphere. Models sampled for background estimation are sampled at receptor locations to account for the rest of the column (detailed in Sect. 2.4). A similar method (the application of a regional Lagrangian model to estimate source–receptor relationships) was developed recently by Wu et al. (2018). In that study, a set of model particles distributed throughout an entire column of air (weighted to appropriately represent the retrieved total column, i.e., the column receptor) was transported and tracked backward in time, and a single surface flux sensitivity array was computed for the total column. In contrast, we here establish source–receptor relationships at discrete altitude intervals and then apply appropriate vertical weighting for the column. This method has a higher computational cost and results in larger output files but adds flexibility in the simulation; for instance, a quick recalculation can be performed if the OCO-2 averaging kernel is modified, and more importantly, our simulated profiles can be used with experimental retrieval products such as partial column retrievals (e.g., Kulawik et al., 2017).

2.3 Transport model and constructing column-weighted surface sensitivity arrays: H

We use output from the Stochastic Time-Inverted Lagrangian Transport (STILT) particle dispersion model (Lin et al., 2003) driven with high-resolution meteorological fields from a customized implementation of the Weather Research and Forecasting (WRF) model (Nehrkorn et al., 2010). We use the WRF-STILT modeling configuration developed for the CarbonTracker–Lagrange modeling framework (Hu et al., 2019).

The high spatial resolution of regional models provides an appropriate framework to investigate a relatively fine-scale structure in atmospheric transport (e.g., mesoscale) and atmospheric signatures of surface flux heterogeneity. WRF model fields are computed at 10 km spatial resolution over temperate North America and 30 km spatial resolution outside of temperate North America. STILT is run offline (i.e., driven by archived hourly WRF output), and trajectories are computed backwards in time from each receptor location. STILT surface sensitivity arrays represent simulated upwind surface flux influence for 10 d prior to each observation at $\mathrm{1}{}^{\circ }\times \mathrm{1}{}^{\circ }$ spatial and 1 h temporal resolution. A library of WRF-STILT surface sensitivity arrays is precomputed, archived, and can be efficiently convolved with independently estimated surface fluxes. This is in contrast to most modern Eulerian model CO₂ simulations, where the transport model needs to be rerun whenever a new surface flux product becomes available. STILT surface sensitivity arrays have units of ppm (µmol m⁻² s⁻¹)⁻¹.

The OCO-2 retrieval is a pressure-weighted mean of ${\mathit{\chi }}_{{\mathrm{CO}}_{\mathrm{2}}}$ obtained at 20 equally spaced pressure boundaries from the surface to the top of the atmosphere. STILT receptors, however, are specified not on the pressure grid used for OCO-2 retrievals but at fixed altitudes above ground level and with high vertical resolution where strong gradients in CO₂ are expected (e.g., near the surface or at the top of the planetary boundary layer). Additionally, $X_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{ret}}$ is computed from a combination of the signal received by the spectrometer and an a priori profile (${\mathit{\chi }}_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{pri}}$). This approach constrains the uppermost portion of the atmospheric column where the satellite sensor lacks sensitivity. The relative weights of the received signal and ${\mathit{\chi }}_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{pri}}$ are described by the column averaging kernel (a; O'Dell et al., 2012), which is computed during the retrieval and archived along with $X_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{ret}}$ and ${\mathit{\chi }}_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{pri}}$. Thus, the first step in creating column-weighted surface sensitivity arrays (H) is to interpolate ${\mathit{\chi }}_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{pri}}$ and a onto the STILT grid. Then, a pressure weighting function (Eq. 2) is applied to appropriately weight the surface sensitivity obtained from all receptors for each column. The upper 150 hPa of the atmosphere is considered as part of the lateral boundary condition (the column-weighted background; Sect. 2.4), as sensitivity to recent surface flux is assumed to be zero at these pressure levels, and the WRF-STILT framework has not been optimized for upper atmospheric simulations.

2.4 Column-weighted background: $X_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{bkg}}$

The background or lateral boundary condition is an essential component of regional models required to isolate changes in CO₂ from surface fluxes within the model domain. Boundary values need to represent synoptic variability and contributions from surface fluxes outside the model domain and may contribute significantly to uncertainty in modeled $X_{{\mathrm{CO}}_{\mathrm{2}}}$ (Feng et al., 2019). Here, we combine WRF-STILT back-trajectories with 4D global mole fraction fields from simulations that were optimized using global in situ measurements. From each receptor, 500 back-trajectories (simulating air parcels) are released and tracked backwards in time until the point at which they exit the WRF domain or for the duration of the simulation (10 d). At that coordinate (longitude, latitude, altitude, and time) a global 4D mole fraction field is sampled, and that value is considered as the background value for that particle. Background values for all 500 particles are then averaged to calculate the estimated background value for a given receptor. For background estimation, the WRF domain is subset to only include continental North America plus margins along the coast. These strategies minimize transport-related errors in the trajectories that inflate with increasing distance from the receptor and time of release. The background values for each receptor in a column are summed according to the pressure weighting function (Eq. 2).

Table 1 Global CO₂ 4D mole fraction fields used in this study. Comparisons with the Global Monitoring Laboratory (GML)'s GlobalView Plus version 6.0 data product are also presented. These comparisons are provided for all assimilable observations and all assimilable observations between 4 and 8 km. For more information on the models, see the text.

In total, there are four different global 4D mole fraction fields that are informed by in situ measurements and routinely updated and sampled as background. These are two versions of CarbonTracker (CT2016 and CT2019B; Peters et al., 2007; Jacobson et al., 2020) from NOAA's Global Monitoring Laboratory (GML), the Copernicus Atmosphere Monitoring Service reanalysis (CAMS; Chevallier et al., 2019) produced by the European Centre for Medium-Range Weather Forecasts (ECMWF), and the Jena CarboScope model from the Max Planck Institute for Biogeochemistry (Rödenbeck et al., 2020). We evaluate each model against all designated assimilable CO₂ data from NOAA GML's GlobalView Plus version 6.0 data product (Table 1; Schuldt et al., 2020). Assimilable data include assimilated and withheld data. Withheld data are qualitatively equivalent to assimilated data (i.e., they pass the same quality flags) but are excluded and used to evaluate model results in CarbonTracker (Andrew R. Jacobson​​​​​​​, personal communication). Here, over 150 000 ground, tall tower, and aircraft in situ observations spanning 2014–2015 over North America and the eastern North Pacific Ocean are used. Comparisons with these observations are provided for all assimilable observations and assimilable observations between 4 and 8 km.

2.5 Surface fluxes of CO₂: s

We sample optimized and imposed flux fields from regional and global inverse models. Net non-fire terrestrial ecosystem exchange (i.e., CO₂ fluxes from photosynthesis and respiration from autotrophic and heterotrophic sources or s_bio) are from NOAA's CarbonTracker–Lagrange (Hu et al., 2019). Briefly, CarbonTracker–Lagrange is a regional atmospheric inverse model in which biospheric fluxes for North America are optimized using surface sensitivity arrays from high-resolution WRF-STILT simulations and North American measurements of CO₂ from GlobalView Plus v2.1 (Cooperative Global Atmospheric Data Integration Project, 2016), which is composed largely of data from NOAA's Global Greenhouse Gas Reference Network and from Environment and Climate Change Canada. Observations include flask air measurements from the near-surface and aircraft and quasi-continuous in situ measurements primarily made on towers. The inversions were run with three different prior estimates of s_bio. These included two versions of the Carnegie–Ames–Stanford approach (CASA; Potter et al., 1993) biogeochemical model runs (CASA GFED-CMS and CASA GFEDv4.1) and the Combined Simple Biosphere/Carnegie–Ames–Stanford approach terrestrial carbon cycle model (SiBCASA; Schaefer et al., 2008). Prior error covariance parameters were derived from a maximum likelihood estimation (MLE) with fixed correlation scales of 1000 km and 7 d and also for optimized correlation scales, for each model run, resulting in six different posterior estimates of s_bio. Non-biospheric fluxes include imposed biomass burning (s_bmb) and fossil fuel (s_ff) fluxes and optimized ocean fluxes (s_ocn) from CT2016. We use the mean of two fossil fuel emission products (Miller and ODIAC data sets) and fire emission products (GFED4.1s and GFED-CMS) used in CT2016. All fluxes are $\mathrm{1}{}^{\circ }\times \mathrm{1}{}^{\circ }$ spatial and 3 h temporal resolution.

3 Results and discussion

3.1 Comparing simulated and satellite retrievals

For all retrievals selected over the spatiotemporal domain of this study, the impact of the OCO-2 bias correction is 2.01±0.87 ppm, and the mean difference between seasons is 0.5 ppm (1.76 in winter and spring and 2.31 ppm in summer; blue distributions in Fig. 3). Across seasons, the difference in residuals between $X_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{ret},\mathrm{bc}}$ and $X_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{sim}}$ is significantly lower than that between $X_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{sim}}$ and $X_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{ret}}$: ${\mathit{\mu }}_{\mathrm{sim}-\mathrm{ret}}=\mathrm{2.23}\pm \mathrm{1.36}$ ppm, whereas ${\mathit{\mu }}_{\mathrm{ret},\mathrm{bc}-\mathrm{sim}}=-\mathrm{0.22}\pm \mathrm{1.91}$ ppm, i.e., the OCO-2 bias correction brings the distribution of OCO-2 $X_{{\mathrm{CO}}_{\mathrm{2}}}$ substantially closer to the in situ data-constrained synthetic columns ($X_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{sim}}$), as expected. Residuals are lowest in the Northern Hemisphere summer months of June, July, and August (${\mathit{\mu }}_{\mathrm{ret},\mathrm{bc}-\mathrm{sim}}=\mathrm{0.2}\pm \mathrm{1.36}$ ppm; Fig. 3d) and highest in the winter and spring (${\mathit{\mu }}_{\mathrm{ret},\mathrm{bc}-\mathrm{sim}}=-\mathrm{0.61}$ ppm; Fig. 3a). Apart from the summer, mean $X_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{ret},\mathrm{bc}}$ over North America is lower than $X_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{sim}}$.

Figure 3 Kernel density distributions for residuals between simulated and satellite retrievals, grouped seasonally. Blue distributions show the impact of the OCO-2 bias correction on land nadir retrievals. The orange and green curves show the difference between simulated retrievals and retrievals before and after bias correction, respectively. Printed numbers report the mean and standard deviation of residuals. Uncertainties in $X_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{sim}}$ are discussed in Sect. 3.3.

3.2 Spatial patterns

To examine spatial patterns of simulated and retrieved soundings, we first sort retrievals in $\mathrm{2}{}^{\circ }\times \mathrm{2}{}^{\circ }$ bins and then average all retrievals (and simulations) in each bin for each season. The spatial extent of OCO-2 soundings varies seasonally (Fig. 4a, e, i, m). The northern extent of valid retrievals follows the solar declination, as the OCO-2 spectrometer is unable to retrieve a signal over land that is dark or blanketed by snow. Consequently, between September and March, there are few soundings north of the USA–Canada border (∼ 49^∘ N). Conversely, all of North America is observable between March and August.

Figure 4 Spatial patterns of $X_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{ret},\mathrm{bc}}$ (a, e, i, m), $X_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{sim}}$ (b, f, j, n), and $\mathrm{\Delta }{X}_{\mathrm{sim}-\mathrm{ret},\mathrm{bc}}$ (c, g, k, o) and impact of recent surface flux (${\mathrm{\Delta }}_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{flux}})$ on the $X_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{sim}}$ (d, h, l, p) for soundings grouped seasonally and plotted on a $\mathrm{2}{}^{\circ }\times \mathrm{2}{}^{\circ }$ grid. Units for all maps are in parts per million (ppm). Dotted lines in (c), (g), (k), and (o) are drawn along 40^∘ N and 100^∘ W.

Across seasons, $X_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{sim}}$ and $X_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{ret},\mathrm{bc}}$ exhibit broadly similar spatial patterns (Fig. 4a–n). However, residuals between the two reveal that $X_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{ret},\mathrm{bc}}$ is usually lower than $X_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{sim}}$, except in the summer, when a majority of retrievals in the eastern half of the continent (right of the dotted lines in Fig. 4o) are higher than the simulations. Such coherent spatial differences are not evident when examining the mean bias over the continent (Fig. 3). Importantly, the magnitude of this difference is similar to that of the recent flux signals in the total column (Fig. 4d, h, l, p). For example, the mean residuals in the northeastern quadrant during the summer (Fig. 4o) are 0.79±1.63 ppm, significantly higher than that for the entire domain (0.09±1.38 ppm) and ∼ 45 % of the impact of recent surface flux on $X_{{\mathrm{CO}}_{\mathrm{2}}}$, i.e., ${\mathrm{\Delta }}_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{flux}}$.

3.3 Examining bias and uncertainty in $X_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{sim}}$

To establish whether differences between $X_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{ret},\mathrm{bc}}$ and $X_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{sim}}$ presented above (Figs. 3, 4) are due to residual biases in the former, we characterize systematic errors and errors arising from unresolved variability in model fields used to construct $X_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{sim}}$. Cross-model standard deviations of four background models and ${\mathrm{\Delta }}_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{flux}}$ from six biospheric flux ensembles on $X_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{sim}}$ are shown in Fig. 5. Cross-model standard deviation in background fields is highest in the northwest during the summer (0.36 ppm) but usually less than 0.3 ppm. Over the entire spatiotemporal domain of the study, the standard deviation in estimates of ${\mathrm{\Delta }}_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{flux}}$ and $X_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{bkg}}$ is 0.26±0.14 and 0.28±0.15 ppm, respectively. Model spread in ${\mathrm{\Delta }}_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{flux}}$ is largest along the Pacific coast of Mexico, a region that is relatively less well constrained by the in situ network. The standard deviation in ${\mathrm{\Delta }}_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{flux}}$ is highest in the southeastern quadrant in the summer (0.32 ppm) but usually between 0.1 and 0.3 ppm. Uncertainty from model spread in flux and background is propagated in comparisons of $X_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{sim}}$ and $X_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{ret},\mathrm{bc}}$ presented earlier (Fig. 3).

Figure 5 Cross-model standard deviations flux impact (a, c, e, g) and background models (b, d, f, h). A total of six biospheric flux ensembles and four model fields for background are used in this study.

Figure 6 Simulated and observed vertical profiles from independent in situ aircraft observations collected over 2007–2015 (a, c, e, g), showing simulated enhancements or depletion in CO₂ as a result of surface flux (pink hexagons). Additionally, the difference between observations and CT2016 background is plotted (green triangles). We report the vertically resolved and pressure-weighted bias (mean and standard deviation) between simulations and observations (b, d, f, h). Finally, a map shows the location of aircraft profiles (blue circles in panel i) The size of each data point indicates a relative number of samples. Additionally, AirCore profile locations are also shown (orange squares in panel i). AirCore data are used to evaluate biases in boundary conditions in the upper 350 hPa of the column and are presented in Table 2.

Systematic error or bias in $X_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{sim}}$ can arise from errors in the estimation of background and surface flux, both of which are linked by an atmospheric transport model (Fig. 1). We use the same transport model used by Hu et al. (2019) to generate source–receptor relationships and background (Sect. 2.3 and 2.4). The accuracy of the six-member ensemble of fluxes we use is also dependent upon the accuracy of WRF-STILT used in Hu et al. (2019). Thus, potential biases in WRF-STILT form a common thread for error propagation. The combined accuracy of fluxes and transport is evaluated by examining the aircraft vertical profiles of CO₂ collected under NOAA GML's aircraft program (Sweeney et al., 2015) not assimilated by Hu et al. (2019). These data were obtained from NOAA GML's GlobalView Plus version 6.0 data product (Schuldt et al., 2020). We simulate all independent aircraft observations over North America for 2007–2015 (the entire spatiotemporal range of  Hu et al., 2019) using existing WRF-STILT source–receptor relationships. To ensure consistency with Hu et al. (2019), we perform this evaluation with the same background conditions as in that study (i.e., CT2016). Aircraft profiles are sorted in 1 km altitude bins from the surface to 8 km a.s.l. (above sea level) and separated by season. Aircraft vertical profiles of CO₂ (after removing the influence of background) and surface fluxes propagated with WRF-STILT show the net release of CO₂ in non-summer months and net uptake from photosynthesis depleting near-surface CO₂ during the summer (green triangles and pink hexagons in Fig. 6a, c, e, and g, respectively). The difference between independent, unassimilated observations and simulations show that bias at any given altitude level for any season is usually less than 0.5 ppm (Fig. 6b, e, h, j). Bias is also usually the largest near the surface (except in spring). The pressure-weighted partial column (from the surface to 8 km a.s.l.) mean bias (μ_{sim − obs}) ranges from −0.12 in autumn to 0.18 ppm in the spring and is comparable to the typical measurement uncertainty within the in situ Global Greenhouse Gases Research Network of ∼ 0.15 ppm as derived from long-term comparisons of differences between different within-network sampling and analysis approaches for CO₂ (e.g., Andrews et al., 2014; Lan et al., 2017; Sweeney et al., 2015). Low partial column bias relative to independent vertical profile CO₂ data show that errors in WRF-STILT transport contribute very minimally to bias in $X_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{sim}}$.

Table 2 Systematic bias between global 4D CO₂ fields and AirCore profiles above 8 km. Across-model mean and standard deviation is weighted with the pressure weighting function and OCO-2 averaging kernel. All values are in parts per million. The location of AirCore profiles is shown as orange squares in Fig. 6i.

To examine potential systematic errors in the highest part of the column, (above 8 km; i.e., the upper troposphere and the stratosphere), we compare 4D model fields used in this study with CO₂ profiles collected by NOAA GML's AirCore program (Karion et al., 2010; Baier et al., 2021, data version v20200210). There is a considerably larger spread in the models' ability to replicate AirCore observations (Table 2) than surface and aircraft observations (Table 1 and Fig. 6). Bias between model and observations varies seasonally and is highest in the spring and lowest in the summer. However, this region forms the upper ∼ 350 hPa of the column. This has two implications. First, the contribution to total column bias reduces by ∼ 65 %. Second, the ACOS convolution equation (Eq. 1) requires that modeled estimates of $X_{{\mathrm{CO}}_{\mathrm{2}}}$ be appropriately weighted with the averaging kernel (a_i). a_i determines the balance between information obtained by the retrieval and that contained in the OCO-2 prior (${\mathit{\chi }}_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{pri}}$), and tends to decrease from approximately near the surface to <0.5 in the upper third of the column. Thus ∼ 50 % of the error in the highest ∼ 350 hPa of the simulated column is due to errors in ${\mathit{\chi }}_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{pri}}$. Consequently, the total contribution of model bias presented in Table 2 contributes to only ∼50 % of the bias in $X_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{sim}}$ in the uppermost third of the column. The other source of error above 8 km is the error in ${\mathit{\chi }}_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{pri}}$, but this information is currently unavailable in the OCO-2 v10 lite files and is thus ignored in our error analysis. However, errors in ${\mathit{\chi }}_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{pri}}$ are identical for $X_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{sim}}$ and $X_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{ret},\mathrm{bc}}$ and do not affect comparisons between the two.

Combined with the standard deviation across 24 flux–background ensemble members that comprise $X_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{sim}}$, comparisons with independent unassimilated aircraft (Fig. 6) and AirCore profiles (Table 2) allow us to comprehensively assess uncertainties associated with $X_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{sim}}$. We estimate the combined uncertainty from surface flux, background estimation, and transport as follows:

$$\begin{array}{}\text{(5)}& \begin{array}{rl}& {\mathit{\sigma }}_{X_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{sim}}}={\mathit{\mu }}_{\text{model-aircraft}}+{\mathit{\mu }}_{\text{model-AirCore}}\\ & \pm \sqrt{({\mathit{\sigma }}_{\text{model-aircraft}}{)}^{\mathrm{2}}+({\mathit{\sigma }}_{\text{model-AirCore}}{)}^{\mathrm{2}}+({\mathit{\sigma }}_{\mathrm{flux}+\mathrm{bkg}\mathrm{ensembles}}{)}^{\mathrm{2}}},\end{array}\end{array}$$

where

$$\begin{array}{}\text{(6)}& \begin{array}{rl}& {\mathit{\sigma }}_{\mathrm{flux}+\mathrm{bkg}\mathrm{ensembles}}=\\ & \sqrt{{\mathit{\sigma }}_{\mathrm{flux}\mathrm{ensembles}}^{\mathrm{2}}+{\mathit{\sigma }}_{\mathrm{bkg}\mathrm{ensembles}}^{\mathrm{2}}+\mathrm{2}{\mathit{\sigma }}_{\mathrm{flux},\mathrm{bkg}}}.\end{array}\end{array}$$

The first two terms in Eq. (5) represent the pressure-weighted mean partial column bias between modeled and observed CO₂ vertical profiles from aircraft (Fig. 6) and AirCore (Table 2), respectively. The sum of these two terms provides an estimate of bias or systematic error. The third term in Eq. (5) represents unresolved variability.

Table 3 Total uncertainty estimates for $X_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{sim}}$, the mean difference between $X_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{sim}}$ and $X_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{ret},\mathrm{bc}}$, and bias in $X_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{ret},\mathrm{bc}}$. Bias is obtained as the difference between the first two terms in the middle and left columns. All values are in parts per million. The number on the right in the third column indicates the spatial variability in $X_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{ret},\mathrm{bc}}$ bias. Finally, we also present the same quantity for the previous version of OCO-2 retrievals, i.e., the v9 data product. Note that, for v9, only land nadir retrievals are analyzed.

We find that unresolved variability due to model spread in $X_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{sim}}$ is ∼ 0.35 ppm (second term under ${\mathit{\sigma }}_{X_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{sim}}}$ in Table 3), similar across seasons, and significantly lower than the uncertainty of $X_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{ret},\mathrm{bc}}$ as reported in the OCO-2 v10 lite files (∼0.6 ppm). Systematic error or bias in $X_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{sim}}$ (first term under ${\mathit{\sigma }}_{X_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{sim}}}$ in Table 3) shows significant variability across seasons. Comparing this to $\mathrm{\Delta }{X}_{\mathrm{ret},\mathrm{bc}-\mathrm{sim}}$ allows an estimation of mean bias in $X_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{ret},\mathrm{bc}}$ over North America (Table 3). We find that $X_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{ret},\mathrm{bc}}$ bias ranges from −0.62 ppm in autumn to 0.14 ppm in summer. This indicates that $\mathrm{\Delta }{X}_{\mathrm{ret},\mathrm{bc}-\mathrm{sim}}$ in autumn (Figs. 3a and 4c) is entirely due to a residual bias in $X_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{ret},\mathrm{bc}}$ but almost entirely due to bias in ${\mathit{\sigma }}_{X_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{sim}}}$ in the spring (Figs. 3c and 4k). During summer the mean bias in $X_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{sim}}$ is 0.05 ppm, and consequently, the mean bias in $X_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{ret},\mathrm{bc}}$ is 0.25 ppm. However, the mean conceals larger regional differences. In the northeastern quadrant, for instance, $\mathrm{\Delta }{X}_{\mathrm{ret},\mathrm{bc}-\mathrm{sim}}$ (Fig. 4o) of 0.79 ppm translates to a high bias of 0.84 ppm during the summer. Finally, we find that the OCO-2 v10 bias correction shows an improvement over the v9 bias correction, particularly in the winter.

3.4 Evaluating the OCO-2 bias correction

Bias in OCO-2 v10 $X_{{\mathrm{CO}}_{\mathrm{2}}}$ is a combination of footprint bias (eight coincident across track retrievals; C_f) and feature biases (related to surface or atmospheric parameters, e.g., aerosol optical depth; C_p). Finally, a global scaling factor (C₀) obtained from comparisons with TCCON retrievals is used to empirically link retrievals to the WMO scale.

To examine residual feature biases, we perform simple linear regressions between parameters used in the OCO-2 v10 bias correction with $\mathrm{\Delta }{X}_{\mathrm{ret},\mathrm{bc}-\mathrm{sim}}$. These parameters are ΔP_frac (ppm), which accounts for fractional change in $X_{{\mathrm{CO}}_{\mathrm{2}}}$ due to difference in prior and retrieved surface pressure (Kiel et al., 2019), CO₂-grad  del, defined as the difference of the difference between retrieved and prior CO₂ at the surface and at 0.7 times the surface pressure, dws, which is the total retrieved optical depth associated with aerosols from dust, water cloud, and aerosol, and aod_fine, the aerosol optical depth from sulfate and organic carbon (O'Dell et al., 2018; Osterman et al., 2020). Additionally, we perform simple linear regressions with altitude and surface albedo. All parameters are available in the OCO-2 v10 lite files. We find no significant correlations observed between $\mathrm{\Delta }{X}_{\mathrm{ret},\mathrm{bc}-\mathrm{sim}}$ and any parameters suggesting that there are no regional-scale parametric biases over North America for our study period that are not already removed by the OCO-2 global bias correction. These figures are shown in the Supplement.

4 Implications for carbon flux estimation

The impact of recent surface flux on $X_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{sim}}$ is small (e.g., right column in Fig. 4). The interquartile range of ${\mathrm{\Delta }}_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{flux}}$ over the entire spatiotemporal domain is less than 1 ppm, implying that the imprint of recent surface flux on the total column is roughly half the magnitude of the OCO-2 bias correction (blue curves in Fig. 3). Moreover, only around 2 % of simulations in the summer (when surface fluxes are highest) are associated with absolute ${\mathrm{\Delta }}_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{flux}}$ magnitudes higher than 4 ppm, indicating that recent surface flux rarely accounts for more than a ∼ 1 % change in $X_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{sim}}$. In autumn, ${\mathrm{\Delta }}_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{flux}}$ is of the same magnitude as bias in $X_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{ret},\mathrm{bc}}$ (Table 3). OCO-2 retrievals in autumn 2014 are therefore unlikely to provide reliable estimates of North American surface flux. During the summer, OCO-2 has a high bias of 0.25 ppm over the continent, but this bias may be significantly larger in the eastern half of the domain (Fig. 4o). In the northeastern quadrant, for example, a potential bias of 0.84 ppm constitutes ∼ 47 % of the mean ${\mathrm{\Delta }}_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{flux}}$ of −1.79 ppm. Considering that the vast majority of current inverse modeling or data assimilation systems used for CO₂ flux estimation are designed to correct errors in an a priori estimate, the effective flux signal is considerably smaller than shown above (right column in Fig. 4). In fact, when projected onto the total column, we find that the difference between a biospheric prior flux model (CASA-CMS) and flux optimized using in situ observations over the continent by the CarbonTracker–Lagrange inversion system (using the same prior flux estimate; Hu et al., 2019) for January and July 2015 are 0.15±0.38 and 0.09±1.06 ppm, respectively. Flux adjustment impacts on $X_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{sim}}$ are usually indistinguishable from 0 in January. While a gradient of ∼ 3 ppm is visible across the continent in July (Fig. 7), less than a third of simulated retrievals show differences between prior and optimized flux greater than ±1 ppm, when biospheric uptake over North America is strongest. Errors in terrestrial biosphere models of CO₂ flux translate to similarly small impacts on $X_{{\mathrm{CO}}_{\mathrm{2}}}$, posing questions on the utility of these data in evaluating terrestrial biosphere models.

Figure 7 Difference between prior and optimized biospheric flux from CarbonTracker–Lagrange shows small differences when projected onto the total column (ppm).

5 Conclusions

In this study, we compare 1 year of $X_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{ret},\mathrm{bc}}$ over North America from NASA's OCO-2 (v10; land nadir and glint retrievals) satellite against synthetic columns that are constructed using a high-resolution regional model of atmospheric transport and driven by fluxes and background that are optimally consistent with in situ measurements of CO₂ dry air mole fraction, which are rigorously calibrated to the WMO CO₂ X2007 scale. Although $X_{{\mathrm{CO}}_{\mathrm{2}}}$ from OCO-2 and from the posterior of in situ data inversions has been compared previously, and used in its bias correction (O'Dell et al., 2018; Kiel et al., 2019), this is the first such evaluation at the regional scale that uses high-resolution atmospheric transport. We use a suite of optimized non-fire net ecosystem exchange fluxes and background fields to assess under in $X_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{sim}}$. Potential systematic errors in fluxes, transport, and background fields are evaluated by comparisons with vertical gradients of atmospheric CO₂ from independent aircraft and AirCore vertical profiles. $X_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{sim}}$ is associated with errors arising from unresolved variability in model fields and systematic bias. The first of these results in an uncertainty of ∼ 0.35 ppm. Bias or systematic error in $X_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{sim}}$ is found to vary seasonally and ranges from −0.01 ppm in autumn to 0.50 ppm in the spring. Bias is highest in the upper 350 hPa of the column (Table 2), a region that is most poorly constrained by atmospheric measurements. However, the effect of this bias is relatively small in the total column comparisons.

Comparisons with $X_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{sim}}$ show that the OCO-2 v10 global bias correction greatly improves the quality of OCO-2 data over North America (Fig. 3). However, generally good agreement between $X_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{ret},\mathrm{bc}}$ and $X_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{sim}}$ at the continental scale masks significant differences at regional scales and for some seasons (Fig. 4). Error analysis of the components of $X_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{sim}}$ (i.e., transport, background, and flux) allows us to better characterize the difference between simulations and retrievals. Differences in $\mathrm{\Delta }{X}_{\mathrm{ret},\mathrm{bc}-\mathrm{sim}}$ are highest in autumn and indicative of a low bias in $X_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{ret},\mathrm{bc}}$ of 0.62 ppm, which is identical to the mean impact of recent surface flux (mean ${\mathrm{\Delta }}_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{flux}}$ over the continent is 0.64 ppm) in that season. In winter, a low bias in $X_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{ret},\mathrm{bc}}$ of 0.36 ppm is roughly 50 % of the mean ${\mathrm{\Delta }}_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{flux}}$ of 0.71 ppm. In summer, we find spatially coherent regional patterns in $\mathrm{\Delta }{X}_{\mathrm{ret},\mathrm{bc}-\mathrm{sim}}$. $\mathrm{\Delta }{X}_{\mathrm{ret},\mathrm{bc}-\mathrm{sim}}$ is highest in the northeastern quadrant of North America (Fig. 4o) at −0.81 ppm, 50 % of the mean expected ecosystem flux impact over this region. Since inverse models of CO₂ flux usually optimize a prior flux estimate, the surface flux signal (i.e., difference between prior and optimized flux) in $X_{{\mathrm{CO}}_{\mathrm{2}}}$ is minuscule (Fig. 7), significantly smaller than the magnitude of the OCO-2 v10 bias correction, and translates to extremely strenuous requirements on the quality of space-based retrievals. The OCO-2 community has worked diligently to reduce uncertainty on satellite retrievals (e.g., O'Dell et al., 2012, 2018; Kiel et al., 2019; Wunch et al., 2017; Kulawik et al., 2019); bias in v10 retrievals over North America has reduced in autumn and winter of 2014–2015 compared to the v9 data product (Table 3), but further improvement is necessary for both existing satellite data sets and planned missions that will provide this quantity in order to accurately constrain surface fluxes in a changing climate. Finally, we argue that a greatly expanded global reference network of calibrated in situ vertical profile measurements is necessary to reliably detect and correct systematic errors in satellite $X_{{\mathrm{CO}}_{\mathrm{2}}}$.