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

Year-to-year variability in CO2 fluxes can yield insight into climate-carbon cycle relationships, a fundamental yet uncertain aspect of the terrestrial carbon cycle. In this study, we use global observations from NASA’s Orbiting Carbon Observatory-2 (OCO-2) satellite for years 2015–2019 and a geostatistical inverse model to evaluate 5 years of interannual variability (IAV) in CO2 fluxes and its relationships with environmental drivers. OCO-2 launched in late 2014, and we specifically evaluate IAV during the time period when OCO-2 observations are available. We then compare inferences from OCO-2 with state-of-the-art process-based models (terrestrial biosphere model, TBMs). Results from OCO-2 suggest that the tropical grasslands biome (including grasslands, savanna, and agricultural lands within the tropics) makes contributions to global IAV during the 5 year study period that are comparable to tropical forests, a result that differs from a majority of TBMs. Furthermore, existing studies disagree on the environmental variables that drive IAV during this time period, and the analysis using OCO-2 suggests that both temperature and precipitation make comparable contributions. TBMs, by contrast, tend to estimate larger IAV during this time and usually estimate larger relative contributions from the extra-tropics. With that said, TBMs show little consensus on both the magnitude and the contributions of different regions to IAV. We further find that TBMs show a wide range of responses on the relationships of CO2 fluxes with annual anomalies in temperature and precipitation, and these relationships across most of the TBMs have a larger magnitude than inferred from OCO-2. Overall, the findings of this study highlight large uncertainties in process-based estimates of IAV during recent years and provide an avenue for evaluating these processes against inferences from OCO-2.

inferences from OCO-2 with state-of-the-art process-based models (terrestrial biosphere model, TBMs). Results from OCO-2 suggest that the tropical grasslands biome (including grasslands, savanna, and agricultural lands within the tropics) makes contributions to global IAV during the 5 year study period that are comparable to tropical forests, a result that differs from a majority of TBMs. Furthermore, existing studies disagree on the environmental variables that drive IAV during this time period, and the analysis using OCO-2 suggests that both temperature and precipitation make comparable contributions. TBMs, by contrast, tend to estimate larger IAV during this time and usually estimate larger relative contributions from the extra-tropics. With that said, TBMs show little consensus on both the magnitude and the contributions of different regions to IAV. We further find that TBMs show a wide range of responses on the relationships of CO 2 fluxes with annual anomalies in temperature and precipitation, and these relationships across most of the TBMs have a larger magnitude than inferred from OCO-2. Overall, the findings of this study highlight large uncertainties in process-based estimates of IAV during recent years and provide an avenue for evaluating these processes against inferences from OCO-2.

Introduction
Interannual variability (IAV) in CO 2 fluxes is of critical importance in understanding the global carbon cycle. The magnitude of IAV in global terrestrial fluxes is large, comparable to that of land fluxes in an average year (Peylin et al 2013). An investigation of IAV also offers an opportunity to explore the sensitivity of the carbon cycle to changing environmental conditions and may therefore yield a better understanding of climate-carbon cycle relationships. Insight into contemporary IAV may therefore help inform projections of future CO 2 budgets under climate change (e.g. Cox et al 2013, Friedlingstein et al 2014, Hoffman et al 2014. IAV in CO 2 fluxes at global scales is relatively well known; however, the contributions of different ecoregions to global IAV remain highly uncertain (e.g. Baker et al 2006, Peylin et al 2013, Deng et al 2014. In situ observations of atmospheric CO 2 have been extensively used to estimate global-and regional-scale IAV (e.g. Bousquet et al 2000, Gurney et al 2003, Rödenbeck et al 2003, Bruhwiler et al 2011, Peylin et al 2013, Shiga et al 2018, Hu et al 2019, Keppel-Aleks et al 2014. However, in situ measurement sites are unevenly distributed across the globe; therefore, inverse models using in situ observations are arguably not sensitive to CO 2 fluxes in undersampled regions like the tropics, where IAV is the largest (e.g. Baker et al 2006, Peylin et al 2013, Piao et al 2020. The advent of space-based observations of CO 2 provides greater coverage than in situ measurements and opens new window into the study of IAV, though these observations are available for a more limited time period (e.g. Houweling et al 2004, Chevallier et al 2007. For example, the Greenhouse Gases Observing Satellite (GOSAT), launched in 2009, is the first Earth-orbiting satellite that is designed to collect atmospheric CO 2 with sufficient accuracy to estimate surface CO 2 fluxes (Yokota et al 2009). Existing studies have leveraged GOSAT observations to estimate IAV (e.g. Guerlet et al 2013, Byrne et al 2019, Liu et al 2021. However, in a recent study, Byrne et al (2019) argued that current GOSAT observations can only be used to constrain IAV across continent-sized regions or larger spatial scales.
It is even more challenging to understand the relationships between IAV and underlying environmental drivers across different regions of the globe. Here we define 'environmental drivers' as meteorological variables or characteristics of the physical environment that may correlate with net ecosystem exchange (NEE) and can be quantified by measurements or models (e.g. precipitation, air temperature). Specifically, at small, local scales (∼1 km 2 ), eddy flux tower measurements have been used as one of the most important approaches in quantifying IAV and its relationship with underlying environmental drivers (e.g. Baldocchi et al 2001, Barford et al 2001, Suyker et al 2001, Schwalm et al 2007, Jensen et al 2017. Due to the limited footprints of flux towers (∼ 1 km 2 ), extrapolation of local eddy covariance measurements in space and time is necessary to obtain larger, regionalto global-scale estimate of CO 2 fluxes (e.g. Beer et al 2010, Jung et al 2011, Tramontana et al 2016. This upscaling can be challenging to do, and recent studies have argued that upscaling efforts tend to underestimate IAV relative to atmospheric inverse models (e.g. Byrne et al 2019, Jung et al 2020, Piao et al 2020. Furthermore, process-based terrestrial biosphere models (TBMs) have been widely used to simulate net carbon exchange (NEE) and IAV, and these models provide a means to attribute changes in net carbon uptake to specific environmental drivers (e.g. Tarnocai et al 2009, Medvigy et al 2010, Belshe et al 2013. With that said, TBMs do not show consensus on the magnitude of IAV in different ecoregions of the globe (e.g. Shiga et al 2018). Even when TBMs do agree on the magnitude of IAV, they often yield similar patterns for different reasons, i.e. TBMs display very different sensitivities of IAV to environmental drivers (Piao et al 2013, Huntzinger et al 2017. NASA's Orbiting Carbon Observatory 2 (OCO-2) satellite has the potential to provide additional, new information to investigate IAV and its relationships with environmental drivers (e.g. Eldering et al 2017). OCO-2 was launched in late 2014, so it does not provide as long of a record as many in situ or eddy flux measurement sites. With that said, its unprecedented global coverage and density of observations provides a new opportunity to examine IAV, albeit across a shorter window of time. Indeed, previous studies have used the global set of OCO-2 observations to estimate regional and global CO 2 fluxes (e.g. In this study, we employ OCO-2 observations (years 2015-2019) and a geostatistical inverse model (GIM) to examine the contributions of different regions to global IAV in CO 2 fluxes across the 5 year study period. We also explore the relationships between this IAV and anomalies in environmental drivers across the globe. We first present an overview of the GIM approach. We then present results on 5 years of IAV for different global biomes and evaluate the relationships between recent IAV and

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

Geostatistical inverse model
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 2 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: The fluxes s (dimensions m × 1), when run through an atmospheric model, h (s), should match the observations z (dimensions n × 1) within a specific error ϵ (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': 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, Xβ is the 'determinstic model' . The unknown vector ζ (dimensions m × 1) is the stochastic component and is estimated at the model grid scale. The posterior flux estimate ( s, dimensions m × 1) is a sum of the deterministic model (Xβ) and the stochastic components (ζ). We simultaneously estimate the fluxes (s) and the coefficients (β) via minimizing the GIM cost function (e.g. Michalak et al 2004): The cost function includes two covariance matrices; R (dimensions n × n) and Q (dimensions m × m). The covariance matrix R describes z − h (s), 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 (ζ) 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 5 ) to estimate ∼6.0 × 10 6 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 ( s) using a reduced rank algorithm described in Saibaba and Kitanidis (2015) and , an approach described in detail in supporting information text S1.2.
We further compare the GIM flux estimates against aircraft-based observations of CO 2 and against CO 2 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.

Model selection
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 (β) 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 2 fluxes for seven global biomes, and we therefore use seven biomes in this study.
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 . 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) Figure 2. IAV during the 5 year study period from a suite of 16 TBMs in TRENDY (blue), from the ensemble mean of TBMs (black), and from the estimate using OCO-2 (red). Error bars indicate associated uncertainties (one standard deviation; see supporting information text S6). Most TBMs exhibit larger IAV compared to the estimate using OCO-2. Furthermore, a large number of TBMs (nine out of 16) estimate the largest IAV from tropical forests, whereas the findings using OCO-2 suggest that the tropical grasslands and tropical forests biomes make comparable contributions to global IAV, at least during the 5 year study period. and inverse modelling studies (e.g. Gourdji et al 2012, Miller et al 2013, Fang and Michalak 2015, Miller et al 2016. 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.

State-of-the-art process-based models
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 2 , climate, and land use; NBP from TRENDY represents net CO 2 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.

The contributions of different regions to global IAV during the study period
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 2 , at least during the time period for which satellite observations are available.
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 from the estimate using OCO-2 (red). In each panel, we order the models from the smallest contribution to the largest contribution. Compared to the estimate using OCO-2, most TBMs (12 out of 16) estimate lower relative contributions from tropical grasslands to global IAV during the study period, and approximately half of TBMs display lower relative contributions from tropical forests. By contrast, most TBMs suggest higher relative contributions from the extra-tropics to global IAV. 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.
Specifically, the estimate using OCO-2 suggests that the tropical grasslands and tropical forests biomes make comparable contributions to global IAV 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  explored drought impact on CO 2 fluxes over the contiguous U.S. (CONUS), and suggested that the regional flux anomaly (departure from a 6 year mean) over the CONUS droughtimpacted regions are in opposite directions of the global atmospheric CO 2 growth rate anomaly. The overall contribution of the extra-tropics is a sum of both positive and negative contributions of smallerscale 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.

Connections between IAV during years 2015-2019 and environmental drivers
There is a growing consensus that tropical regions are the largest contributor to global IAV (e.g. Bousquet  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 2 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 2 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. between years 2015 and 2017 contribute most to the grid-scale differences in CO 2 fluxes between the 2 years ( figure 4(b)). By contrast, the stochastic component ζ (figure 4(f)) describes relatively little yearto-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.

Relationships with environmental drivers help explain differences between OCO-2 and TBMs
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).
We find that CO 2 fluxes in most TBMs tend to have stronger relationship with annual anomalies in Figure 5. CO2 fluxes in most TBMs (blue) show stronger relationships with annual anomalies in temperature and precipitation than inferred from OCO-2 (red) (a)-(d). This conclusion is particularly evident if we exclude models that show low correlation between anomalies in environmental drivers and flux anomalies (R 2 < 0.1), denoted in light blue above. In each panel, we order the models from the weakest relationship to the strongest relationship. For each TBM and the OCO-2 estimate, we fit a linear regression model between flux anomalies (departure from the 5 year mean) and annual anomalies in environmental drivers using model grid-scale datasets within each biome, and the slope of the linear regression model therefore indicates the relationships with the environmental drivers. Here we show the estimated slopes for various TBMs and OCO-2. Note the climate forcing data used in TRENDY models are from the CRUJRA datasets (see section 2.4), and we use environmental drivers from the CRUJRA data to examine the relationships here for these TBMs. 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).

Conclusions
In this study, we employ 5 years of OCO-2 satellite observations and a GIM to evaluate year-to-year variability in the global carbon cycle and its relationships with environmental drivers; we then compare IAV inferred from OCO-2 against 16 state-of-the-art TBMs.
This analysis provides an avenue for evaluating IAV in process-based flux estimates using satellite observations of CO 2 . Despite the global observational coverage provided by OCO-2, the inverse modeling results presented here show substantial uncertainties (e.g. figure 2), pointing to the limitations of current space-based CO 2 monitoring. With that said, the uncertainties in inverse model are still smaller than the large range of estimates from the TBMs, implying that satellite observations can provide an important tool to evaluate IAV and even the environmental factors that correlate with this IAV. Furthermore, a longer record of satellite observations from OCO-2 (or comparable instruments) in the future would facilitate an examination of longer time periods, and this study is a first step into exploring IAV using OCO-2 satellite observations.
Overall, we find that state-of-the-art TBMs show little consensus on the magnitude of IAV at either regional or global scales, at least for the 5 year period examined in this study. Compared to these models, our analysis using OCO-2 indicates (a) a larger contribution of tropical grasslands, savanna, and agricultural lands to global IAV during the study period, (b) a smaller contribution from the extra-tropics, and (c) a slightly smaller overall magnitude of IAV. These differences, particularly across the tropics, appear to stem from disagreement in the relationships within each model between carbon fluxes and annual anomalies in temperature and precipitation. This study hence reinforces the need to rigorously characterize the relationships between environmental drivers and IAV within TBMs, not just the estimated magnitude.

Data availability
Data information of the ObsPack data product is available at www.esrl.noaa.gov/gmd/ccgg/obspack/; data information of the TRENDY v8 is available at http://sites.exeter.acuk/trendy/.
The data that support the findings of this study are openly available at the following URL/DOI: ftp:// ftp.cira.colostate.edu/ftp/BAKER/.