A global surface CO₂ flux dataset (2015–2022) inferred from OCO-2 retrievals using the GONGGA inversion system

Accurate assessment of the size and distribution of carbon dioxide (CO₂) sources and sinks is important for efforts to understand the carbon cycle and support policy decisions regarding climate mitigation actions. Satellite retrievals of the column-averaged dry-air mole fractions of CO₂ (XCO₂) have been widely used to infer spatial and temporal variations in carbon fluxes through atmospheric inversion techniques. In this study, we present a global spatially resolved terrestrial and ocean carbon flux dataset for 2015–2022. The dataset was generated by the Global ObservatioN-based system for monitoring Greenhouse GAses (GONGGA) atmospheric inversion system through the assimilation of Orbiting Carbon Observatory-2 (OCO-2) XCO₂ retrievals. We describe the carbon budget, interannual variability, and seasonal cycle for the global scale and a set of TransCom regions. The 8-year mean net biosphere exchange and ocean carbon fluxes were −2.22 ± 0.75 and −2.32 ± 0.18 Pg C yr⁻¹, absorbing approximately 23 % and 24 % of contemporary fossil fuel CO₂ emissions, respectively. The annual mean global atmospheric CO₂ growth rate was 5.17 ± 0.68 Pg C yr⁻¹, which is consistent with the National Oceanic and Atmospheric Administration (NOAA) measurement (5.24 ± 0.59 Pg C yr⁻¹). Europe has the largest terrestrial sink among the 11 TransCom land regions, followed by Boreal Asia and Temperate Asia. The dataset was evaluated by comparing posterior CO₂ simulations with Total Carbon Column Observing Network (TCCON) retrievals as well as Observation Package (ObsPack) surface flask observations and aircraft observations. Compared with CO₂ simulations using the unoptimized fluxes, the bias and root mean square error (RMSE) in posterior CO₂ simulations were largely reduced across the full range of locations, confirming that the GONGGA system improves the estimates of spatial and temporal variations in carbon fluxes by assimilating OCO-2 XCO₂ data. This dataset will improve the broader understanding of global carbon cycle dynamics and their response to climate change. The dataset can be accessed at https://doi.org/10.5281/zenodo.8368846 (Jin et al., 2023a).

1 Introduction

Atmospheric carbon dioxide (CO₂) concentrations are rapidly rising, mainly because of increases in anthropogenic emissions. Land and oceans can absorb substantial amounts of CO₂ and thus mitigate global warming. During the past decade (2012–2021), approximately one-fourth of total CO₂ emissions were absorbed by the land and oceans, respectively (Friedlingstein et al., 2022). However, there are large uncertainties in estimates of the size, spatial distribution, and interannual variability in land and ocean fluxes (Piao et al., 2009b; Eldering et al., 2017; Hauck et al., 2020; Piao et al., 2020). Accurate estimates of these fluxes at global and regional scales are essential for improving overall knowledge regarding the current status of the carbon cycle and projecting long-term changes (Zscheischler et al., 2017).

There are many methods for the estimation of global and regional carbon budgets, including the inventory method, the eddy covariance method, the ecosystem process modeling method, and the atmospheric inversion method (Piao et al., 2022). The first three methods upscale the site-level ground observations using statistical or process-based models; they are usually regarded as bottom-up approaches. In contrast, atmospheric inversion infers carbon fluxes by combining information from atmospheric CO₂ concentrations, prior flux estimates, and atmospheric transport (Bousquet et al., 2000; Gurney et al., 2002), which is regarded as a top-down approach. Atmospheric inversion is appropriate for assessments of global and regional carbon fluxes because spatiotemporal variations in atmospheric CO₂ concentrations contain the signatures of sources and sinks at large spatial scales. However, inversion accuracy is limited by the number and distribution of atmospheric CO₂ observations, uncertainties regarding the atmospheric transport model and the CO₂ emission inventories (such as fossil fuel combustion emissions), and insufficient knowledge of prior flux uncertainties (Liu et al., 2021; Piao et al., 2022).

Currently, atmospheric inversions use either ground-based or space-based observations. Ground-based in situ and flask observations have higher precision, but they are unevenly distributed. Most ground-based observations are mainly concentrated in North America and Europe (Peters et al., 2007; Chevallier et al., 2010; Lauvaux et al., 2016). Inversions using in situ and flask observations can consistently constrain surface CO₂ fluxes at the global scale and at some regional scales (for well-sampled continents), but their uncertainty rapidly increases at the sub-continental scale or when considering continents with sparse observations (Peylin et al., 2013; Byrne et al., 2017; Crowell et al., 2019). For example, there are only eight sites in the Chinese mainland under the World Meteorological Organization (WMO)–Global Atmosphere Watch (GAW) program (Wang et al., 2020b), and Chinese land sinks constrained by in situ CO₂ observations can differ by up to an order of magnitude (Chen, 2021; Wang et al., 2022a, b). The space-based column-averaged CO₂ dry-air mole fraction (XCO₂) retrievals serve as an emerging data stream for atmospheric inversions. Satellite XCO₂ retrievals have broader spatial coverage than in situ and flask observations; accordingly, they fill observational gaps over areas with few stations. The two most widely used satellites dedicated to measure CO₂ are Greenhouse gases Observing SATellite (GOSAT) (Yokota et al., 2009) and Orbiting Carbon Observatory-2 (OCO-2) (Crisp et al., 2004). GOSAT retrievals have been used in multiple inversions and were shown to be able to reduce the uncertainty in flux estimates in regions where surface CO₂ observations are sparse (Takagi et al., 2011; Basu et al., 2013; Chevallier et al., 2014). The OCO-2 team updates satellite retrievals roughly once per year. Refinements of instrument error characterization, retrieval algorithms, and bias correction procedures have led to substantial improvements in the accuracy and precision of satellite-retrieved XCO₂ data through these updates (O'Dell et al., 2018; Kiel et al., 2019); the single sounding random error in official OCO-2 retrievals is now better than 1 ppm (Eldering et al., 2017; Wunch et al., 2017). These improvements in XCO₂ retrievals have a transformative effect on satellite-based estimates of global carbon fluxes (O'Dell et al., 2018; Miller and Michalak, 2020). For example, the OCO-2 version 7 retrievals – the basis of early inversion studies using OCO-2 data – are fit to constrain land carbon fluxes at continental and hemispheric scales (Miller et al., 2018; Crowell et al., 2019). Chevallier et al. (2019) showed that the OCO-2 version 9 retrievals have a similar performance in terms of constraining carbon fluxes to the inversions that use observations from surface stations when the inversed fluxes and CO₂ concentrations are compared with independent aircraft data. More recently, the OCO-2 team has released the retrieval product for version 11r. The effectiveness and potential applications of these updated satellite retrievals in efforts to infer surface CO₂ fluxes require continuous and persistent investigation.

In this study, we used the GONGGA (Global ObservatioN-based system for monitoring Greenhouse GAses) inversion system (Jin et al., 2023b) to generate a global dataset of terrestrial ecosystem and ocean carbon fluxes from 2015 to 2022 by assimilating OCO-2 XCO₂ retrievals (v11r). Here, we present the prior and posterior global 3-hourly gridded terrestrial ecosystem and ocean carbon fluxes at a spatial resolution of 2° latitude × 2.5° longitude. Gridded fluxes from fossil fuel emissions and biomass burning emissions are also available for inferring the total fluxes.

This paper is organized as follows: Sect. 2 describes the methods and data used; Sect. 3 describes the format and content of the dataset; Sect. 4 analyzes the key characteristics of global and regional carbon cycles; Sect. 5 evaluates posterior fluxes using TCCON and ObsPack observations; Sect. 6 discusses the factors affecting the inversion results; Sect. 7 introduces data availability; and Sect. 8 summarizes the paper.

2 Methods and data

2.1 The GONGGA inversion system

GONGGA is an atmospheric inversion system that constrains gridded carbon fluxes with atmospheric CO₂ observations and transport simulations (Jin et al., 2023b). The assimilated observations are OCO-2 v11r XCO₂ retrievals (OCO-2/OCO-3 Science Team et al., 2022), and the transport model is GEOS-Chem v12.9.3 (Suntharalingam et al., 2004; Nassar et al., 2010, 2013). The spatial resolution of GEOS-Chem is 2° latitude × 2.5° longitude, with 47 layers in the vertical direction from the surface to the top of the atmosphere. The model is driven by Modern-Era Retrospective analysis for Research and Applications 2 (MERRA-2) meteorological data provided by the Goddard Earth Observing System (GEOS) of the National Aeronautics and Space Administration (NASA) Global Modeling and Assimilation Office (Gelaro et al., 2017). Four types of carbon fluxes are used to drive the atmospheric CO₂ simulations, including NEE (net ecosystem exchange, i.e., the balance of photosynthesis and respiration) from terrestrial ecosystems, atmosphere–ocean carbon exchange, fossil fuel carbon emissions, and biomass burning carbon emissions. NEE and ocean carbon fluxes are optimized by GONGGA, whereas fossil fuel emissions and biomass burning emissions are assumed to be well-known and not optimized, which is a usual convention in global atmospheric inversions (Peters et al., 2007; Jiang et al., 2022).

GONGGA uses the nonlinear least-squares four-dimensional variational data assimilation (NLS-4DVar) method (Tian and Feng, 2015; Tian et al., 2018) to minimize the following cost function:

$$\begin{array}{}\text{(1)}& \begin{array}{rl}J\left(\mathit{x}\right)& ={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\left(\mathit{x}-{\mathit{x}}_{\mathrm{a}}\right)}^{\mathrm{T}}{\mathbf{B}}^{-\mathrm{1}}\left(\mathit{x}-{\mathit{x}}_{\mathrm{a}}\right)\\ & +{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\left(\mathit{y}-H\left(\mathit{x}\right)\right)}^{\mathrm{T}}{\mathbf{R}}^{-\mathrm{1}}\left(\mathit{y}-H\left(\mathit{x}\right)\right),\end{array}\end{array}$$

where x is the state vector that contains the variables to be optimized and x_a is its prior estimate; y gathers the XCO₂ retrievals; B is the prior error covariance matrix and R is the observation error covariance matrix; and H(⋅) is the observation operator, which relies on GEOS-Chem simulations and sampling of modeled atmospheric CO₂. Firstly, the atmospheric transport model is used to simulate gridded CO₂ concentrations driven by surface fluxes. Then, the simulated gridded CO₂ profiles are interpolated horizontally by inverse distance weighting and vertically by linear interpolation. Thirdly, the interpolated CO₂ profiles are used to construct the simulated XCO₂ using the following equation:

$$\begin{array}{}\text{(2)}& {\mathrm{XCO}}_{\mathrm{2}}^{\mathrm{m}}={\mathrm{XCO}}_{\mathrm{2}}^{\mathrm{a}}+{\mathit{h}}^{\mathrm{T}}\mathbf{A}\left({\mathit{x}}_{\mathrm{CO}\mathrm{2}}-{\mathit{x}}_{\mathrm{CO}\mathrm{2},\mathrm{a}}\right),\end{array}$$

where ${\mathrm{XCO}}_{\mathrm{2}}^{\mathrm{m}}$ is the modeled XCO₂ and x_CO2 is the interpolated CO₂ profile from the GEOS-Chem simulation. ${\mathrm{XCO}}_{\mathrm{2}}^{\mathrm{a}}$, h, A, and x_CO2,a are the prior value of XCO₂, the pressure weighting function, the averaging kernel matrix, and the prior CO₂ vertical profile, respectively, provided by the OCO-2 Lite file.

The optimization algorithm NLS-4DVar is a hybrid method that combines the advantages of the conventional four-dimensional variational (4DVar) method and ensemble Kalman filter (EnKF), achieving high inversion accuracy with low computational cost and complexity (Tian and Feng, 2015; Tian et al., 2018). GONGGA adopts a novel dual-pass inversion strategy, successively optimizing initial CO₂ concentrations and surface carbon fluxes within each inversion window of 14 d, which distinguishes model–data mismatches caused by errors from initial CO₂ concentrations and surface fluxes (Jin et al., 2023b). Note that during the flux optimization, the state vector x gathers gridded scaling factors for NEE and ocean carbon fluxes. The spatial resolutions of the optimization for both initial CO₂ concentrations and fluxes are 2° latitude × 2.5° longitude, the same as the transport model resolution. The temporal resolution of the optimization is 14 d, indicating that the fluxes within each 14 d window are uniformly adjusted by the same scaling factor. In this study, GONGGA was run from 6 September 2014 to 31 December 2022. The 2014 results were regarded as spin-up, whereas the 8-year results spanning 2015–2022 comprised the dataset.

2.2 Prior CO₂ fluxes and uncertainties

The prior CO₂ fluxes include NEE, ocean-atmosphere carbon fluxes, fossil fuel emissions, and biomass burning emissions. The prior NEE was simulated by ORCHIDEE-MICT (Guimberteau et al., 2018). The prior ocean carbon fluxes were from the CT2022 pCO₂-Clim prior data, which were derived from the Takahashi et al. (2009) climatology of seawater pCO₂. Fossil fuel emissions were from the monthly Global Carbon Budget Gridded Fossil Emissions Dataset (GCP-GridFED; version 2023.1) (Jones et al., 2021). Biomass burning emissions were from the Global Fire Emissions Database (GFED; version 4.1s) 0.25° × 0.25° monthly data scaled with daily factors (Randerson et al., 2017; van der Werf et al., 2017). For the estimation of prior flux uncertainties, we first used a prior perturbation ensemble to approximate the prior scaling factor error covariance matrix and then calculated the prior flux uncertainties through the matrix multiplication between the scaling factor error covariance matrix and prior fluxes. The posterior flux uncertainties were calculated in the same manner using the ensemble of posterior scaling factors and prior fluxes. The difference between the prior and posterior flux uncertainties was regarded as the difference in the perturbation ensemble. For detailed steps, see Text S1 in the Supplement.

2.3 OCO-2 column CO₂ observations

We used OCO-2 Level 2 Lite v11r XCO₂ products (O'Dell et al., 2012, 2018; Gunson and Eldering, 2020) retrieved by the Atmospheric Carbon Observations from Space (ACOS) algorithm (Connor et al., 2008) to constrain the surface carbon fluxes. The OCO-2 satellite carries high-resolution spectrometers that return high-precision measurements of reflected sunlight received within the CO₂ and O₂ bands in the short-wave infrared spectrum (Crisp et al., 2012). The OCO-2 spacecraft flies in a 705 km altitude sun-synchronous orbit with a 16 d (233 orbits) ground track repeat cycle. OCO-2 has a footprint of 1.29 × 2.25 km² in nadir mode and acquires eight cross-track footprints, creating a swath width of 10.3 km.

Before assimilation, the XCO₂ retrievals were filtered with the xco2_quality_flag parameter provided by the OCO-2 Lite products; xco2_quality_flag = 0 (1) denotes good (bad) retrieval quality. Only retrievals with good quality were selected. We applied a data-thinning algorithm (Liu and Rabier, 2002; Campbell et al., 2017; Reale et al., 2018) to reduce the potential impacts of correlated errors in adjacent soundings. We set the threshold of the number of daily observations to 20 000. If the number of good retrievals exceeded the threshold within a single day, excess data were removed. For example, if there were 60 000 good retrievals in 1 d, one of every three sequential retrievals was selected according to sounding ID. Before data thinning, there were 203 368 424 XCO₂ retrievals with good quality from 6 September 2014 to 31 December 2022. After data thinning, 40 337 763 XCO₂ retrievals, which is about a fifth of total good retrievals, were actually assimilated in the inversion. Furthermore, OCO-2 retrievals were scaled to the official World Meteorological Organization (WMO) X2019 standards following instructions provided by the National Oceanic and Atmospheric Administration (NOAA; https://docs.google.com/document/d/e/2PACX-1vQ0JqK72fAOThaJwJyILLgfOE2qpHYdgNsIYAs6T2cM GumwVliSK7lurIYKCMOFgz1fyxuKYwlm5FEx/pub, last access: 12 September 2023).

2.4 Evaluation of posterior fluxes

Generally, it is difficult to directly verify the posterior fluxes because of the lack of direct flux observations that exhibit a footprint size comparable with the spatial resolution of global inversion models (typically several hundred kilometers). Instead, we compared the simulated CO₂ concentrations driven by posterior fluxes with atmospheric CO₂ observations to achieve indirect verification (e.g., Wang et al., 2019; Wu et al., 2020; Liu et al., 2021). In this study, we performed these comparisons using observations from TCCON version GGG2020 (Laughner et al., 2024) and ObsPack (CO₂ GLOBALVIEWplus v8.0 and NRT v8.1) datasets (Cox et al., 2022; Di Sarra et al., 2023).

2.4.1 TCCON XCO₂ retrievals

TCCON is a network of ground-based Fourier transform spectrometers that record direct solar spectra in the near-infrared spectral region. From these spectra, accurate and precise column-averaged CO₂ abundances are retrieved and reported (Wunch et al., 2011). TCCON XCO₂ retrievals are estimated to have precisions better than 0.25 % (i.e., ∼ 1 ppm) (Wunch et al., 2011). These retrievals have been used as primary validation data for several satellite missions, including GOSAT and OCO-2 (Wunch et al., 2011, 2017). Here, we used GGG2020 version data (Wunch et al., 2015). There are 27 TCCON sites with observations covering the inversion period (Table 1); the site locations are shown in Fig. 1a.

Table 1 Geographic locations and references of TCCON sites used for validation. Sites are listed according to latitude from north to south.

2.4.2 ObsPack CO₂ observations

ObsPack is a framework that combines atmospheric greenhouse gas observations from various sampling platforms (e.g., surface, aircraft, towers, or ships) and strategies (e.g., flask or in situ), ensuring consistent data quality (Masarie et al., 2014). In this study, we used the surface flask and aircraft observations from obspack_co2_1_GLOBALVIEWplus_v8.0_2022-08-27 in 2015–2021 and obspack_co2_1_NRT_v8.1_2023-02-08 in 2022, both of which are established according to the WMO X2019 scale (Cox et al., 2022). Surface flask observations are usually made on a weekly basis. During the 2015–2022 period, surface flask observations from 57 sites with parameter CT_assim = 1 or 2 were used for evaluation (Fig. 1b). Observations may be reported by multiple institutes at a single site. Here, we only used data from the NOAA laboratory and ignored duplicate records from other sites. Based on the spatial distribution of surface flask sites, we evaluated terrestrial carbon fluxes in six regions: North America, South America, Europe, Africa, East Asia, and Australia (Fig. 1b).

Aircraft observations contain data from the Comprehensive Observation Network for TRace gases by AIrLiner (CONTRAIL) program, the Atmospheric Tomography Mission (ATom), and several localized measurements concentrated in North America (Fig. 1b). The CONTRAIL program is Japan's unique aircraft observation project that measures atmospheric CO₂ concentrations using instruments on board Japan Airlines (JAL) commercial airliners. In the CO₂ GLOBALVIEWplus v8.0 ObsPack dataset, the CONTRAIL program contains aircraft measurements between Japan and Australia from 2015 to 2021. The ATom was a NASA Earth Venture Suborbital-2 mission that investigated the greenhouse gases in the atmosphere over the Pacific and Atlantic oceans from August 2016 to May 2018. We used CO₂ observations above the planetary boundary layer (altitude > 1 km) for evaluations to avoid effects on local emissions related to aircraft ascent and descent.

Figure 1 Spatial distributions of (a) TCCON sites and (b) ObsPack sites used for flux evaluations. Rectangles show the ranges of the six regions used for comparisons with surface flask observations.

3 Dataset description

Here, we present a global dataset that contains surface carbon fluxes from 2015 to 2022. The flux files contain NEE, ocean–atmosphere carbon fluxes, fossil fuel emissions, and biomass burning emissions. The NEE and ocean–atmosphere carbon fluxes include prior and posterior estimates. The corresponding gridded uncertainties in NEE and ocean–atmosphere fluxes are also included in the flux files. The global gridded fluxes are 3-hourly with a resolution of 2° latitude × 2.5° longitude.

4 Results

4.1 Global carbon budget

Here, we present the five major components of the global carbon budget, including the fossil fuel CO₂ emissions (E_FOS), biomass burning emissions (E_FIRE), atmospheric CO₂ concentration growth rate (G_ATM), ocean–atmosphere carbon fluxes (F_OCEAN), and NEE (Fig. 2). During the 2015–2022 period, E_FOS was 9.71 ± 0.20 Pg C yr⁻¹, with a minimum of 9.44 Pg C yr⁻¹ in 2020 and a maximum of 9.94 Pg C yr⁻¹ in 2022; E_FIRE was 1.86 ± 0.22 Pg C yr⁻¹, with a minimum of 1.47 Pg C yr⁻¹ in 2022 and a maximum of 2.14 Pg C yr⁻¹ in 2019. Over these 8 years, NEE exhibited a substantial mean sink with considerable interannual variability, estimated to be the standard deviation across years (−4.08 ± 0.53 Pg C yr⁻¹). The sinks were significantly weaker in 2015 and 2016 than in other years. The annual mean NEE in 2015–2016 was −3.35 Pg C yr⁻¹, whereas the NEE in 2017–2022 was −4.33 Pg C yr⁻¹. The reduced NEE during 2015–2016 was mainly related to the El Niño event, which caused substantial carbon release in the tropics (Wang et al., 2013; Liu et al., 2017; Piao et al., 2020; Dannenberg et al., 2021). Compared with NEE, interannual variation in the atmosphere–ocean fluxes was much smaller (−2.32 ± 0.18 Pg C yr⁻¹). From 2015 to 2022, the terrestrial biosphere (NEE+E_FIRE) and ocean absorbed approximately 23 % and 24 % of total fossil fuel CO₂ emissions, respectively, resulting in a G_ATM of 5.17 ± 0.68 Pg C yr⁻¹.

We compared the GONGGA-estimated global carbon budget with results from the measurements and Global Carbon Budget 2023 (hereafter referred to as GCB2023) (Friedlingstein et al., 2023). The G_ATM directly estimated from atmospheric CO₂ concentration measurements provided by the NOAA Earth System Research Laboratories' Global Monitoring Laboratory (ESRL/GML) was 5.25 ± 0.61 Pg C yr⁻¹ during 2015–2023 (Lan et al., 2023), which corroborates the GONGGA estimate. We also compared net biosphere exchange (NBE, i.e., the net carbon flux of all the land–atmosphere exchange processes except fossil fuel emissions, calculated as NEE+E_FIRE) and ocean sink estimated from the GONGGA inversion with GCB2023. Note that the GCB2023 estimations represent the carbon accumulated in the land and ocean reservoirs. We followed GCB2023's definitions and adjusted riverine CO₂ transport from the net atmosphere–surface CO₂ exchange over land (NBE) and ocean (F_OCEAN). Specifically, pre-industrial lateral carbon transport through the land–ocean aquatic continuum (LOAC) of 0.65 ± 0.35 Pg C yr⁻¹ (Regnier et al., 2022) was subtracted from −NBE to represent land carbon sink and added to −F_OCEAN to represent ocean carbon sink. During 2015–2022, the mean of corrected land carbon sink from GONGGA was 1.57 ± 0.67 Pg C yr⁻¹ and the mean of corrected ocean sink was 2.97 ± 0.18 Pg C yr⁻¹. GCB2023's estimate of ocean sink was 2.88 ± 0.07 Pg C yr⁻¹ based on global ocean biogeochemistry models and surface ocean fCO₂-observation-based products. The land carbon sink from GCB2023 was 2.00 ± 0.62 Pg C yr⁻¹ from the dynamic global vegetation models (DGVMs) and was 1.55 ± 0.77 Pg C yr⁻¹, calculated as the residual sink from the global budget of fossil fuel emissions, atmospheric growth rate, and ocean sink (Friedlingstein et al., 2023). As the estimate of land carbon sink from DGVMs introduces a budget imbalance in GCB2023, our estimates are well consistent with GCB2023's estimates based on ocean models and the residual land sink and close the global budget.

Figure 2 Global carbon budget estimated by GONGGA and atmospheric CO₂ growth rate from NOAA during 2015–2022.

4.2 Global distribution and regional fluxes

Figure 3 shows the global distributions of GONGGA annual mean NBE and ocean carbon fluxes during 2015–2022. Terrestrial carbon sinks were mainly in temperate North America, central South America, southern Africa, Europe, boreal Asia, India, eastern China, and most of Australia. Terrestrial carbon sources mainly occurred over western America, the eastern Amazon, central Africa, Southeast Asia, the southeastern coast of Australia, and New Zealand. The ocean sources mainly occurred over tropical oceans and the high-latitude Southern Ocean; the equatorial Pacific was the most prominent source area. Sinks mainly occurred over mid-latitude regions of both hemispheres and the high-latitude Northern Ocean. Generally, NBE had a more complex spatial distribution and higher uncertainty, compared with ocean carbon fluxes. Therefore, we explored the distribution and attribution of NBE over 11 TransCom land regions (Fig. 4) (Gurney et al., 2004).

Figure 3 GONGGA-estimated global distributions of annual mean (2015–2022) NBE and ocean carbon fluxes.

Figure 4 Spatial distributions of 11 TransCom land regions.

Here, we present the GONGGA-estimated annual mean (2015–2022) NBE for 11 TransCom land regions and their comparison with OCO-2 model intercomparison project (MIP) v10 inversions (Fig. 5). OCO-2 MIP v10 (Baker et al., 2023; Byrne et al., 2023) includes an ensemble of 14 atmospheric inversions over the period of 2015–2020, assimilating OCO-2 v10r retrievals, and each of them is characterized by distinct transport models, data assimilation algorithms, and prior fluxes (Table S1 in the Supplement). We used OCO-2 MIP v10 results from the inversions that assimilate land nadir and land glint (LNLG) satellite retrievals, and those assimilate in situ (IS) measurements. Here, in situ inversions are used to provide a baseline against satellite-driven results.

For the 11 TransCom regions, we estimated that Europe had the strongest terrestrial carbon sink, followed by Boreal Asia, Temperate Asia, Temperate North America, Temperate South America, Southern Africa, Boreal North America, and Australia, whereas Tropical South America, Northern Africa, and Tropical Asia were terrestrial carbon sources. All GONGGA and OCO-2 MIP LNLG and IS inversions consistently indicated that Europe was the largest terrestrial sink. GONGGA showed good agreement with OCO-2 MIP inversions for most regions, and divergences occurred mainly in Boreal North America and Northern Africa. The difference between GONGGA and OCO-2 MIP inversions may be related to the prior NBE adopted and retrieval pre-processing methods utilized. In Boreal North America, GONGGA's prior emerged as a carbon source, whereas OCO-2 MIP's prior was a carbon sink (Fig. S1). After assimilating OCO-2 retrievals, GONGGA and OCO-2 MIP consistently showed Boreal North America was a carbon sink, but the sink in GONGGA was smaller than OCO-2 MIP. The same situation happened in Northern Africa. Both GONGGA's prior and OCO-2 MIP's prior estimated Northern Africa as a terrestrial carbon sink, but the sink from GONGGA was stronger than that from OCO-2 MIP (Fig. S1 in the Supplement). Constrained by OCO-2 retrievals, both GONGGA and OCO-2 MIP estimated it as a carbon source, and the source from GONGGA was weaker than that from OCO-2 MIP, aligning with the sizes of their prior sinks. In addition, the impact of prior fluxes may be amplified by the insufficient coverage of OCO-2 retrievals. For example, in Boreal North America, satellites cannot measure XCO₂ in dark high-latitude areas in winter. In Northern Africa, OCO-2 also has difficulties in accurately measuring XCO₂ over the desert because of its high albedo, demonstrated by its high proportion of bad retrievals (xco2_quality_flag = 1) (Zhang et al., 2016). Therefore, the posterior fluxes in these two regions were more dependent on the prior fluxes during the period with few OCO-2 retrievals. Notably, even in OCO-2 MIP inversions, the ensemble spread was prominent, indicating the difficulty of inversion in these regions using current satellite or in situ observations (Table S2).

Figure 5 GONGGA annual mean (2015–2022) NBE for 11 TransCom land regions and comparisons with OCO-2 MIP LNLG and IS inversions. Error bars represent standard deviations in the annual mean budget across the whole period.

4.3 Interannual variability and seasonal cycle

Here, we analyzed the interannual variability (IAV) and seasonal cycle of NBE at global and regional scales. We divided the globe into three large latitude bands: northern extratropics (30–90° N), tropics (30° S–30° N), and southern extratropics (90–30° S). The global net terrestrial carbon flux has a prominent year-to-year variability (Friedlingstein et al., 2022). We computed the standard deviation of global NBE to represent its magnitude of IAV, which amounted to 0.63 Pg C yr⁻¹ during the 2015–2022 period. The variations in NBE in northern extratropics, tropics, and southern extratropics were quite different (Fig. 6). We calculated the contribution of each latitude band to the global IAV using Eq. (1) from Ahlström et al. (2015). The contribution of the tropics to the global NBE IAV was 100.8 %, whereas the contributions of the northern and southern extratropics were −13.2 % and 12.4 %, respectively. A positive (negative) score here indicates the variation is in the same phase as (opposite phase to) the global IAV. The scores from our estimate indicated that the global IAV arises from the tropics. Given the short time series of the inversion, the latitudinal contributions in this study are suggestive but not statistically conclusive. The dominant role of tropical terrestrial ecosystems in the signal of the global carbon cycle IAV is consistent with previous results based on multiple observations and models (Baker et al., 2006; Rödenbeck et al., 2018b, a; Jung et al., 2020). Piao et al. (2020) reviewed and analyzed the regional contribution to global net terrestrial carbon flux IAV from 1980 to 2017 with process-based land carbon cycle models, atmospheric inversion models, and FLUXCOM data products. The contributions of the tropics to the global IAV obtained by these three methods were 83.4 %, 71.7 %, and 69.7 %, respectively. In addition to the short time series, the inclusion of the 2015–2016 strong El Niño event in the period is an important reason for the large contribution score of the tropics in our estimate. Climatic variations are the main factors that drive the IAV of the net terrestrial carbon flux (Braswell et al., 1997; Zeng et al., 2005; Raupach et al., 2008; Liu et al., 2017). El Niño is the major climatic mode that modulates global temperature, precipitation, and solar radiation (Gu and Adler, 2011); thus, it drives the IAV of the carbon cycle (Bacastow, 1976; Rayner et al., 2008). The characteristics of hot and dry climate conditions in El Niño years are the primary reasons for the lower net carbon uptake or net carbon release by terrestrial ecosystems (Jones et al., 2001; Piao et al., 2009a), which is particularly evident in the tropics (Fig. S2) (Liu et al., 2017; Jin et al., 2023b). During 2015–2016, tropical land released 0.66 Pg C yr⁻¹ CO₂ into the atmosphere, whereas it is a net terrestrial sink in normal years (−0.52 Pg C yr⁻¹).

Figure 6 Annual NBE anomaly over the globe, northern extratropics (30–90° N), tropics (30° S–30° N), and southern extratropics (90–30° S) during 2015–2022. The shaded area represents the uncertainty in NBE in each region.

The shape of the NBE seasonal cycle varies among regions and different years. In the northern extratropics, the size and phase of the seasonal cycle are very similar in all years, with July having the largest sink and northern winter being a carbon source. In the tropics, however, the seasonal cycles have smaller amplitudes and the shapes are distinct in different years. The largest deviations of the tropical seasonal cycle from the 8-year mean estimate are in 2016 (coefficient of determination between annual mean seasonal cycle and the year investigated R²=0.34) and 2019 (R²=0.50); the most prominent deviations occurred during the peak 2015–2016 El Niño between July 2015 and June 2016 as well as 2019 El Niño between April 2019 and July 2019. The shape of the global seasonal cycle is nearly similar to the shape of the northern extratropics (with 103.6 % contribution), whereas the tropics and southern extratropics have opposite phases compared with the global seasonal cycle (with −1.1 % and −2.5 % contribution, respectively). The dominance of the northern extratropics in the global seasonal cycle is consistent with previous findings (Forkel et al., 2016; Piao et al., 2020).

The amplitude is an important index of the seasonal cycle (i.e., seasonal cycle amplitude, SCA). The peak-to-trough amplitude was calculated as the difference between the maximum and minimum monthly NBE in each year. The 8-year mean SCAs of NBE for the globe, northern extratropics, tropics, and southern extratropics were 3.55, 3.50, 0.43, and 0.12 Pg C month⁻¹, respectively. The larger mean amplitude in northern land ecosystems, compared with other regions, was mainly related to the strong seasonality of gross primary production and ecosystem respiration (Randerson et al., 1997).

Figure 7 Seasonal cycle of NBE for the globe, northern extratropics, tropics, and southern extratropics during 2015–2022.

5 Dataset evaluation

5.1 Comparison with TCCON observations

In this section, we compare the simulated monthly mean XCO₂ driven by the posterior CO₂ fluxes with the retrievals from 27 TCCON sites during 2015–2022 (Table 1). The global mean root mean square error (RMSE) and bias between posterior simulated and TCCON retrieved XCO₂ were 0.81 and 0.24 ppm, respectively. Through the assimilation of OCO-2 retrievals, the atmospheric CO₂ simulations were considerably improved compared with prior simulations, which exhibited 1.15 ppm RMSE and 0.51 ppm bias at the global scale. At most sites, posterior RMSE was < 1 ppm and bias was in the range of −0.5 to 1 ppm (Fig. 8). The maximum simulation deviation occurred at Eureka station (unless otherwise stated, “simulations” hereafter refers to posterior simulations, which means the simulation is driven by posterior fluxes), where an overestimation of simulated XCO₂ was observed in winter. This overestimation was also evident at Ny Ålesund and Sodankylä, which are located in the high latitudes of the Northern Hemisphere. Prior simulations generally overestimated CO₂ concentrations, particularly in winter (Fig. S3). Positive deviations were adequately mitigated at most sites after the inversion. However, for the sites mentioned above, considering the lack of satellite retrievals in winter at high northern latitudes, the posterior flux may be poorly constrained and is thus similar to the prior flux. Additionally, the coarse spatial resolution of the transport model is another challenge for the detection of sub-grid variations in XCO₂. For example, Edwards station and Pasadena station are close to each other, and they are located in the same grid cell of the transport model. The simulated XCO₂ time series at these two sites are similar, and the minor difference mainly arises from the interpolation process (Fig. S4). In contrast, the TCCON XCO₂ retrievals are considerably higher at Pasadena station than at Edwards station, with a multi-year mean difference of 0.84 ppm.

Figure 8 Spatial distributions of (a) root mean square error (RMSE) and (b) bias between the posterior monthly XCO₂ simulations and corresponding observations at each TCCON site (simulations minus observations; unit: ppm).

Figure 9 Time series of monthly averaged observations and posterior simulations at each TCCON site.

5.2 Comparison with ObsPack observations

Here, we compare posterior CO₂ simulations with ObsPack surface flask and aircraft observations. The global mean RMSE and bias between surface flask observations and corresponding simulations were 1.76 and −0.33 ppm, respectively. For most surface flask sites located on the ocean and in tropical and southern extratropical terrestrial regions, RMSE was < 2.0 ppm; bias was in the range of −0.5 ppm to 0.5 ppm. The high RMSE values mainly occurred over northern middle latitudes, particularly over Europe and East Asia. Jiang et al. (2022) used GOSAT XCO₂ retrievals to estimate global CO₂ fluxes and also found that posterior CO₂ concentrations could differ from surface observations, mainly in the northern extratropics. Because of limitations regarding the coarse resolutions of global transport models and thus differences in representativeness between simulated CO₂ concentrations and actual observations over land, some sites have significant data–model mismatches. For example, at the three sites with posterior RMSE values exceeding 4.0 ppm, the observed atmospheric CO₂ concentrations had strong temporal fluctuations, which were presumably caused by localized and short-term surface fluxes (Fig. S5).

Figure 10 Spatial distributions of (a) RMSE and (b) bias between the posterior monthly XCO₂ simulations and corresponding observations at each surface flask site (posterior simulations minus observations; unit: ppm).

To decrease mismatches in temporal and spatial representativeness between observations and simulations, we compared the monthly observed and simulated CO₂ concentrations in six land regions (Fig. 11). Apart from RMSE and bias, we further present the random error here, which was calculated as the standard deviation of the differences between simulated and observed CO₂ concentrations (Rastogi et al., 2021). The monthly simulations closely agreed with the observations. The bias was in the range of −0.44 to −1.27 ppm, random error was in the range of 0.39 to 1.65 ppm, and RMSE was in the range of 0.58 to 2.08 ppm. The simulation deviations remained higher for North America, Europe, and East Asia, compared with other regions. In these three regions, there was a significant difference in terms of comparisons with TCCON and ObsPack surface flask observations. Mainly positive bias arose from TCCON evaluations and negative bias arose from ObsPack evaluations. This discrepancy may be related to the nature of the two types of observations. TCCON observations are column-averaged atmospheric CO₂ concentrations, whereas ObsPack observations are surface atmospheric CO₂ concentrations. The opposite signs of bias between the two comparisons may be related to the imperfect simulation of the vertical mixing of GEOS-Chem (Schuh et al., 2019).

Figure 11 Time series of monthly averaged ObsPack surface flask observations and corresponding posterior simulations for the six sub-regions.

For aircraft observations, we calculated the mean statistics of each grid cell (Fig. 12). The simulations closely agreed with the aircraft observations. For most grid cells, the RMSE was < 2.0 ppm; the bias was between −1.0 and 1.0 ppm. The simulated deviations over Boreal North America and Temperate North America were generally larger than over the ocean, similarly to the surface flask results. We also compared the vertical distribution of modeled CO₂ against the observations. Figure S6 shows that the random errors were typically within 1 ppm, showing a good agreement between the simulations and observations. However, large biases of up to 2 ppm were seen in the high latitudes and above 9 km, which is consistent with the comparisons against the TCCON retrievals (Fig. 8).

Figure 12 The (a) random error and (b) bias between posterior CO₂ simulations and aircraft observations at each grid cell (posterior simulations minus observations; unit: ppm).

6 Discussion

Regarding the regional carbon budget, we found that fire emission, although it was not optimized in the inversion, largely impacted the net CO₂ fluxes from terrestrial ecosystem, i.e., NBE, in equatorial regions and Australia. With the frequent occurrence of wildfires in these regions, carbon emissions from wildfires may exceed regional NEE and make these regions net carbon sources (Fig. S7). For the past few decades, ∼ 50 % of fire-related carbon emissions and ∼ 70 % of global burned areas occurred across African subtropical savannah systems (Giglio et al., 2013; Andela and van der Werf, 2014). In the Amazon, despite the decline in deforestation rate during 2003–2015, carbon emissions from drought-induced fires increased very quickly (Aragão et al., 2018). Southeast Australia also experienced intensive and geographically extensive wildfires during the 2019–2020 summer season, and the fires released substantial amounts of CO₂ into the atmosphere (Wang et al., 2020a; Byrne et al., 2021; van der Velde et al., 2021). As a result, the 8-year mean biomass burning emissions in Tropical South America, Northern Africa, and Tropical Asia amounted to 0.17, 0.33, and 0.13 Pg C yr⁻¹, and were 6.2, 1.2, and 1.4 times higher than regional NEE, respectively, resulting in net carbon sources in these regions. The increasing fire emissions thus present a great challenge to climate mitigation efforts.

The processing of XCO₂ uncertainties also had an impact on the inversion results. We performed three sensitivity inversions with different XCO₂ uncertainties. The XCO₂ uncertainties were inflated 2 and 4 times in the first (E1) and second (E2) test, respectively. In the third test (E3), the XCO₂ uncertainties were increased by 5 ppm. The three sensitivity tests adopted the same configuration as the reference inversion in this study only except for the XCO₂ uncertainties. The distributions of different XCO₂ uncertainties were shown in Fig. S8. At the global scale, the inverted annual NBE and F_OCEAN from the original inversion, E1, and E2 were very close, but E3 had a different partitioning between land and ocean fluxes than the other inversions, which amounted to about 0.2 Pg C yr⁻¹ (Fig. S9). When it comes to the regional scale, the differences were larger in some regions and years but were still broadly consistent with the reference inversion (Fig. S10). These sensitivity tests highlighted the fact that the inversion results were indeed impacted by the assumption regarding XCO₂ uncertainty and careful assessment of uncertainties in satellite XCO₂ retrievals was necessary for accurate estimates of global and regional carbon fluxes.

In the current version of GONGGA, we assimilated the OCO-2 v11r Lite XCO₂ dataset. A recent paper found that the v11r Lite product has a bias of −0.4 to −0.8 ppm across regions north of 60° N due to the variations in the digital elevation model (DEM) used in the retrieval algorithm (Jacobs et al., 2024), and this bias introduces a ∼ 100 Tg C shift in the partitioning of carbon fluxes for the latitudinal bands. A preliminary test of GONGGA using the latest v11.1r Lite product showed the inverted terrestrial carbon sink tends to be 90 to 140 Tg C yr⁻¹ lower north of 60° N than when using the v11r Lite product, which is consistent with the previous findings. In addition, some parts of GONGGA's inversion algorithm, such as the data selection, were partly different from those proposed by the OCO-2 Science Team (Baker et al., 2022; Peiro et al., 2022; Byrne et al., 2023), but GONGGA's inversion results were broadly consistent with the ensemble of OCO-2 MIP inversions and GCB2023 and gave reasonable estimates of global and regional carbon budgets within the uncertainties. In the future, GONGGA will regularly publish new versions of inverted fluxes using the latest OCO-2 data on an annual basis. These updates will align with the latest suggestions from the OCO-2 Science Team, enabling the ongoing monitoring of CO₂ fluxes.

7 Data availability

The dataset is available at https://doi.org/10.5281/zenodo.8368846 (Jin et al., 2023a). As the satellite XCO₂ retrievals, prior carbon fluxes, and meteorological data are being persistently improved and updated, we plan to update the dataset annually in the future, aiming to support scientific research and policy-making.

8 Summary

Here, we present a global resolved surface carbon flux dataset for the 2015–2022 period. The dataset includes 3-hourly gridded (2° latitude × 2.5° longitude) NEE and ocean carbon fluxes (prior and posterior) together with prescribed fossil fuel emissions and biomass burning emissions. The dataset was generated by the GONGGA inversion system constrained by OCO-2 XCO₂ retrievals. We analyzed the key characteristics of the global and regional carbon cycles in terms of carbon budget, interannual variability, and seasonal cycle. The global annual estimate from GONGGA was consistent with the estimate from GCB2022. Regional fluxes were analyzed based on TransCom partitions. The strongest carbon sinks were observed in Europe, followed by Boreal Asia and Temperate Asia. We validated posterior fluxes by comparing posterior simulated CO₂ concentrations with TCCON XCO₂ retrievals as well as ObsPack surface flask and aircraft observations. Both evaluations demonstrated the optimization of posterior fluxes through assimilation of OCO-2 satellite retrievals.