Characteristics of interannual variability in space-based XCO₂ global observations

Atmospheric carbon dioxide (CO₂) accounts for the largest radiative forcing among anthropogenic greenhouse gases. There is, therefore, a pressing need to understand the rate at which CO₂ accumulates in the atmosphere, including the interannual variations (IAVs) in this rate. IAV in the CO₂ growth rate is a small signal relative to the long-term trend and the mean annual cycle of atmospheric CO₂, and IAV is tied to climatic variations that may provide insights into long-term carbon–climate feedbacks. Observations from the Orbiting Carbon Observatory-2 (OCO-2) mission offer a new opportunity to refine our understanding of atmospheric CO₂ IAV since the satellite can measure over remote terrestrial regions and the open ocean, where traditional in situ CO₂ monitoring is difficult, providing better spatial coverage compared to ground-based monitoring techniques. In this study, we analyze the IAV of column-averaged dry-air CO₂ mole fraction (XCO₂) from OCO-2 between September 2014 and June 2021. The amplitude of the IAV, which is calculated as the standard deviation of the time series, is up to 1.2 ppm over the continents and around 0.4 ppm over the open ocean. Across all latitudes, the OCO-2-detected XCO₂ IAV shows a clear relationship with El Niño–Southern Oscillation (ENSO)-driven variations that originate in the tropics and are transported poleward. Similar, but smoother, zonal patterns of OCO-2 XCO₂ IAV time series compared to ground-based in situ observations and with column observations from the Total Carbon Column Observing Network (TCCON) and the Greenhouse Gases Observing Satellite (GOSAT) show that OCO-2 observations can be used reliably to estimate IAV. Furthermore, the extensive spatial coverage of the OCO-2 satellite data leads to smoother IAV time series than those from other datasets, suggesting that OCO-2 provides new capabilities for revealing small IAV signals despite sources of noise and error that are inherent to remote-sensing datasets.

1 Introduction

Increasing atmospheric CO₂ concentration from anthropogenic emissions is the major driver of the observed warming of Earth's climate since the industrial revolution (IPCC, 2023). Although CO₂ accumulation in the atmosphere generally is ∼45 % of anthropogenic emissions on a multiyear average (Ciais et al., 2013; Friedlingstein et al., 2019), the growth rate shows substantial interannual variability (Conway et al., 1994). The difference between emissions and the atmospheric CO₂ growth rate results from net CO₂ uptake by oceans and terrestrial ecosystems (Prentice et al., 2001; Doney et al., 2009), and the fluctuations reflect variations in the strength of those sinks due to climate variations (Peters et al., 2017; Friedlingstein et al., 2019). Much research has suggested that interannual variability (IAV) in the growth rate is predominantly due to variations in terrestrial ecosystem carbon uptake (Marcolla et al., 2017), even though the average uptake is roughly comparable between land and ocean (Le Quéré et al., 2009). Existing atmospheric CO₂ observations from surface flask sampling and in situ networks have been used to estimate global- and regional-scale interannual variability in CO₂ fluxes (Gurney et al., 2008; Peylin et al., 2013; Keppel-Aleks et al., 2014; Piao et al., 2020). We note, however, that the surface-observing network is located primarily at land and coastal sites, and more subtle ocean-flux signals may be obscured by the large IAV in terrestrial fluxes.

Previous analyses of surface CO₂ IAV have shown a strong relationship with the phase and intensity of El Niño–Southern Oscillation (ENSO) (Le Quéré et al., 2009; Schwalm et al., 2011). ENSO variations originate from coupled ocean–atmosphere dynamics that are reflected in large wind and sea surface temperature anomalies over the central and eastern Pacific Ocean. ENSO affects the climate of much of the tropics and subtropics via atmospheric teleconnections on timescales of 2–7 years (Timmermann et al., 2018). On land, suppressed precipitation and high temperature associated with positive phases of ENSO (El Niño conditions) suppress CO₂ uptake by tropical ecosystems while promoting fires that further reduce the CO₂ uptake by land (Feely et al., 2002; McKinley et al., 2004; Piao et al., 2009; Wang et al., 2014). Although of smaller magnitude, the equatorial Pacific Ocean experiences weakening of the easterly trade winds and suppression of ventilation of deep, cold, carbon-rich waters to the surface during an El Niño, reducing the efflux of natural CO₂ to the atmosphere (Patra et al., 2005; Chatterjee et al., 2017).

Chatterjee et al. (2017) were able to directly observe the ocean-flux-driven signal on atmospheric CO₂ from El Niño for the first time using XCO₂ (column-averaged dry-air CO₂ mole fraction) observed over the ocean by NASA's OCO-2 satellite. Space-based observations from OCO-2, which was launched in July 2014, provide novel opportunities to characterize the patterns of IAV in XCO₂ in areas that were previously not directly observed by existing monitoring networks. The IAV in XCO₂ is being used implicitly for flux attribution in inverse modeling studies (Nassar et al., 2011). These exciting results, however, must be tempered by an awareness that atmospheric CO₂ IAV is a relatively small signal. For example, IAV in the surface network is about 1 ppm in scale compared to a seasonal amplitude of around 10 ppm at northern high latitudes. OCO-2 measures column-averaged CO₂, so its measurements are sensitive to variations in the boundary-layer mole fraction, which is in direct contact with the land or atmospheric fluxes but also variations in the free troposphere and stratosphere, where flux signals are generally smaller than those observed at the surface (Olsen and Randerson, 2004). Furthermore, variations in the free troposphere are expected to have relatively long correlation length scales due to efficient mixing, making it important to consider the spatial scales at which XCO₂ observations provide unique information. This is especially important in light of analysis which suggests that the error variance budget in OCO-2 observations is large and contains a substantial spatially coherent signal (Baker et al., 2022; Torres et al., 2019; Mitchell et al., 2023).

In this paper, we analyze XCO₂ from OCO-2 to characterize spatiotemporal patterns in IAV at a near-global scale, over both land and ocean, and relate XCO₂ variations to ENSO conditions. We contextualize the information contained in OCO-2 observations by comparing them with space-based GOSAT and ground-based Total Carbon Column Observing Network (TCCON) XCO₂ and with surface measurements of CO₂. Finally, we use these comparisons to emphasize the spatial scales at which the IAV signal emerges from instrumental noise.

2 Data and methods

2.1 Datasets

2.1.1 OCO-2 observatory

We analyzed IAV in dry-air column-averaged mole-fraction XCO₂ inferred from OCO-2 satellite observations. The OCO-2 observatory was launched in July 2014 and has measured passive, reflected solar near-infrared CO₂ and O₂ absorption spectra using grating spectrometers since September 2014 (Eldering et al., 2017). XCO₂ data are retrieved from the measured spectra using the Atmospheric CO₂ Observations from Space (ACOS) optimal-estimation algorithm, which is a full-physics algorithm that takes into account XCO₂ and other physical parameters, including surface pressure, surface albedo, temperature, and water vapor profile in its state vector (O'Dell et al., 2018). The satellite flies in a polar and sun-synchronous orbit that repeats every 16 d with three different observing modes of OCO-2, namely, nadir (land only, views the ground directly below the spacecraft), glint (over ocean and land, views just off the peak of the specularly reflected sunlight), and target (typically for comparison with specific ground-based or airborne measurements) (Crisp et al., 2012, 2017). We use the version 10 OCO-2 level-2 bias-corrected XCO₂ data product from the Goddard Earth Sciences Data and Information Services Center (GES DISC) archive (https://disc.gsfc.nasa.gov/datasets/OCO2_L2_Lite_FP_10r/summary, last access: 1 March 2022), which has been validated with collocated ground-based measurements from TCCON, discussed in more detail in Sect. 2.2. After filtering and bias correction, the OCO-2 XCO₂ retrievals agree well with TCCON in the nadir, glint, and target observation modes and generally have absolute median differences of less than 0.4 ppm and root-mean-square differences of less than 1.5 ppm (O'Dell et al., 2018; Wunch et al., 2017).

2.1.2 TCCON

We corroborate patterns of XCO₂ IAV from OCO-2 with those from TCCON, a ground-based network of Fourier-transform spectrometers (FTSs) that measure direct solar-absorption spectra in the near infrared (Wunch et al., 2011). Retrievals of XCO₂ and other gases are computed using the GGG2014 version of the TCCON standard retrieval algorithm (Wunch et al., 2015), a nonlinear least-squares spectral-fitting algorithm. The TCCON retrievals are tied to the World Meteorological Organization (WMO) X2007 CO₂ scale via calibration with aircraft and AirCore profiles above the TCCON sites (Karion et al., 2010; Wunch et al., 2010). This ensures an accuracy and precision of ∼0.6 ppm (1-sigma) throughout the network (Washenfelder et al., 2006; Messerschmidt et al., 2010; Deutscher et al., 2010; Wunch et al., 2010). TCCON has been used widely as a validation standard by providing independent measurements to compare with multiple satellite XCO₂ retrievals, including OCO-2. In previous work, Sussmann and Rettinger (2020) demonstrated a concept to retrieve annual growth rates of XCO₂ from TCCON data, which are regionally to hemispherically representative in spite of the nonuniform sampling in time and space inherent to the ground-based network. In our study, we focus on IAV in the XCO₂ time series from 26 TCCON sites (Table 1 and Fig. 1) that have at least 3 years of observational coverage within the period from September 2014 to June 2021. These TCCON data have been filtered using the standard filter that is based on a measure of cloudiness and that limits the solar-zenith angle. Data are publicly available from the TCCON GGG2014 data archive (https://tccondata.org/, last access: 2 December 2020) hosted by the California Institute of Technology.

Table 1 TCCON column-averaged dry-air mole fractions of CO₂ (GGG2014 data).

Figure 1 Map showing the locations and the abbreviations of the TCCON sites.

2.1.3 Marine-boundary-layer observations

To explore differences in surface and column-averaged CO₂ IAV, we analyze IAV in the surface CO₂ mole fraction at marine-boundary-layer (MBL) sites in the NOAA (National Oceanic and Atmospheric Administration) cooperative sampling network (https://gml.noaa.gov/dv/site/?program=ccgg, last access: 6 May 2023). At these sites, boundary-layer CO₂ is measured using weekly flask samples (Masarie and Tans, 1995; Dlugokencky et al., 2021). MBL sites are typically far away from anthropogenic sources and regions of active terrestrial exchange, so they provide an estimate for large-scale patterns in the global background CO₂ concentration. The surface MBL dry-air mole-fraction data have an accuracy level of about 0.1 ppm. In this study, we select 16 sites with at least 80 % data coverage for the approximately 7-year period overlapping with OCO-2 (Table 2 and Fig. 2), and the data are aggregated into four north–south zones for comparison with OCO-2 XCO₂: Northern Hemisphere and Southern Hemisphere tropical zones (0–20^∘) and Northern Hemisphere/Southern Hemisphere extratropical zones (20–60^∘). Each belt contains at least three MBL sites. Higher latitudes (60–90^∘) are not considered in this comparison due to the gaps remaining in the OCO-2 XCO₂ record at high latitudes during wintertime and shouldering seasons.

Table 2 Marine-boundary-layer stations within the NOAA Earth System Research Laboratory CO₂ sampling network.

Figure 2 Map showing the locations and the abbreviations of the marine-boundary-layer stations within the NOAA Earth System Research Laboratory CO₂ sampling network.

2.1.4 GOSAT

We compare patterns of XCO₂ IAV from OCO-2 with those from GOSAT. Also known as Ibuki, GOSAT is the world's first satellite dedicated to greenhouse gas monitoring, measuring global total column CO₂ and CH₄ since 2009 with the Thermal and Near infrared Sensor for carbon Observation (TANSO) FTS on board for greenhouse gas monitoring using three SWIR bands and one TIR band (Cogan et al., 2012; Yoshida et al., 2013). Column-averaged dry mole fractions are obtained at a circular footprint of approximately 10.5 km. GOSAT has a regional bias of approximately 0.3 and 1.7 ppm single observation error versus TCCON (Kulawik et al., 2016). We utilize the FTS SWIR level-3 data global monthly 2.5^∘ resolution mean CO₂ mixing ratio products from 2009 June to 2021 December to generate IAV and make comparisons with OCO-2. Level-3 products are generated by interpolating, extrapolating, and smoothing the FTS SWIR column-averaged mixing ratios of CO₂ and applying the geostatistical calculation technique kriging method. GOSAT observation datasets are available to the public at the NIES GOSAT website (https://www.gosat.nies.go.jp/en/about_5_products.html, last access: 5 August 2022).

2.1.5 Multivariate ENSO index (MEI)

We use the bimonthly MEI (downloaded from the Physical Sciences Laboratory: https://psl.noaa.gov/enso/mei/, last access: 5 April 2023) to explore the relationship between CO₂ IAV and ENSO. The MEI is the time series of the leading combined empirical orthogonal function of five different variables (sea level pressure, sea surface temperature, zonal and meridional components of the surface wind, and outgoing longwave radiation) over the tropical Pacific basin. Positive values in the MEI indicate El Niño conditions, while negative values indicate La Niña conditions, and the magnitude reflects the relative strength. Unlike other ENSO indices which use only one climate metric (e.g., the sea level pressure difference between Tahiti and Darwin or the sea surface temperature anomaly within a predefined box), the MEI provides for a more complete and flexible description of the ENSO phenomenon than traditional single-variable ENSO indices and has less vulnerability to errors (Wolter and Timlin, 2011).

2.2 Methods

2.2.1 Spatial aggregation

We aggregate daily XCO₂ observations from the version 10 OCO-2 level-2 lite product to the monthly scale, exploring patterns of IAVs at three spatial scales: grid-cell level, zonal averages over 5^∘ of latitude, and broad zonal belts. Aggregating soundings reduces random noise in the observations, mitigates the impact of data gaps due to cloud cover, and partly mitigates effects from low winter sunlight levels in polar regions. For grid-cell-level analysis, we aggregate data equatorward of 45^∘ to $\mathrm{5}{}^{\circ }\times \mathrm{5}{}^{\circ }$ bins since these data are not limited by polar night or degraded by high solar-zenith angles during winter. Poleward of 45^∘ in both hemispheres, we aggregate the satellite observation to a latitude–longitude resolution of $\mathrm{5}{}^{\circ }\times \mathrm{10}{}^{\circ }$ to compensate for fewer and noisier soundings in these latitudes, especially during winter and its shoulder seasons. Within each $\mathrm{5}{}^{\circ }\times \mathrm{5}{}^{\circ }$ or $\mathrm{5}{}^{\circ }\times \mathrm{10}{}^{\circ }$ grid cell, only months that have more than five soundings are included in the analysis. Our criteria for aggregation are based on sensitivity experiments in which we modulated the grid-cell resolution from $\mathrm{1}{}^{\circ }\times \mathrm{1}{}^{\circ }$ to $\mathrm{15}{}^{\circ }\times \mathrm{15}{}^{\circ }$ (Figs. S1 and S2 in the Supplement) and varied the threshold on the required number of soundings within a month from 1 to 25 (Figs. S3–S5 in the Supplement). Our goal was to reduce noise but maintain high spatial coverage (Figs. S6 and S7 in the Supplement). The $\mathrm{5}{}^{\circ }\times \mathrm{5}{}^{\circ }$ and $\mathrm{5}{}^{\circ }\times \mathrm{10}{}^{\circ }$ aggregations strike the necessary balance of reducing noise (evidenced by the smoother IAV amplitude fields as aggregation increases in Fig. S1) but maintaining spatial information by not oversmoothing (evidenced by the fact that the aggregation occurs at spatial scales finer than the “elbow” where correlations among 1^∘ grid cells stop changing with separation distance in Fig. S8 in the Supplement)

In our analysis, we also aggregate data to zonal averages. At intermediate spatial scales, we average all data around the 5^∘ latitude bins described above. For comparison with TCCON and MBL data, which are spatially sparse, we further aggregate XCO₂ data into four broad zonal belts – each of which contains at least one TCCON or three MBL stations – (delineated in Tables 1 and 2) to assess IAV patterns among the datasets. Keppel-Aleks et al. (2014) showed that drivers of IAV (i.e., temperature, drought stress, or fire) could be attributed when surface CO₂ were aggregated into similar broad zonal belts, whereas process-level attribution was not possible with global averaging. We therefore analyze broad zonal belts to gain a large-scale understanding of how three CO₂ datasets are similar and where differences lie.

2.2.2 Deriving interannual variations

We use a consistent process to calculate IAV (Eq. 1) from the raw OCO-2, TCCON and MBL time series. The methodology is based on approaches used in Keppel-Aleks et al. (2014) and NOAA curve fitting methodology (Thoning et al., 1989). We decompose the raw time-series data into a long-term trend (which is a function of location (x,y) and time (t)), a seasonal cycle (which is a function of location and calendar month (m)), and IAV anomalies using Eq. (1):

$$\begin{array}{}\text{(1)}& \begin{array}{rl}\text{IAV}(x,y,t)=& \text{raw}(x,y,t)-\text{trend}(x,y,t)\\ & -\text{seasonal}(x,y,m).\end{array}\end{array}$$

We first fit a third-order polynomial to the raw time series to calculate the observed trend at each location (Fig. 3a). After removing the trend calculated at each grid cell (Fig. 3b), we calculate a mean seasonal cycle by taking the mean value of all January, February, etc., data (Fig. 3c). Particularly at high latitudes, some months are systematically undersampled. For these grid cells, we must have at least 2 years with sufficient observations to calculate a climatological mean for that month; otherwise, that calendar month is assumed to have insufficient data to infer the IAV. Finally, we remove the mean seasonal cycle from the detrended time series at each grid cell to obtain the IAV anomaly time series (Fig. 3d). Given the short data record, we quantify the uncertainty in our calculation of the climatological seasonal cycle as the standard error for each calendar month (blue shading in Fig. 3c), and this uncertainty is propagated to the corresponding IAV time series (Fig. 3d). We fit a third-order polynomial to the raw time series since the GOSAT, MBL, and TCCON time series extend over a decade in length. We confirm that the use of a third-order polynomial versus a second-order polynomial does not remove the IAV signal from the shorter OCO-2 time series (Fig. S9 in the Supplement).

Figure 3 Methodology to calculate the CO₂ interannual variability time series using OCO-2 XCO₂ data at the 5^∘ grid cell at 20^∘ N, 155^∘ W, which contains Moana Loa as an example. (a) 5^∘-resolution monthly mean raw OCO-2 XCO₂ and the associated third-order polynomial trend. (b) Detrended monthly XCO₂ after removing the long-term trend with a repeating 12-month annual cycle obtained from calculating the mean for each month. The light-blue shading gives the uncertainty of the seasonal cycle, which is derived by calculating the standard deviation across all Januaries, Februaries, etc. (c) 12-month mean annual cycle together with the uncertainty range plotted in panel (b). (d) Resulting interannual variability when the mean annual cycle is removed from detrended time series.

3 Results

3.1 Spatiotemporal variations based on OCO-2 observation

When averaged into broad zonal belts representing the tropics and mid-latitudes, the OCO-2 XCO₂ IAV time-series anomalies range between −0.5 and 0.75 ppm (Fig. 4a). All latitude bands show increasing IAV during positive MEI (El Niño) and decreasing IAV during negative MEI (La Niña), although the phasing varies among latitudes. During the strong 2015–2016 El Niño, which began around March 2015 and reached its peak at the start of 2016, XCO₂ showed the largest IAV. The Southern Hemisphere extratropical region (Fig. 4d) has a larger and more rapid response in the IAV associated with ENSO compared to other zones, especially for the smaller El Niño that peaked at the beginning of 2020. At this time, the XCO₂ IAV time series (Fig. 4d) had an anomaly nearly twice as large as that of other latitude belts (Fig. 4a–c). During both El Niño events, the IAV time series in the Northern Hemisphere tropics zone peaks nearly 6 months after the maximum MEI value.

Figure 4 IAV time series averaged for zonal bands between 60^∘ N and 60^∘ S from four different observing strategies: space-based OCO-2 XCO₂ (black), surface CO₂ observations from NOAA's marine-boundary-layer (MBL) sites (blue), ground-based TCCON XCO₂ (red), and space-based GOSAT XCO₂ (grey). (a) Temperate Northern Hemisphere (20–60^∘ N), (b) tropical Northern Hemisphere (0–20^∘ N), (c) tropical Southern Hemisphere (0–20^∘ S), (d) temperate Southern Hemisphere (20–60^∘ S). For all panels, the background shading indicates the multivariate ENSO index (MEI), which is positive during El Niño phases.

We assess the spatial correlation patterns with no time lag, 3-month lag, and 6-month lag between the IAV time series and the MEI (Fig. 8a). The XCO₂ IAV time series have a strong correlation coefficient with the MEI at both the Southern Hemisphere and Northern Hemisphere low latitudes from 0 to 30^∘ N at lag 0, whereas in the Northern Hemisphere extratropics, the maximum positive correlation occurs at month 4 (Fig. 8b). The positive correlation between the MEI and the IAV time series is gradually attenuated, with no clear correlation at 6 months' lag (Fig. 8c).

The differences in temporal phasing between the broad zonal belts (Fig. 4a) associated with El Niño events can be linked to transport of El Niño-driven CO₂ flux anomalies away from the tropics when zonal means are calculated from OCO-2 observations at 5^∘ latitude resolution (Fig. 5). For the two El Niño periods in 2015–2017 and late 2018 to 2021, high IAV values originate in the tropics, and a smooth transition to high IAV values is seen at higher latitudes as time progresses (Fig. 5a). We note that fluxes outside the tropics may also be influenced by ENSO-related climate variability, yet the transport of tropical-driven anomalies appears to dominate. This 7-year study period also captures the half-year lags for atmospheric transport or climate–ecological teleconnections that impact XCO₂ variations in the far north. While the OCO-2 patterns largely conform to the variability expected based on ENSO and are in broad agreement with other observational networks, there are some anomalies that cannot be explained, such as the high XCO₂ in early 2020 around 60^∘ S (Fig. 5a). Even with more aggressive data filtering, this episode persists, requiring more investigation of unknown geophysical drivers of high XCO₂ or potential retrieval issues that could cause a high bias.

Figure 5 Hovmöller diagram showing zonal-mean OCO-2 XCO₂ IAV time series for 5^∘ latitude bins (a) and the zonal standard deviation of XCO₂ IAV (b), which gives an estimate of coherence in the IAV patterns among grid cells in the 5^∘ zonal belt.

We quantify coherence in CO₂ IAV within a latitude circle by taking the standard deviation across grid-cell-level IAV anomalies within each 5^∘ latitude zone. The standard deviation among grid cells is highest in the far north, with values as high as 1 ppm poleward of 45^∘ N and as low as 0.2 ppm in the southern tropical bands (Fig. 5b), indicating that IAV is less spatially coherent in the Northern Hemisphere. This may be consistent with studies that show greater IAV in terrestrial ecosystem fluxes (concentrated in the Northern Hemisphere) (Zeng et al., 2005) relative to ocean fluxes or may reflect the fact that our IAV time series also retains the imprint of sampling, measurement, and retrieval errors, which become more pronounced at higher latitudes. In general, there is no time-dependent or ENSO-related pattern for the longitudinal variation of IAVs (no obvious changes during the two El Niño periods), which suggests that the variation within each 5^∘ band may be approximately stable and does not change substantially with interannual climate events.

The XCO₂ IAV amplitude (the standard deviation of the IAV time series) is notably larger over continental grid cells compared to ocean grid cells (Fig. 6). In both hemispheres, the IAV amplitude over subtropical ocean basins is less than 0.4 ppm, while the IAV amplitude over tropical land in Southeast Asia, the Congo forests, and the Amazon basin is about 1 ppm. At higher latitudes, the XCO₂ IAV amplitude can exceed 1.2 ppm above deciduous and boreal forests in North America and Eurasia. Higher values over land likely occur due to the active CO₂ exchange between the terrestrial ecosystem and the atmosphere, but we cannot rule out that retrievals over land show more variance due to complex topography, albedo, etc., which are elements of the retrieval state vector. Nevertheless, over land areas with low carbon exchange (e.g., Australia, the Middle East, the Sahara), the XCO₂ IAV amplitude is nearly at the same low level as the ocean basins. It is worth noting that, for high-latitude regions, including both northern continents and the Southern Ocean, OCO-2 does not obtain observations over a full calendar year (stippled grid cells in Fig. 6) due to polar nights, low light levels, and high solar-zenith angles. The XCO₂ IAV amplitudes are less zonally coherent through these regions than those in the tropics and at the mid-latitudes for both land and ocean. When averaging all ocean or land grid cells around a latitude circle, the zonal-mean IAV amplitude over the ocean ranges from 0.3 to 1.0 ppm, while the land IAV amplitude ranges from 0.4 to 1.1 ppm (Fig. 9). Both the land and ocean profiles have similar north–south patterns, with a higher IAV amplitude in the Northern Hemisphere and a lower IAV amplitude in the Southern Hemisphere and small IAV amplitudes in the subtropics of both hemispheres, with more scatter among land grid cells than the ocean (Figs. 5b and 9), suggesting either the influence of local-flux IAV on land or greater error associated with retrievals on land. We note better coherence between the XCO₂ IAV time series of each local grid cell and that of zonal-mean XCO₂ IAV time series for the ocean, with correlation coefficients of approximately 0.8. In contrast, land grid cells are generally correlated with the zonal mean at around 0.4 to 0.6 (Fig. 10).

Figure 6 OCO-2 XCO₂ IAV amplitude, determined as the standard deviation of the IAV time series. Data equatorward of 45^∘ are averaged at 5^∘ by 5^∘ resolution, and data poleward of 45^∘ are averaged at 5^∘ by 10^∘ resolution. Shaded regions indicate grid cells that lack mean annual cycle data for at least 2 calendar months due to polar night or related retrieval challenges.

3.2 OCO-2 XCO₂ IAV compared to GOSAT XCO₂ IAV

We carried out comparisons between the global spatiotemporal pattern of XCO₂ IAV between OCO-2 and GOSAT, since GOSAT has data beginning in 2009. The XCO₂ time series from OCO-2 provides higher coverage over midlatitude oceans and tropical rainforests (stippling in Figs. 6 and 7). The IAV amplitude of OCO-2 is generally smaller than that of GOSAT worldwide (Figs. 6 and 7), which may be due to greater data volume and reduced noise in the OCO-2 dataset (Wu et al., 2020). OCO-2 and GOSAT zonal-mean IAV time series generally share the same feature from 2014 to 2021 (Fig. 4a–d), with an increasing trend during El Niño and a decreasing trend during La Niña; however, the GOSAT XCO₂ shows a delayed response at the northern midlatitudes, by almost 9 months, to the strong 2015 El Niño compared to the other datasets. Generally, GOSAT IAV time series are nosier, from month to month, compared to those from OCO-2.

Figure 7 Similar to Fig. 6 but based on GOSAT data.

Figure 8 Correlation coefficient between local grid-cell OCO-2 XCO₂ IAV time series and MEI (a) for synchronous time series, (b) with 3-month lags, and (c) with 6-month lags.

Figure 9 Latitudinal profile for the zonal mean of IAV amplitude and the standard deviation among land (green) or ocean (blue) grid cells in each latitude band (shaded area). Individual points represent all grid cells valid in the IAV record within a certain zonal band.

Figure 10 Correlation coefficient between local grid-cell IAV time series and the corresponding 5^∘ zonal-mean OCO-2 XCO₂ IAV time series.

3.3 XCO₂ IAV compared to surface and TCCON ground-based sites

Given that the small IAV signal (up to 1 ppm over land and smaller over the ocean) is similar in magnitude to noise and systematic bias in OCO-2 soundings (Torres et al., 2019), we corroborate patterns of IAV from OCO-2 with other datasets. The OCO-2 IAV time series in broad latitudinal belts share similarities with those of TCCON XCO₂ and MBL surface CO₂ ground-based IAV time series, with all time series showing similar relationships with the MEI. Especially striking is that all time series capture the lagged response in the Northern Hemisphere midlatitude belt to the strong 2015/16 El Niño (Fig. 4a–d). Although the patterns are similar, the magnitude of IAV at the MBL sites is almost double the IAV in the OCO-2 XCO₂ time series. Given that the atmospheric boundary layer, where surface observations are made, is on average 10 % of the total column, this suggests that much IAV in total column observations is present within the free troposphere. For TCCON, the amplitude of IAV is similar to that of OCO-2, since both methods capture total column variations. We note that the zonal IAV time series for MBL and TCCON appear to have more high-frequency variations than those from OCO-2 (Figs. S10–S12 in the Supplement), which likely stems from the fact that the zonal composites are developed from sparse ground-based sites (between 1 and 12 observatories) within each latitude belt, whereas the satellite measures at all longitudes within a belt though with more limited time resolution. The zonal-mean OCO-2 observations are correlated with MBL sites within the same latitude band with R between 0.5 and 0.75 (diagonal elements in Fig. 13b). Correlations between zonal TCCON and OCO-2 observations range between 0.15 and 0.55 (Table S1 in the Supplement). The correlations are weakest in the northern tropical band, where TCCON data were unavailable during the strong El Niño (Fig. 3c). It is noteworthy that OCO-2 zonal averages are more correlated among different latitudes than are MBL or TCCON observations (off-diagonal elements in Fig. 13c–e). The greater correlation across latitudes for OCO-2 compared to MBL sites is likely due to the sensitivity of the OCO-2 XCO₂ observations to the free troposphere, where meridional transport is more rapid than at the surface. While TCCON data are also sensitive to the free troposphere, we hypothesize that the zonal-belt averages for TCCON, constructed from only a few sites, are more affected by noise, both instrumental and geophysical, and thus show lower coherence than the OCO-2 XCO₂ averages constructed from the whole latitudinal bands.

We further compared the IAV from OCO-2 XCO₂ with TCCON stations at the site level (Fig. 12). Across all the sites, the IAV amplitude generally shows good agreement and lies between 0.4 and 1.2 ppm. We note a slight low-IAV amplitude in OCO-2 relative to TCCON for all five sites in the Southern Hemisphere which lie below the one-to-one line. Low OCO-2 IAV amplitudes may be due to the fact that a $\mathrm{5}{}^{\circ }\times \mathrm{5}{}^{\circ }$ grid cell encompassing these near-coastal locations includes both land and ocean OCO-2 soundings and may be due to specific sources of variance from retrieval bias affected by surface type for OCO-2 (e.g., Fig. 9). It is also worth noting that OCO-2 looks at a region of 5^∘ by 5^∘ grid cells (or 5^∘ by 10^∘ at higher latitudes) around TCCON sites, so there are different signals affecting the variance between the two types of observations.

We derive the regression slopes and correlation coefficient R between OCO-2 and monthly averaged TCCON IAV through bootstrapping linear-regression fitting techniques to investigate the coherence between IAV signals from space-based and in situ ground-based observations. We compute the linear regression 1000 times by iteratively resampling the IAV time series with replacement and calculate the 95 % significance level for regression slopes based on the histogram of the sample distributions during the bootstrapping (Fig. S13 in the Supplement). Despite having similar IAV amplitudes, the IAV time series from OCO-2 are only moderately correlated with those from TCCON (Fig. 11). The regression slopes range from 0.1 to 0.6, and R values are generally around 0.1–0.5, indicating that less than 25 % of the IAV in OCO-2 is explained by IAV measured by TCCON. These R values are, as expected, smaller than the zonal averages shown in Fig. 11b, which average some of the site-level noise for TCCON and grid-cell-level noise for OCO-2. The detailed XCO₂ IAV time series of each site (Fig. S10) for OCO-2 and TCCON show that the IAV time series in the Northern Hemisphere are more variable, which can partly explain the hemispheric difference in amplitude, slope, and correlation coefficients.

Figure 11 Latitudinal profile of regression slope (panel a) and correlation coefficient (R, panel b) of OCO-2 versus TCCON XCO₂ IAV. The slope and R values are based on using monthly XCO₂ IAV. The error bars result from a Monte Carlo bootstrapping approach. The colors represent the number of months of data which are used for the regression calculation given gaps in both the OCO-2 and TCCON datasets.

Figure 12 Comparison of OCO-2 and TCCON XCO₂ IAV amplitude at individual sites. Colors reflect site latitudes. The grey dashed line is the one-to-one identity line. The grey solid line is the error bar of the IAV amplitude.

Figure 13 Correlation coefficient (R) between mean CO₂ time series using three observing strategies. Panel (a) shows the correlation between zonal-mean OCO-2 XCO₂ IAV and zonal-mean marine-boundary-layer CO₂. Panel (b) shows the correlation between zonal-mean XCO₂ IAV from OCO-2 and TCCON. Panels (c–e) show the correlation in zonal-mean IAV time series across four latitude bands for a single observing strategy. Panel (c) shows OCO-2 XCO₂, panel (d) shows MBL CO₂, and panel (e) shows TCCON XCO₂. For panels (c–e), the diagonal elements are 1 by construction. Zonal bands include the tropical (0–20^∘) and Northern Hemisphere/Southern Hemisphere temperate (20–60^∘) zones.

4 Discussion

We use 7 years of OCO-2 total column carbon dioxide observations from late 2014 to mid 2021 to illustrate the global temporal–spatial patterns of atmospheric XCO₂ interannual variations. OCO-2 and GOSAT showed reasonable agreement (Fig. 4) in Northern Hemisphere and Southern Hemisphere tropical zones (0–20^∘), although there were some notable phase differences during the strong 2015 El Niño for GOSAT compared to the other time series in both the northern and southern extratropical regions. In contrast, OCO-2 shows good temporal agreement with the ground-based observations from the MBL and TCCON. The temporal agreement of the OCO-2 and TCCON XCO₂ IAV time series and the MBL surface CO₂ IAV time series in broad zonal belts improves our confidence that we can quantify IAV time series from the satellite record. We note that amplitude differences remain among the time series, owing to two major factors: first, compared to MBL surface observations, we expect XCO₂ time series to have smaller amplitudes of variability since it integrates over the entire atmospheric column (Olsen and Randerson, 2004), and second, the fact that the OCO-2 time-series averages around a full-latitude circle rather than a few discrete sites reduce some of the IAV contained in site-level records. From the space-based and ground-based detection, we are able to characterize the global response of OCO-2 and TCCON XCO₂ or MBL surface CO₂ IAV to ENSO and track the CO₂ IAV against the positive/negative phase of ENSO, together with the transport of the signal from south to north (Fig. 4). All the datasets show consistent patterns in the response to the El Niño periods, although we note that the IAV amplitude is a factor of almost 2 smaller in the column-averaged mole fraction compared to the boundary-layer CO₂, which reflects the fact that IAV variations emerge due to surface fluxes in the lower part of the atmosphere (Olsen and Randerson, 2004) but are efficiently transported into the free troposphere, which comprises the bulk of the column. When taken together, the use of surface and column data may allow better separation of transport-driven versus local flux-driven variations at the interannual timescale. In the future, as partial column retrievals (e.g., Kulawik et al., 2017) mature, intercomparisons of the lowermost tropospheric partial columns may provide a useful bridge between variations in surface MBL observations and total column observations.

Our results, however, underscore the difficulty in detecting IAV signals from remote sensing of XCO₂ – while Northern Hemisphere seasonal amplitudes are typically 10 ppm in scale (Basu et al., 2011), the magnitude of OCO-2-detected XCO₂ IAV is almost an order of magnitude smaller (less than 0.4 ppm over the ocean and about 1 ppm over continents). The magnitude of IAV is therefore comparable to other components of the XCO₂ variance budget; for instance, Torres et al. (2019) show random noise in individual OCO-2 soundings of about 0.3 ppm in the Southern Hemisphere and about 0.7 ppm in the Northern Hemisphere and spatially coherent errors in the retrievals ranging from 0.3 to 0.8 ppm (Torres et al., 2019). Moreover, the uncertainty which originally comes from the varying climatological seasonal cycle can also reach a level of 0.5 ppm (Fig. 3d). Therefore, robust partitioning of IAV from the observed XCO₂ signal at a given location requires a comprehensive variance budget (Mitchell et al., 2023), and efforts to infer interannual variations in fluxes from OCO-2 must take grid-cell-level variance into account or leverage zonally averaged data, which are characterized by greater separation between IAV signal and noise.

Our analysis shows that proper spatial averaging of the monthly XCO₂ signal can mitigate the imprint of random noise and systematic effects from weather systems at submonthly timescales. Based on sensitivity tests, we recommend averaging low- to mid-latitude XCO₂ (equatorward of 45^∘) to $\mathrm{5}{}^{\circ }\times \mathrm{5}{}^{\circ }$ bins and a $\mathrm{5}{}^{\circ }\text{ latitude}\times \mathrm{10}{}^{\circ }\text{ longitude}$ grid cell poleward of 45^∘, ensuring that each grid-cell aggregates at least five soundings within a month. At these levels of spatial averaging, the XCO₂ IAV amplitude was comparable to that of the co-located ground-based XCO₂ IAV amplitude measured by TCCON (Fig. 12). However, the moderate to low correlation between the IAV time series from each monitoring platform reveals the discrepancies of the two measurements in sampling, detection, or retrieval, suggesting that one or both is still convolving another source of variance with the calculated IAV signal. Based on the good agreement between the two time series in broad zonal belts, we expect that random noise in both observations may degrade the comparison.

The smaller coherence in the IAV time series in nearby land and ocean grid cells may be due to larger error over land or may reflect the fact that XCO₂ observations over land contain information about heterogeneous local-flux IAV. Complete analysis of the variance budget for OCO-2 observations (Mitchell et al., 2023) will elucidate the likely imprint of each process. When using IAV time series for flux inference, it will be crucial to account for non-flux imprints such as imprint from atmospheric transport, random errors, systematic errors, and remote geophysical coherence on the time series (e.g., Torres et al., 2019; Mitchell et al., 2023), since spurious attribution of IAV will lead to biased fluxes.

5 Conclusions

We examined IAV in OCO-2 data to determine whether the small variations that result from interannual flux variations can be detected in light of other sources of variance in the space-based dataset. Our results show that zonal averages reveal relationships with ENSO that are consistent with those from an established ground-based monitoring network. Zonal averages greatly reduce random noise in XCO₂ compared to $\mathrm{5}{}^{\circ }\times \mathrm{5}{}^{\circ }$ averages. In general, OCO-2 can successfully monitor CO₂ IAV over both land and ocean, contributing important spatial coverage beyond inferences of IAV from existing ground-based networks.