Five years of variability in the global carbon cycle: comparing an estimate from the Orbiting Carbon Observatory-2 and process-based models

We use 10 s averages of bias-corrected total column CO₂ (XCO₂) generated from version 9 OCO-2 satellite observations. The 10 s averaged XCO₂ retrievals have previously been used in the OCO-2 inverse model inter-comparison project (e.g. Basu et al 2018, Chevallier et al 2019, Crowell et al 2019). Specifically, we include both land nadir and land glint retrievals in the model for years 2015 through 2019.

We couple a GIM and a global adjoint model (version v35n of the GEOS-Chem adjoint, Henze et al 2007) to estimate five years of global CO₂ fluxes and associated uncertainties at a daily temporal resolution and a spatial resolution of 4° (latitude) by 5° (longitude). We define negative values of fluxes as net carbon uptake by the land, and positive values thus represent a net release from the land to the atmosphere. In this section, we provide an overview of the approach, and the supporting information text S1 (available online at stacks.iop.org/ERL/16/054041/mmedia) provides additional description of detailed model setup and specific equations.

In a GIM, we estimate surface fluxes that will best match atmospheric observations using an atmospheric transport model:

$$\begin{equation}z = h\left( s \right) + \varepsilon \end{equation} \tag{ 1 }$$

The fluxes s (dimensions m × 1), when run through an atmospheric model, $\,h\left( {{\boldsymbol{s}}} \right)$, should match the observations $z$ (dimensions n × 1) within a specific error $\epsilon$ (dimensions n × 1).

A unique aspect of the GIM is that we can incorporate an array of environmental drivers to help describe the fluxes (s) instead of prescribing a traditional prior flux model. The GIM will scale the environmental drivers to best reproduce the atmospheric observations, and this component of fluxes is referred to as 'the deterministic model'. Furthermore, the GIM also estimates flux patterns that cannot be explained by the environmental drivers but are still evident in the atmospheric observations, and this component of fluxes is referred to as the 'stochastic component':

$$\begin{equation}s = {\mathbf{X}}\beta + \zeta \end{equation} \tag{ 2 }$$

where X is a matrix of environmental drivers (dimensions m × p; see section 2.3), and the drift coefficients β (dimensions p × 1) scale these environental drivers. Collectively, ${\mathbf{X}}{{\boldsymbol{\beta }}}$ is the 'determinstic model'. The unknown vector ${{\boldsymbol{\zeta }}}$ (dimensions m × 1) is the stochastic component and is estimated at the model grid scale. The posterior flux estimate ($\widehat {{\boldsymbol{s}}}$, dimensions m × 1) is a sum of the deterministic model $\left( {{\mathbf{X}}\beta } \right)$ and the stochastic components (${{\boldsymbol{\zeta }}}$).

We simultaneously estimate the fluxes (${{\boldsymbol{s}}}$) and the coefficients ($\beta$) via minimizing the GIM cost function (e.g. Michalak et al 2004):

$$\begin{align}{L_{s,\beta }} &= \frac{1}{2}{\left( {z - h\left( s \right)} \right)^T}{{\mathbf{R}}^ - }^1\left( {z - h\left( s \right)} \right) \nonumber\\ &\quad+ \frac{1}{2}{\left( {s - {\mathbf{X}}\beta } \right)^T}{{\mathbf{Q}}^ - }^1\left( {s - {\mathbf{X}}{\beta} } \right).\end{align} \tag{ 3 }$$

The cost function includes two covariance matrices; R (dimensions n × n) and Q (dimensions m × m). The covariance matrix R describes ${{\boldsymbol{z}}} - h\left( {{\boldsymbol{s}}} \right)$, referred to here as the model-data mismatch errors. These errors include errors from the atmospheric measurements z and from the transport model $h()$. The covariance matrix Q prescribes the variances and spatiotemporal covariances of the stochastic component (${{\boldsymbol{\zeta }}}$) and includes both diagonal and off-diagonal elements. supporting information text S1.1 describes the estimation of the covariance matrix parameters in detail.

In total, we assimilate five years of OCO-2 observations (i.e. n = ∼3.7 × 10⁵) to estimate ∼6.0 × 10⁶ unknown fluxes (m) at the model grid scale. We do so by minimizing equation (3), a process described in the supporting information text S1. We also estimate the posterior uncertainties in the flux estimate ($\widehat {{\boldsymbol{s}}}$) using a reduced rank algorithm described in Saibaba and Kitanidis (2015) and Miller et al (2020), an approach described in detail in supporting information text S1.2.

We further compare the GIM flux estimates against aircraft-based observations of CO₂ and against CO₂ observations from the Total Carbon Column Observing Network (TCCON; Wunch et al 2011), as an external means to evaluate the inverse model estimates using OCO-2. supporting information text S7, figures S3–S10, and tables S3 and S4 display the details of these comparisons.

We group the globe into seven biomes (figure 1(b)); for each biome, we consider an array of environmental drivers to include in the GIM (i.e. in X) available from NASA's Modern-Era Retrospective Analysis for Research and Applications, Version 2 (MERRA-2, Rienecker et al 2011). Specifically, these environmental drivers include daily 2 m air temperature, daily precipitation, 30 d average precipitation, photosynthetically active radiation (PAR), surface downwelling shortwave radiation, soil temperature at 10 cm depth, soil moisture at 10 cm depth, specific humidity, and relative humidity. We also include a non-linear function of air temperature (refer to hereafter as scaled temperature) from the Vegetation Photosynthesis and Respiration Model (Mahadevan et al 2008); in brief, this function characterizes the non-linear relationships between temperature and photosynthesis (e.g. Raich et al 1991). Note that we estimate a different coefficient ($\beta$) for each environmental driver in each biome and each year of the study period. Miller and Michalak (2020) argued that current OCO-2 observations are best-equipped to constrain terrestrial CO₂ fluxes for seven global biomes, and we therefore use seven biomes in this study.

Zoom In Zoom Out Reset image size

We also incorporate additional fluxes in X that do not necessarily map onto any environmental drivers, including fossil fuel emissions from the Open-source Data Inventory for Anthropogenic CO₂ monthly fossil fuel emissions (ODIAC, Oda et al 2018), oceanic fluxes from the Estimating the Circulation and Climate of the Ocean consortium (ECCO-Darwin; Brix et al 2015, Carroll et al 2020), and biomass burning fluxes from Global Fire Emissions Database (GFED) version 4.1 (Giglio et al 2013). We find that the estimated coefficients (β) for ODIAC, ECCO-Darwin, and GFED are very uncertain if we attemp to constrain those coefficients (β) for each of these sources in each biome. Hence, we instead estimate a single coefficient for all three sources for the entire globe.

Furthermore, we include a constant column of ones in X for each biome and each year (including for the oceans). These columns help describe the mean behavior of the fluxes in each biome, and the environmental variables in X describe deviations from that mean behavior. All existing GIM studies to date include constant components within X, similar to the setup here (e.g. Gourdji et al 2008, 2012, Shiga et al 2018).

We utilize a model selection approach to objectively determine the optimal subset of environmental drivers that can best reproduce OCO-2 observations. We specifically employ a type of model selection known as the Bayesian Information Criterion (BIC), a commonly used statistical approach in regression analyses (e.g. Kass and Raftery 1995, Raftery 1995) and inverse modelling studies (e.g. Gourdji et al 2012, Miller et al 2013, Fang and Michalak 2015, Miller et al 2016, 2018). In brief, the BIC rewards combinations of environmental drivers that better reproduce OCO-2 observations and penalizes combinations with too many environmental drivers to prevent overfitting; the best combination of environmental drivers is the combination with the lowest score. Miller et al (2018), Miller and Michalak (2020), and supporting information text S2 describes the specific setup for the BIC in greater detail.

We compare the inverse modeling estimates of IAV during the 5 year study period against an ensemble of 16 TBMs participating in TRENDY (v8; Sitch et al 2015, Friedlingstein et al 2019). Specifically, we collect net biospheric production (NBP) from TRENDY S3 simulations, which are forced over years 1901–2018 with changing CO₂, climate, and land use; NBP from TRENDY represents net CO₂ fluxes that take into account photosynthesis, plant and soil respiration, fire, land use change, and any carbon fluxes that are in and out of the ecosystem ( https://sites.exeter.acuk/trendy ). We compare against 5 years of model outputs from TRENDY to match the number of years in the inverse model. At the time of writing, TRENDY models are available through year 2018, so we use model outputs for years 2014–2018 (see table S1 for a full list of TBMs; Le Quéré et al 2018, Friedlingstein et al 2019, Piao et al 2020). Note that the climate forcing data used in TRENDY models are from Climatic Research Unit (CRU) and Japanese Reanalysis (JRA) datasets (CRUJRA; Harris et al 2014, Kobayashi et al 2015, Harris et al 2020).

We find large differences in the magnitude of IAV (defined here as the standard deviation of annual flux totals) between the results from OCO-2 and the TBMs (figure 2). Most TBMs simulate larger IAV compared to the estimate using OCO-2 across the five-year study period. Furthermore, TBMs show little consensus on IAV across different biomes and the globe. With that said, the ensemble mean of these TBMs across different biomes and the globe are generally within the uncertainty bounds of IAV estimated using OCO-2, but numerous individual TBMs still fall outside the uncertainty bounds; this result implies that there is an opportunity to evaluate estimates of IAV within individual TBMs using current satellite observations of CO₂, at least during the time period for which satellite observations are available.

Zoom In Zoom Out Reset image size

The estimate using OCO-2 further indicates that tropical biomes dominate global IAV and account for 87.8% of the variability during the 5 year study period. By contrast, most TBMs estimate a lower relative contribution (%) to global IAV from the tropics during this time, especially from the tropical grasslands biome (figure 3). Using OCO-2, we also find large variations in the carbon cycle associated with the 2015–2016 El Niño and subsequent recovery; these perturbations dominate global IAV during 2015–2019 and account for the very large contribution of the tropics to IAV in the OCO-2 estimate (figures 1(a) and 2). Note that we estimate regional contributions (figure 3) to global IAV using a contribution framework developed by Ahlström et al (2015) (supporting information text S5); in brief, we weigh the flux anomaly (i.e. departure from a 5 year mean) from each individual region by its similarity to the global anomaly, and therefore enable a direct comparison of their relative importance to global IAV.

Zoom In Zoom Out Reset image size

Specifically, the estimate using OCO-2 suggests that the tropical grasslands and tropical forests biomes make comparable contributions to global IAV during 2015–2019 (figure 3). Note here the tropical grasslands biome in Olson et al (2001; figure 1(b)) broadly encompasses grasslands, savanna, and agricultural land ecosystems (see figure 1 in Teluguntla et al 2015 or figure 16.1 in Thenkabail et al 2011) within the tropics. This finding differs from the TBMs, a majority of which (nine out of 16) estimate the largest contribution to global IAV from tropical forests (figure 3); a smaller number of TBMs (five out of 16) estimate the largest contribution from tropical grasslands.

Conversely, results from OCO-2 indicates that the extra-tropics accounts for a small fraction of recent, global IAV (12.2%). This result also disagrees with the estimates from the TBMs; most TBMs (12 out of 16) estimate larger IAV across the extra-tropics compared to the estimate using OCO-2 (figure 3). In fact, using OCO-2 we find a negative contribution from temperate grasslands to global IAV, which indicates the regional anomaly from temperate grasslands is in the opposite direction from the global anomaly. A previous study (Liu et al 2018) explored drought impact on CO₂ fluxes over the contiguous U.S. (CONUS), and suggested that the regional flux anomaly (departure from a 6 year mean) over the CONUS drought-impacted regions are in opposite directions of the global atmospheric CO₂ growth rate anomaly. The overall contribution of the extra-tropics is a sum of both positive and negative contributions of smaller-scale regions (e.g. Ahlström et al 2015). The negative contribution from temperate grasslands, therefore, acts to reduce the overall contributions of the extra-tropics to global IAV during years 2015–2019.

There is a growing consensus that tropical regions are the largest contributor to global IAV (e.g. Bousquet et al 2000, Peylin et al 2005, Gurney et al 2008, Jung et al 2017), but there is an ongoing debate over whether temperature or precipitation is a stronger correlate with this IAV in tropical regions (e.g. Zhao et al 2010, Wang et al 2014, Ahlström et al 2015, Jung et al 2017, Humphrey et al 2018). The findings using OCO-2 suggest that both environmental drivers play a key role in describing CO₂ variability across multiple scales—both in describing daily, grid-scale variability in CO₂ fluxes and in describing IAV, at least during the years when OCO-2 observations are available.

At the daily, model grid scale, a combination of temperature and precipitation can describe a substantial portion of the variability in CO₂ fluxes over the tropical grassland and tropical forest biomes (i.e. ∼78%–87% and ∼71%–83% of the flux variance across the study years, respectively). We define grid-scale variability as any spatiotemporal patterns in CO₂ fluxes that manifest at the resolutions of the GEOS-Chem model (daily, 4° (latitude) × 5° (longitude)) during the 5 year study period. Using the statistical model selection, we only select scaled temperature and precipitation in tropical biomes, further indicating the explanatory power of these environmental drivers. Figures 4(c)–(e) illustrate the sets of environmental drivers that are selected across individual biomes and their contributions to IAV in the inverse modeling estimate for years 2015–2019 (supporting information text S2, S3 and table S5 describe the full results of model selection for all biomes.).

Zoom In Zoom Out Reset image size

An analysis of the 2015–2016 El Ni$\tilde n$o further illustrates the correlation between these environmental drivers and IAV in the posterior flux estimate across tropical regions. The 2015–2016 El Ni$\tilde n$o induced anomalously dry and warm environmental conditions in the tropics (e.g. Jiménez-Muñoz et al 2016) that altered global biospheric CO₂ fluxes (e.g. Liu et al 2017). We find that differences in temperature and precipitation in the deterministic component of the flux estimate (X β, figures 4(c)–(e)) between years 2015 and 2017 contribute most to the grid-scale differences in CO₂ fluxes between the 2 years (figure 4(b)). By contrast, the stochastic component ${{\boldsymbol{\zeta }}}$ (figure 4(f)) describes relatively little year-to-year differences in the estimated fluxes. This result indicates that most year-to-year differences in the fluxes, as estimated using OCO-2, can be described by temperature and precipitation across tropical regions.

We also perform additional analysis to explore the impact of possible legacy effects, and disturbance (i.e. fire and land cover change) on IAV during the 5 year study period. We have limited ability in detecting legacy effects on IAV using five years of OCO-2 observations. Furthermore, fire and land cover change appear to play a small role in IAV for a 5 year period of this study. We discuss these topics in detail in supporting information text S3 and S4.

To better understand the differences of estimated IAV between OCO-2 and the TBMs, we further examine the relationships between IAV and annual anomalies in key environmental drivers during the 5 year study period. We explicitly examine the relationships with annual anomalies in temperature and precipitation across the tropics, a region that makes the largest contributions to global IAV during this time period (figure 5).

Zoom In Zoom Out Reset image size

We find that CO₂ fluxes in most TBMs tend to have stronger relationship with annual anomalies in temperature and precipitation than estimated using OCO-2 (figure 5). This analysis likely helps to explain the differences of IAV between the TBMs and OCO-2. Indeed, TBMs that have stronger relationships with anomalies in temperature and precipitation relative to the estimate using OCO-2 tend to also display a larger magnitude of IAV over tropical biomes (figure 2).

Furthermore, some TBMs show different relationships with environmental drivers but still arrive at the same IAV. For example, VISIT shows strong relationship with temperature across tropical grasslands while LPX-Bern shows strong relationship with precipitation; however, estimated IAV over tropical grasslands between the two models are in close agreement. This finding further underscores the importance of rigorously characterizing drivers of IAV in TBMs, not just the magnitude of IAV. TBMs that describe very different relationships with temperature and precipitation would likely yield divergent projections of future carbon budgets under climate change (e.g. Huntzinger et al 2017).