Forests dominate the interannual variability of the North American carbon sink

Understanding what drives the interannual variability (IAV) of the land carbon sink is crucial for improving future predictions of this important, yet uncertain, component of the climate system. While drivers of global and hemispheric-scale net ecosystem exchange (NEE) IAV have been investigated, our understanding of the drivers of NEE IAV at regional scales (e.g. sub-continental, biome-level) is quite poor. Here we explore the biome-level attribution and drivers of North American NEE using inverse estimates derived from a dense network of atmospheric CO2 observations. We find that deciduous broadleaf and mixed forests are the primary regions responsible for North American NEE IAV, which differs from the ecoregions identified for the globe and Northern Hemisphere. We also find that a suite of terrestrial biosphere models (TBMs) do not agree on the dominant biome contributing to NEE IAV, with TBMs falling along an apparent spectrum ranging between those with IAV dominated primarily by forested ecosystems to those with IAV dominated by non-forested ecosystems. Furthermore, this regional trade-off in TBM NEE IAV is found to be linked to differing regional responses to environmental drivers among TBMs. This work displays the importance of extra-tropical forests in driving continental NEE IAV and also highlights the challenges and limitations of using TBMs to inform regional-scale carbon flux dynamics.


Introduction
The terrestrial biosphere absorbs approximately a quarter of anthropogenic carbon dioxide (CO 2 ) emissions, thereby acting as a natural buffer against rising atmospheric CO 2 levels and the resulting changes in climate (IPCC 2013). This net sink of CO 2 can vary considerably from year to year, however, from over double the long term global average to becoming a net source of CO 2 to the atmosphere (Le Quéré et al 2016). Understanding the contemporary interannual variability (IAV) in the land sink is a crucial step toward better predicting the future buffering capabilities of the terrestrial biosphere and the subsequent fate of the climate system. The uncertainty in estimates of the future land sink by state-of-the-science models remains too large to discern whether, at a global scale, lands will be a net source or sink of CO 2 in the future (Friedlingstein et al 2006(Friedlingstein et al , 2014. The uncertainty in carbon flux dynamics at finer regional-to-continental scales has been even harder to address (Schimel 2007). Understanding regional-to-continental scale carbon cycle-climate responses is arguably even more critical as they offer relevant information for land management and climate-change mitigation decisions (Michalak et al 2011) and are key for improving the representation of climate controls in carbon cycle model development (Schimel 2007).
North America plays an integral role in the global carbon cycle, representing a quarter of the global land sink (Le Quéré et al 2016). North America spans primarily boreal and temperate zones and consists of a variety of biomes, including evergreen needleleaf forests, deciduous broadleaf forests, mixed forests, grasslands, croplands, and shrublands. At the plot level (∼1 km 2 ), flux tower studies have indicated that air temperature, precipitation, solar radiation, and disturbance ( using TBMs and inverse models, yet they have not examined sub-continental or regional (down to ∼100 km 2 ) NEE IAV nor have they explored regional attribution and drivers of NEE IAV. While TBMs can provide estimates of carbon fluxes at regional scales and have been the primary tool for exploring regional carbon flux IAV dynamics (Ahlström et al 2015, Jung et al 2017, TBMs have been shown to perform poorly in terms of predicting IAV at the plot level ∼ 1 km 2 (Keenan et al 2012, Raczka et al 2013). Even at larger scales, the use of TBMs to examine IAV could be problematic as a given TBM may be simulating historical IAV accurately for the 'wrong reasons,' i.e. with differing responses to environmental conditions (Fang et al 2017, Huntzinger et al 2017. Hence, a need exists to explore North American NEE IAV at the intermediate or regional scale with a suitable large scale observational constraint.
One reason for the gap in understanding IAV at intermediate scales is the difficulty in estimating carbon fluxes at these scales. This gap in understanding is not specific to North America, but rather emblematic of regional scale IAV, which is generally poorly understood (Niu et al 2017). Analyzing multiple flux towers within a given biome has emerged as a key tool for providing information on carbon flux dynamics representative of a specific biome or ecosystem (e.g. Baldocchi et al 2017). However, due to the limited footprint (∼ 1 km 2 ) of individual flux towers, interpolation and extrapolation is needed to examine the integrated carbon flux response of a region or biome, which can pose its own set of challenges (e.g. Hoffman et al , Kumar et al 2016. TBMs can provide estimates of carbon fluxes at intermediate scales, although these estimates are hard to validate due to the lack of an observational constraint at compatible scales. Regional inverse models leverage atmospheric CO 2 observations and an atmospheric transport model to investigate surface fluxes at regional scales. Atmospheric data have been used to infer regional carbon fluxes through inverse modeling for a variety of studies (e.g. Gourdji et al 2012, Lauvaux et al 2012, Miller et al 2013, Schuh et al 2013, Alden et al 2016, yet these studies have typically focused on just one or two years. This limitation is primarily related to the relatively recent development of more dense atmospheric observation networks, which previously were only supported via short-term measurement campaigns (e.g. Miles et al 2013). In North America, the number of continuous in-situ CO 2 measurement sites has expanded from nine sites in 2004 to ∼35 by 2007 (Shiga et al 2013, ObsPack 2016. Thus, the dense multi-year atmospheric CO 2 monitoring network across North America offers an observational constraint that enables the investigation of regional carbon flux IAV.
The main aim of this study is to explore the regional attribution and drivers of IAV in the annual North American carbon sink from a top-down atmosphericdata perspective. We estimate North American CO 2 NEE for the period 2007-2012 using a network of atmospheric CO 2 observations and a geostatistical inversion, quantify NEE IAV at the continental and biome scales, and identify the biome-level contributions to NEE IAV. We compare results from the atmospheric inversion to estimates from an ensemble of TBMs to explore and highlight any difference between these two approaches. Finally, we investigate possible reasons for the large spread in TBM based IAV contributions.

Methods and data
To analyze flux estimates at the regional-scale, we implement an atmospheric inverse approach specifically designed to account for the challenges of a regional analysis. Peylin et al (2013) list several recommendations for exploring regional carbon flux dynamics. We address three pertinent topics here: (1) the network of observations must be sufficiently dense and constant, (2) fluxes must be estimated at relatively high resolution, and (3) a realistic awareness of the impact of prior fluxes must be taken into account. The set-up of the regional inversion utilized in this study addresses each of these goals.
Six years (2007-2012) of continuous atmospheric CO 2 concentration observations are used from a network of towers across North America. CO 2 concentrations are averaged over a three hour interval centered on 3 pm and filtered in the same manner as Gourdji et al (2012) and Fang and Michalak (2015) to account for various data anomalies (e.g. outliers, transport model issues, +/−30 ppm spikes from the background). CO 2 concentration data are obtained from the ObsPack GLOBALVIEWplus v2.1 product (Masarie et al 2014, ObsPack 2016) with the addition of the site at Har-vard Forest, Massachusetts (HFM) (Urbanski et al 2007) and several sites that were part of the Mid-Continent Intensive project (Miles et al , 2013 located in Austin Cary Memorial Forest, Gainesville, Florida (AAC); Chestnut Ridge, Tennessee (ACR); Canaan Valley, West Virginia (ACV); Mead, Nebraska (AME); Missouri Ozark, Missouri (AOZ) (Stephens et al 2011); Centerville, Iowa (RCE); Galesville, Wisconsin (RGV); Kewanee, Illinois (RKW); Mead, Nebraska (RMM); Round Lake, Minnesota (RRL) ; and Rosemount, Minnesota (KCMP) (Griffis et al 2008). Over the six years, 21 sites (∼60% of the data) have a 6 year record, seven sites (∼15% of the data) have a 5 year record, four sites (∼9% of the data) have a 4 year record, nine sites (∼10% of the data) have a 3 year record, and 14 sites (∼6% of the data) have a 1 or 2 year record (for additional site information see table S1 available at stacks.iop.org/ERL/13/084015/mmedia). While the network has fluctuated over time, a majority of the sites and ∼75% of the data (Andrews et al 2014, ObsPack 2016) span at least 5 years. Additionally, the integrated sensitivity or footprint of the network covers a relatively consistent area from year to year (figure S1) and sensitivity tests (supplementary materials) showed results were not impacted by the network composition.
We use the Stochastic Time Inverted Lagrangian Transport model (STILT) (Lin et al 2003) driven by winds from the Weather Research and Forecasting model (WRF) (Skamarock andKlemp 2008, Nehrkorn et al 2010) using an inner 10 km grid over the continental US (∼25 • -55 • N; 135 • -65 • W) and an outer 30 km grid surrounding the continent (∼10 • -80 • N; 170 • -50 • W) to derive the footprint or sensitivity [ppm/( mol m −2 s −1 )] of each observation to upwind surfaces fluxes at a spatial resolution of 1 • × 1 • and a 3 hourly temporal resolution. These footprints are derived from a footprint library available as a part of the NOAA CarbonTracker Lagrange regional inversion framework (www.esrl.noaa.gov/gmd/ccgg/ carbontracker-lagrange/). NEE fluxes are estimated at a relatively high spatial and temporal resolution (1 • × 1 • in space, 3 hourly in time) for each individual year using the geostatistical inverse modeling (GIM) formulation (Michalak et al 2004, Gourdji et al 2012 and the WRF-STILT transport model. GIM incorporates covariates that are selected based on their skill in explaining the atmospheric data (Gourdji et al 2008, Fang and Michalak 2015, eliminating the need for a traditional prior flux estimate based on a TBM. The selection of covariates is informed by a statistical model selection framework . The GIM approach thus removes the dependence of the prior term on any given TBM and, when combined with the dense atmospheric network over North America and high resolution inverse framework, provides a unique platform for investigating continental-scale IAV. Carbon flux output from ten TBMs participating in the Multi-scale Terrestrial Model Intercomparison Project (MsTMIP) (Huntzinger et al 2013) are used to compare with the inverse estimates of NEE. MsTMIP uses a consistent set of environmental drivers (Wei et al 2014a(Wei et al , 2014b for all models as well as a uniform protocol and experimental design so as to focus on structural differences between models. The ten models with output scaled to the monthly temporal and 1 • × 1 • spatial scale used here mirror those used by Fang et al (2017)  The analysis uses the contribution framework developed by Ahlström et al (2015) to explore the portion of North American IAV attributable to a specific region.
where f ij represents the flux anomaly in region j and year t, and F is the North American continental-scale flux anomaly for year t. This formulation weighs a given region's flux anomaly by the similarity to the large-scale North American flux anomaly, thus defining the contribution, c , of a given region j to North American IAV. This framing enables the identification and comparison of the relative importance of a region in driving continental-scale North American IAV. We further calculate the biome-level response or sensitivity of a given flux component annual anomalies to seasonal temperature, precipitation, or radiation anomalies. This is done by first calculating the seasonal anomalies at each 1 • × 1 • grid cell for each driver variable (e.g. temperature, precipitation) as well as the annual flux anomaly. Next, for each biome (defined using a modified IGBP landcover classification map from Frankenberg et al (2011), figure 2(b)), the slope of the linear regression line between the grid cell flux anomalies and the corresponding grid cell anomalies for the given environmental variable is used to quantify the sensitivity. The temperature and precipitation data used here are the same CRU-NCEP data used to drive the MsTMIP simulations covering 1981-2010 (Wei et al 2014a(Wei et al , 2014b (2015)). The continental level IAV found here is comparable in magnitude to the total annual North American net carbon sink for 1990-2009, −0.47 +/− 0.28 PgC yr −1 (King et al 2015), where here and throughout the paper a negative flux represents a net uptake of carbon. Additionally, the average IAV from the 10 MsTMIP TBMs at the continental scale ∼0.20 PgC yr −1 (table S2) is approximately 60% of the IAV estimated here using an atmospheric constraint, which is also consistent with King et al (2015) who found that TBMs exhibited less IAV than atmospheric observation-based estimates.
The biome-level estimates of IAV derived from atmospheric observations are 1.6-4.4 times larger than the average of TBM-based estimates (figure 1, table S2). Previous studies have suggested that TBMs underestimate IAV at the site level (Keenan et al 2012, Raczka et al 2013) and at the continental level (King et al 2015) due to missing or poorly represented processes in TBMs. Here we show that, at the biome level, TBMs on average underestimate IAV in Deciduous broadleaf & mixed forests and Shrublands, by factors of 3.1 and 4.4, respectively. In these two biomes, the largest IAV predicted by a TBM is less than 60% of the corresponding atmospheric observation-based estimate of IAV.
We identify Deciduous broadleaf & mixed forests as the largest overall contributors (∼40%) to North American NEE IAV (figure 2), the largest contributor per unit area, and the biome with the largest NEE IAV magnitude. This is in contrast to the global scale, where Ahlström et al (2015) found that semi-arid ecosystems are the dominant contributor to IAV (39%), according to an ensemble of TBM simulations from the TRENDY model intercomparison project (Sitch et al 2015). Similarly, for the northern hemisphere, Fu et al (2017a) found semi-arid ecosystems (39%) and Grasslands and Croplands (39%) to be the dominant contributors to IAV, based on the Monitoring Atmospheric Composition and Climate (MACC-II) inversion system (www.ecmwf.int/en/research/projects/macc-ii). While Deciduous broadleaf forests have been identified at the site level as an ecosystem with relatively large IAV when compared to other North American biomes (Yuan et al 2009, Baldocchi et al 2017, our results demonstrate that at the North American continental scale these ecosystems are dominating IAV. Addi-tionally, we find that TBMs appear to consistently underestimate the contribution of IAV from Deciduous broadleaf & mixed forests, which may be linked to the documented underestimation of IAV across forested sites by TBMs (Keenan et al 2012, Raczka et al 2013).
The Shrublands biome, while having the second largest overall IAV magnitude (Table S2) shows the lowest relative contribution to North American IAV (figure 2). This is due to the metric defined by Ahlström et al (2015), equation (1), which quantifies the contribution based both on the magnitude of and the 'synchronicity' between regional anomalies and the 'global' anomaly (in this case continental anomaly). Put another way, although Shrublands exhibit high IAV, there is a difference in phasing between Shrublands and the rest of continental North American IAV. Likely Shrublands, which represent 60% of the North American semi-arid ecosystems by area, are more responsive to global drivers of semi-arid ecosystems (e.g. drought (Huang et al 2016)) rather than drivers of the forested biomes that make up a large portion of continental North America. Therefore, these findings may indeed be consistent with previous findings recognizing the role of semi-arid regions on global IAV (Poulter et al 2014, Ahlström et al 2015. Interestingly, the MsTMIP TBMs also show a low contribution of the Shrublands biome to North American IAV, however the low IAV contribution in TBMs is largely due to a smaller magnitude of Shrublands IAV, rather than the phasing difference. We find the spread across IAV contributions from the various TBMs in the MsTMIP ensemble to be quite large (figure 2). Consequently, no dominant biome can be clearly identified for the MsTMIP TBMs. The large spread in biome-level IAV contributions in TBMs can in part be linked to an apparent trade-off between forested (Needleleaf and Deciduous broadleaf & mixed forests) vs. non-forested (Grasslands and Shrublands) biomes (figure 3). For the following analysis, we group the biomes into these two broader categories. This is done based on the notion that the response of TBMs to the climatic controls related to drought stress (e.g. temperature, water availability) in these two broad categories ( TBM responses to anomalies in summertime temperature and precipitation help to explain the IAV trade-off between forested and non-forested regions (figure 3) (similar responses observed for radiation, figure S2 and S3). TBMs with stronger NEE sensitivities to temperature and precipitation in forested biomes (relative to their sensitivities in non-forested biomes) predict that forested regions dominate North American NEE IAV, while TBMs with stronger sensitivities in non-forested biomes have higher NEE IAV contributions from non-forested biomes ( figure 3). This result shows that, for some TBMs, there are clear differ-ences between how forested and non-forested biomes respond to environmental driver anomalies (figure 4), but that there is little agreement about which regions are more or less sensitive to temperature and precipitation anomalies. Hence, the environmental driver sensitivity differences appear to impact the biome-level attribution of NEE IAV in a given TBM, implying a potential indirect tuning effect where a TBM either has a strong NEE sensitivity in forested biomes or a strong NEE  (d)), for forested ecosystems the atmospheric-observation based sensitivities of NEE anomalies appear to be a lower end-member relative to most of the TBMs (figures 4(a) and (c)). One explanation for this is that TBMs may be overly sensitive to summertime environmental driver anomalies in forests. Annual NEE anomalies in forested ecosystems are often more complex than a simple relationship with growing season environmental conditions and can involve a variety of multifaceted responses (e.g. counteracting seasonal anomalies (e.g. As NEE is ultimately the result of two opposing processes, gross primary production (GPP) and ecosystem respiration (RE), we further probed the TBMs to investigate which component drives the trade-off in   relative IAV contributions between the forested and non-forested regions. The correspondence of GPP sensitivity differences with regional NEE IAV contribution differences (figures S3(d)-(f)) is in contrast to the lack of a clear relationship in the corresponding RE sensitivity difference (figure S3(g)-(i)). Thus, we conclude that the regional IAV trade-off is primarily explained by GPP, rather than RE sensitivity differences. This may be somewhat expected as summertime NEE dynamics in the biomes examined are largely driven by GPP, however it suggests that at the component level TBM formulation/parameterization appears to be linked to regional IAV dynamics.

Conclusions
Overall, this work explores North American NEE IAV dynamics at scales that bridge the gap between the site-level and global-level using a dense network of atmospheric CO 2 measurements, which provides a top-down regionally sensitive observational constraint. These findings address one of the fundamental questions regarding the terrestrial carbon cycle, namely quantifying the variability of carbon fluxes from year to year. The regional atmospheric observationconstrained estimates show that North America has larger variability in NEE than TBMs predict, at both the biome and continental level. We also identify Decidu-ous broadleaf & mixed forests as the main sources of North American NEE IAV, which differs from the main global and northern hemisphere drivers. These findings underscore the importance of North American extratropical forests in regulating continental carbon flux variability, and stress the critical need to improve our understanding of the processes driving regional NEE IAV. Across North America, we find that current stateof-the-art TBMs disagree about the main regions responsible for NEE IAV, as well as about the NEE, GPP, and RE responses to climate drivers. In TBMs, a trade-off between whether forested or non-forested regions dominate continental IAV is found to be linked to a given TBM's sensitivity to environmental drivers. TBMs with a stronger sensitivity of NEE to environmental drivers in forested (relative to non-forested) regions have a larger contribution of NEE IAV from forested regions and vice versa. The regional differences in NEE responses to environmental drivers in TBMs has important implications for future predictions of the regional and biome-level behavior of the terrestrial biosphere as temperature and precipitation regimes shift and as extreme events become more likely (Seneviratne et al 2012, Reichstein et al 2013. The ensemble of TBMs examined here, for example, would likely yield divergent predictions about the fate of the carbon sink in forested vs. non-forested regions of North America under changing temperature and precipitation conditions. Hence, more research is needed to understand and observationally constrain the carbon cycle response to climatic conditions at the regional scale. We acknowledge EUMETSAT for the GOME-2 data and Joanna Joiner for the GOME-2 SIF data originally funded in part by NASA Carbon Cycle Science program grant number NNH10DA001N. We gratefully acknowledge the efforts of the PIs of the various towers providing continuous atmospheric CO 2 observations, which were instrumental for these analyses.