Multiscale assessment of North American terrestrial carbon balance

Comparisons of carbon uptake estimates from bottom-up terrestrial biosphere models (TBMs) to top-down atmospheric inversions help assess how well we understand carbon dioxide (CO₂) exchange between the atmosphere and terrestrial biosphere. Previous comparisons have shown varying levels of agreement between bottom-up and top-down approaches, but they have almost exclusively focused on large, aggregated scales (e.g., global or continental), providing limited insights into reasons for the mismatches. Here we explore how consistency, defined as the spread in net ecosystem exchange (NEE) estimates within an ensemble of TBMs or inversions, varies with at finer spatial scales ranging from $\mathrm{1}{}^{\circ }\times \mathrm{1}{}^{\circ }$ to the continent of North America. We also evaluate how well consistency informs accuracy in overall NEE estimates by filtering models based on their agreement with the variability, magnitude, and seasonality in observed atmospheric CO₂ drawdowns or enhancements. We find that TBMs produce more consistent estimates of NEE for most regions and at most scales relative to inversions. Filtering models using atmospheric CO₂ metrics causes ensemble spread to decrease substantially for TBMs, but not for inversions. This suggests that ensemble spread is likely not a reliable measure of the uncertainty associated with the North American carbon balance at any spatial scale. Promisingly, applying atmospheric CO₂ metrics leads to a set of models with converging flux estimates across TBMs and inversions. Overall, we show that multiscale assessment of the agreement between bottom-up and top-down NEE estimates, aided by regional-scale observational constraints is a promising path towards identifying fine-scale sources of uncertainty and improving both ensemble consistency and accuracy. These findings help refine our understanding of biospheric carbon balance, particularly at scales relevant for informing regional carbon-climate feedbacks.

1 Introduction

Reliable estimates of carbon dioxide (CO₂) uptake by the terrestrial biosphere are necessary for understanding both historical and future climate change because the terrestrial biosphere mitigates anthropogenic CO₂ emissions by storing carbon in above- and below-ground biomass. Net ecosystem exchange (NEE) of CO₂ cannot be measured directly at scales greater than a plot (∼ 1 km²); that is, the scale of the footprint of an eddy covariance flux tower (Kljun et al., 2015). Estimates at broader scales therefore rely either on “bottom up” methods such as terrestrial biosphere models (TBMs) that represent process-based understanding of flux drivers (e.g., Hayes et al., 2012; Sitch et al., 2015) or on “top down” methods such as atmospheric inverse models that attribute observed atmospheric variability in CO₂ concentrations to upwind biospheric activity (e.g., Ciais et al., 2011; Gourdji et al., 2012; Gurney et al., 2002; Michalak et al., 2004; Peylin et al., 2013; Shiga et al., 2018a; Thompson et al., 2016).

Understanding of net biospheric carbon uptake has improved through comparisons between TBMs and inversions (Hayes et al., 2012; King et al., 2015), but discrepancies both between these approaches and between specific models using either approach have persisted even in relatively well studied regions such as North America and Europe. The resulting uncertainties in carbon flux estimates limit our ability to anticipate carbon–climate feedbacks, and therefore to assess the impacts of alternative emission pathways and related climate mitigation policies (Friedlingstein et al., 2014; Huntzinger et al., 2017; King et al., 2015). Uncertainties both between modeling approaches and between specific models arise from the specific characteristics of TBMs and inversions.

For TBMs, model performance is limited by incomplete understanding of underlying processes (Hayes et al., 2012; Huntzinger et al., 2012, 2017; Schwalm et al., 2010, 2019; Seiler et al., 2022). One major source of uncertainty is that models may incorporate different key mechanisms or represent them with varying levels of detail. For example, nitrogen limitation likely reduces CO₂ fertilization, but coupled carbon–nitrogen dynamics are not included in all models (Bonan and Levis, 2010; Brovkin and Goll, 2015; Jain et al., 2009; Sokolov et al., 2008; Tharammal et al., 2019; Thornton et al., 2009; Wieder et al., 2015), and permafrost thaw, which can cause release of carbon from high latitudes (Burke et al., 2013; Koven et al., 2013), is not mechanistically resolved in most models. Large discrepancies also arise from different approaches to modeling land use change, vegetation dynamics, and fire (Ahlström et al., 2015; Bastos et al., 2020; Hardouin et al., 2022; Tharammal et al., 2019). Beyond such structural uncertainties, models parameterize processes differently and often use different driving data, which introduces additional uncertainty (Huntzinger et al., 2017; Jung et al., 2007; Lovenduski and Bonan, 2017; Schwalm et al., 2010, 2019). Previous studies have addressed the use of different driver data by establishing model ensembles that adhere to a standardized protocol and utilize consistent forcing data, but doing so does not fully account for observed uncertainties (Huntzinger et al., 2013, 2020; Friedlingstein et al., 2020; Sitch et al., 2015). Understanding and addressing these key sources of uncertainty in TBMs helps to improve our understanding of biospheric carbon exchange

For inversions, uncertainties arise from the limited information content of available atmospheric observations and from choices made in the statistical setup of the model. More specifically, the choice of prior estimates, prior error correlations, observational data, transport model, boundary conditions, data assimilation time period, and model resolution all lead to differences across models (Ciais et al., 2010; Göckede et al., 2010; Kondo et al., 2020; Michalak et al., 2017; Peylin et al., 2013). Paradoxically, even regions with relatively high data availability, such as Europe, can still exhibit large spread in carbon uptake estimates derived from inverse models (Kondo et al., 2020; Monteil et al., 2020). This suggests that other aspects of the inverse modeling framework may be more significant contributors to uncertainty, though there is not a clear consensus on the main source of uncertainty (Gaubert et al., 2019; Kondo et al., 2020). For example, some studies have found transport errors to be a primary source of uncertainty (Peylin et al., 2013; Schuh et al., 2019), while other studies have shown that fossil fuel emission uncertainties play a key role in variability (Gurney et al., 2005; Peylin et al., 2011; Saeki and Patra, 2017). Another limitation is that atmospheric CO₂ observations, and inversions based on these observations, do not directly inform process-level understanding of the controls on carbon uptake to the extent that bottom-up approaches can (Baker et al., 2006; Gurney et al., 2004).

Given the uncertainties inherent to both approaches and their complementary strengths, comparisons between ensembles of TBMs and inversions can be particularly helpful in diagnosing our understanding of the carbon balance of the terrestrial biosphere. If TBMs and inversions yield similar estimates, this agreement increases confidence in the reliability of both model types as sources of realistic information. Previous studies, however, have come to different conclusions about the degree of agreement between bottom-up and top-down estimates (Bastos et al., 2020; Canadell et al., 2011; Kondo et al., 2020; Hayes et al., 2018). Some comparisons have shown agreement between NEE estimates from TBMs and inversions (Ciais et al., 2010; King et al., 2015; Sitch et al., 2015), while others have shown a large amount of variability both within and between bottom-up and top-down ensembles (Bastos et al., 2020; Chevallier et al., 2014; Huntzinger et al., 2012, 2017; Schwalm et al., 2019; Sun et al., 2021). For instance, King et al. (2015) found bottom-up and top-down methods agree on the sign of the North American land sink, albeit with notable spread in estimates, while Kondo et al., found bottom-up and top-down approaches disagree in many countries at the regional scale. Furthermore, Chevallier et al. (2014) showed that an ensemble of inversions had significant disagreement at the hemispheric and regional scales and Huntzinger et al. (2012) reported a lack of consensus among TBMs in determining whether North American land functions as a carbon sink or source.

A common approach for comparing TBMs and inversions is to examine the agreement across the mean of model ensembles (Ciais et al., 2010; Hayes et al., 2012; King et al., 2015). Here we use the term “agreement” between estimates to refer to the degree to which NEE estimates from an ensemble of TBMs differ from NEE estimates from an ensemble of inversions. In addition, analyzing the variability or spread amongst TBMs or inversions provides insights into the degree of consistency in carbon flux estimates. Here we use the term “consistency” of estimates to refer to the degree of variability in NEE estimates within an ensemble of either TBMs or of inversions.

Previous comparisons between bottom-up and top-down estimates have revealed that the agreement between estimates depends on the spatial scale and the region. At the global scale, inversions yield more consistent NEE estimates than TBMs (Friedlingstein et al., 2022), largely due to the global constraint provided by atmospheric CO₂ observations. However, at regional scales, consistency is limited for both TBMs and inversions and it is even difficult to determine whether certain regions are a net sink or source of CO₂ (Ciais et al., 2013; Kondo et al., 2020) due to the uncertainties associated with both approaches outlined above. For North America, agreement between bottom-up and top-down estimates has improved over time, but this apparent agreement is in part due to the large range of estimates (i.e., low consistency) for both TBMs and inversions (Hayes et al., 2012; King et al., 2012, 2015; Pacala et al., 2001). This scale- and region-dependent agreement makes it difficult to determine the optimal path towards reducing uncertainties. This challenge is in part because comparisons of bottom-up and top-down methods are primarily conducted at large aggregated scales – for example, global, hemispheric, and continental scales (Bastos et al., 2020; Ciais et al., 2010; Hayes et al., 2012; Huntzinger et al., 2012; Pacala et al., 2007; Peylin et al., 2010), which aids little in attributing causes of observed mismatches (Bastos et al., 2020; Hayes et al., 2012; Kondo et al., 2020).

A key step forward is to look at agreement beyond large, aggregated scales. Looking across multiple spatial scales provides a more in-depth understanding of the level of agreement between carbon budgets from bottom-up and top-down approaches. Despite this, few studies have taken this approach, especially for sub-continental spatial scales. Agreement between bottom-up and top-down estimates at global and hemispheric scales (Bastos et al., 2020; Kondo et al., 2020; Sitch et al., 2015) is a necessary but not sufficient condition for reconciling differences in carbon budgets at regional scales (Kondo et al., 2020). When large spread in model estimates makes it difficult to determine whether large regions such as Europe, boreal Asia, Africa, South Asia, and Oceania are even net sinks or sources (Kondo et al., 2020), multiscale comparisons may shed new light on the lack of consistency. Gourdji et al. (2012) compared bottom-up and top-down models at sub-continental scales and found better agreement during the growing seasons than in the dormant season, allowing for a more in-depth and focused exploration into the reasons for the observed (dis)agreement. The key insights gained from multiscale comparisons highlight the need for a more comprehensive comparison of bottom-up and top-down NEE estimates across spatial scales (Gourdji et al., 2012; Kondo et al., 2020).

Examining the agreement between bottom-up and top-down methods across spatial scales can also provide insights into the relationship between consistency, agreement, and accuracy in model predictions. Assessing agreement between bottom-up and top-down budget estimates is not necessarily equivalent to determining the accuracy of carbon budgets, however (Knutti et al., 2010; Kondo et al., 2020; Lovenduski and Bonan, 2017). Instead, accuracy should be assessed against observational constraints. There have been efforts, such as the International Land Model Benchmarking Project (ILAMB), to evaluate model accuracy by quantifying agreement between reference datasets and model outputs across multiple statistical metrics (Collier et al., 2018). Model skill scores are useful in assessing agreement between reference data and model data, but it is possible to misinterpret model performance without careful analysis of the metrics that make up the overall skill score (Bonan et al., 2019; Collier et al., 2018). In addition, reference data to which models are compared are often themselves modeled data products (e.g., FLUXCOM is used as reference data for gross primary productivity (GPP); Seiler et al., 2022). Model evaluation against atmospheric CO₂ observations can in principle provide more direct insights into the variability and accuracy of model NEE estimates (Fang et al., 2014; Fang and Michalak, 2015; Sun et al., 2021). By comparing carbon budget estimates to atmospheric CO₂ observations, in addition to comparing these estimates across spatial scales, we can determine the degree to which consistency within an ensemble is representative of accuracy. While better model performance under current conditions does not necessarily indicate better performance under future conditions, this approach helps with model improvement by allowing for quick identification of which models yield more realistic results. Subsequently, in-depth analysis can be done to identify key model characteristics that lead to improved results.

Here, we compare large ensembles of bottom-up and top-down model estimates of NEE for North American across various spatial scales to assess how consistency in model estimates varies across scales and between modeling approaches (Table 1). We expect inversions to be more consistent at larger scales thanks to the constraint provided by atmospheric observations, and TBMs to be more consistent at smaller spatial scales because they are informed by process-based understanding. We then evaluate whether greater consistency corresponds to higher accuracy and lower uncertainty in overall NEE estimates. We determine if the consistency within ensembles and the agreement between top-down and bottom-up approaches are impacted by ensemble subsetting; that is, limiting ensembles to models that can reproduce basic aspects of the variability, magnitude, and seasonality of atmospheric CO₂ observations. We expect consistency to improve when ensembles include only models that agree with basic features of atmosphere CO₂ observations, thereby also increasing the degree to which consistency informs accuracy, or, in other words, making model spread a more apt measure of uncertainty.

Table 1 Names and references for the models included in each ensemble used in this study. The models that meet the variability, seasonality, magnitude, or all three metrics are indicated by an “X”. An asterisk indicates the models that were present in MsTMIP-v2 (CLM4 BG1, ISAM BG1, and VISIT SG3) but are replaced by updated versions in TRENDY-v9; in this case the metrics were evaluated based on TRENDY-v9 versions (see Sect. 2.1.1), but the MsTMIP-v2 versions were used when evaluating the impact of how NEE is defined on consistency and agreement (Sects. 3.1 and 3.3.3, Fig. S7).

2 Data and methods

2.1 Data

2.1.1 Model ensembles

Estimates of NEE from three model ensembles were used (Table 1). Bottom-up estimates came from two TBM intercomparison projects, namely, the Multi-scale Synthesis and Terrestrial Model Intercomparison Project (MsTMIP-v2; Huntzinger et al., 2013, 2018, 2021; Wei et al., 2014a, b) and the Trends in Net Land-Atmosphere Exchange version 9 ensemble (TRENDY-v9; Friedlingstein et al., 2020; Sitch et al., 2015). MsTMIP is a model intercomparison project aimed at exploring the impact of structural differences in models by prescribing a fixed protocol with a semi-factorial design and consistent environmental driver data for an ensemble of models (Huntzinger et al., 2013, 2018, 2021; Wei et al., 2014a, b). We use the SG3 and BG1 simulations from MsTMIP-v2, where atmospheric CO₂ and land-use history are time-varying and nitrogen deposition rates are held constant in SG3, and land-use history, atmospheric CO₂, and nitrogen deposition are time-varying in BG1 (Table 1). MsTMIP-v2 defines NEE to be NEE = R_h+ R_a+ F_disturbance+ F_product− GPP where R_h is heterotrophic respiration, R_a is autotrophic respiration, F_disturbance is the sum of fire and land use change fluxes, and F_product is the decay of harvested wood products. Some models do not include fire, land use, and/or product fluxes, depending on the processes that are represented by the models. TRENDY includes simulations from an ensemble of dynamic global vegetation models (DGVMs) that are run annually as part of the Global Carbon Project yearly evaluation. TRENDY also includes a set of factorial simulations for the historical period (Friedlingstein et al., 2022). Here we use the S3 simulation from TRENDY-v9 where CO₂, climate, and land use forcings are time-varying. TRENDY defines NEE as NEE = R_h+ R_a+ F_fire+ F_harvest+ F_grazing+ F_LUC – GPP, where F_fire are emissions from natural and human-sparked fires, F_harvest is emissions from crop harvest, F_grazing are fluxes from livestock grazing, and F_LUC are fluxes resulting from land use change. Models within TRENDY-v9 vary in terms of which components are included. We detail the impact of varying NEE definitions in bottom-up models on our results in Sect. 3.1. Top-down estimates come from a set of inverse model estimates assembled in support of the REgional Carbon Cycle Assessment and Processes 2 (RECCAP-2) analysis (Ciais et al., 2022). RECCAP-2 is a project aimed at quantifying carbon budgets on regional scales across the globe.

For assessing consistency and agreement at biomes scales we use the same biome map as Shiga et al. (2018b) and Sun et al. (2021) that is based on an International Geosphere-Biosphere Programme (IGBP) land cover classification map. MsTMIP-v2 imposes a consistent biome map for all models (Wei et al., 2014a), while models from TRENDY-v9 use various sets of plant functional types (PFTs), resulting in differences in land cover representations (Seiler et al., 2022). While the TRENDY-v9 models have differences in land cover representation, they do use common land use and land cover change (LULCC) forcing data (Seiler et al., 2022). Though using the same biome map across all models would allow for greater standardization amongst TBMs, doing so is difficult due to model-specific setups. Understanding the impact of the various maps used by models is also difficult as few models provide outputs at the resolution necessary to do a comprehensive analysis, such as evaluating whether specific PFTs are present at in situ observation sites (Seiler et al., 2022). However, Sun et al. (2021) did compare the impact of model-specific biome classifications for four models that provided PFT information at finer resolutions and showed that model-specific biome classification was not a primary driver of inter-model variability.

All models were re-gridded to a $\mathrm{1}{}^{\circ }\times \mathrm{1}{}^{\circ }$ spatial and monthly temporal resolution and then cropped to a uniform North American domain for the study period of 2007–2010. We correct the CLM5.0 timestamp to address an error causing output files to have a start date shifted by 1 month (ESDS 0.1 documentation FAQ, 2024). For TBMs, only land fluxes are included in model output files and thus ocean grid cells are removed prior to regridding. This period was chosen because this was the time frame for which all models overlapped temporally (until 2010) and during which high-resolution atmospheric transport footprints (see Sect. 2.3) were available (starting in 2007) to link model estimates to atmospheric CO₂ observations. We use the model ensemble average when assessing agreement between bottom-up and top-down model ensembles. While this is a simple approach, Schwalm et al. (2015) found that the added complexity of skill-based integration does not materially change flux estimates based on TBM ensembles. The model simulations from MsTMIP-v2 and TRENDY-v9 ensembles were merged to create one TBM superensemble. Because MsTMIP-v2 and TRENDY-v9 have three models in common (CLM, ISAM, and VISIT; Table 1), the simulations from these models included in TRENDY-v9 were used in this analysis because TRENDY-v9 has more recent updates to models than MsTMIP-v2 (e.g., CLM5.0 vs. CLM4).

2.1.2 Atmospheric CO₂ observations

Atmospheric CO₂ observations are from ObsPack CO₂ GLOBALVIEWplus v3.2 (Cooperative Global Atmospheric Data Integration Project, 2017; Masarie et al., 2014) wherein continuous in situ observations were averaged to three-hourly averaged CO₂ measurements. Averaging is centered at 15:00 local time for most sites and 16:00 or 17:00 for a few sites; these times were chosen because afternoon observations are expected to have lower transport model errors stemming from model representation of planetary boundary layer dynamics (Lin et al., 2017). Urban sites were excluded as this analysis focuses on the biospheric signal. From the 44 continuous-monitoring towers for the 2007–2010 period selected here, there are around 57 700 available mid-afternoon observations, among which around 39 300 observations were used in the analysis and around 18 400 (32 %) observations were filtered out (Table S1). Data with extreme outliers and CO₂ enhancements above 30 ppmv are filtered out as described in Fang and Michalak (2015). In addition, data that have a large sensitivity to ocean fluxes (Gourdji et al., 2012) and data with potential transport model errors are also filtered out (Gourdji et al., 2012; Shiga et al., 2018a).

To isolate the biospheric enhancement or drawdown, background CO₂ values and signals from fossil fuel emissions were pre-subtracted. The impact of fossil fuel emissions on available CO₂ observations was calculated using footprints from a Lagrangian atmospheric transport model (see Sect. 2.3) and emissions from the Fossil Fuel Data Assimilation System (FFDAS v2; Asefi-Najafabady et al., 2014), scaled to $\mathrm{1}{}^{\circ }\times \mathrm{1}{}^{\circ }$ spatial resolution and three-hourly temporal resolution to be consistent with the setup of atmospheric transport. The background CO₂ values (or boundary conditions) were calculated similarly to Jeong et al. (2013), where vertical profiles from aircraft data and marine boundary layer data are used to run back trajectories and the endpoints of the back trajectories are sampled to obtain background CO₂ values. Overall, data processing and filtering were done in a similar manner to Shiga et al. (2018a) and Sun et al. (2021).

2.1.3 Absorbed photosynthetically active radiation

We use absorbed photosynthetically active radiation (APAR) as a baseline for assessing model performance. APAR is a first-order driver of gross primary productivity (GPP) (Monteith, 1972). We chose to use APAR as opposed to a GPP product, such as MODIS GPP, because we wanted to use the simplest data-driven baseline possible. We view APAR as a simple and more direct baseline than MODIS GPP because MODIS GPP itself is modeled using multiple parameters and is therefore itself a type of model (Running and Zhao, 2015). Because NEE is the balance between GPP and ecosystem respiration, and ecosystem respiration is highly correlated with GPP (Janssens et al., 2001; Baldocchi, 2008), we expect APAR to explain a portion of the variability in NEE. Given that remotely sensed APAR in and of itself does not incorporate biochemical processes governing gas exchange, which are commonly represented in TBMs, we would therefore expect models to outperform APAR in explaining observed variability in atmospheric CO₂ concentrations. APAR is calculated as the product of MODIS fAPAR (Myneni et al., 2002, 2015) and photosynthetic active radiation (PAR) following Sun et al. (2021). PAR is calculated by rescaling shortwave radiation from the North American Regional Reanalysis dataset (Mesinger et al., 2006) following the empirical relationship from Meek et al. (1984).

2.1.4 Flux tower NEE

We qualitatively compared FLUXNET2015 NEE with $\mathrm{1}{}^{\circ }\times \mathrm{1}{}^{\circ }$ modeled NEE estimates. FLUXNET2015 is a global data product for eddy covariance measurements of carbon, water, and energy exchange between the atmosphere and biosphere (Pastorello et al., 2021). We used the NEE data product with the variable USTAR threshold (VUT), where USTAR (i.e., friction velocity) thresholds vary yearly. We used data from five flux tower locations that are within the same $\mathrm{1}{}^{\circ }\times \mathrm{1}{}^{\circ }$ grid cell as towers from ObsPack CO₂ GLOBALVIEWplus v3.2 that are used to evaluate seasonality (see Sect. 2.3). The three flux towers located in the same $\mathrm{1}{}^{\circ }\times \mathrm{1}{}^{\circ }$ grid cell as the AME tower (Mead, Nebraska; Miles et al., 2012) are located at the University of Nebraska Agricultural Research and Development Center near Mead, Nebraska. The three sites are an irrigated continuous maize site (US-Ne1; Suyker, 2016a), an irrigated maize–soybean rotation site (US-Ne2; Suyker, 2016b), and a rainfed maize–soybean rotation site (US-Ne3; Suyker, 2016c). The two flux towers located in the same $\mathrm{1}{}^{\circ }\times \mathrm{1}{}^{\circ }$ grid cell as the LEF tower (Park Falls, Wisconsin; Andrews et al., 2014) are Park Falls (US-Pfa; Desai, 2016a) and Willow Creek (US-WCr; Desai, 2016b), Wisconsin. Given the scale mismatch between the footprint of flux tower observations (∼ 1 km²) and the resolution of the models examined here ($\mathrm{1}{}^{\circ }\times \mathrm{1}{}^{\circ }$), comparisons are only interpreted qualitatively.

2.2 Determining consistency across spatial scales

Bottom-up and top-down NEE estimates are compared to determine which modeling approach provides the more consistent estimate at various spatial scales. Consistency is quantified as the standard deviation across model estimates in the ensemble of TBMs and across the ensemble of inversions. We assess consistency at nested scales of $\mathrm{1}{}^{\circ }\times \mathrm{1}{}^{\circ }$, $\mathrm{3}{}^{\circ }\times \mathrm{3}{}^{\circ }$, $\mathrm{5}{}^{\circ }\times \mathrm{5}{}^{\circ }$, $\mathrm{7}{}^{\circ }\times \mathrm{7}{}^{\circ }$, and $\mathrm{9}{}^{\circ }\times \mathrm{9}{}^{\circ }$ for all grid cells throughout the North American domain by first calculating the area-weighted average NEE for each model within an ensemble and then calculating the ensemble standard deviation at each scale (Figs. 2, S1). We also assess consistency at the biome and continental scale. We then compare the consistency of TBMs to that of inversions at each scale to determine whether the more consistent approach is scale-dependent. We use an F test to determine whether differences in consistency between the two ensembles is statistically significant (p<0.05).

2.3 Evaluation against atmospheric observations

We assess the degree to which model-simulated NEE estimates can reproduce basic aspects of observed atmospheric CO₂ concentrations using three sets of metrics focusing on variability, magnitude, and seasonality. We compare observed atmospheric CO₂ concentrations with modeled CO₂ concentrations for all models in the TBM and in the inversion ensembles.

The sensitivity of CO₂ enhancements at available observation locations and times to upwind fluxes (ppm [µmol m⁻² s⁻¹]⁻¹) are represented using the Stochastic Time-Inverted Lagrangian Transport (STILT) model (Lin et al., 2003; Nehrkorn et al., 2010) driven by meteorological fields simulated by the Weather Research and Forecasting (WRF) model (Skamarock and Klemp, 2008) for North America and aggregated to a $\mathrm{1}{}^{\circ }\times \mathrm{1}{}^{\circ }$ spatial resolution and three-hourly temporal resolution. Footprints were generated as part of the NOAA CarbonTracker-Lagrange regional inverse modeling framework (Hu et al., 2019; https://gml.noaa.gov/ccgg/carbontracker-lagrange/, last access: 29 November 2023). The footprints used here cover the time period of 2007–2010. The WRF-STILT model has been previously used to assess TBM estimates of CO₂ fluxes (Fang et al., 2014; Fang and Michalak, 2015; Sun et al., 2021) and to quantify greenhouse gas fluxes (Gourdji et al., 2012; Jeong et al., 2013; Miller et al., 2014; Shiga et al., 2018a). Here we use the footprints to translate the space–time patterns of carbon fluxes into their impacts on atmospheric CO₂ observations in order to assess the model estimates' ability to represent specific aspects of atmospheric CO₂ variability in space and time.

We use the coefficient of determination (R²) between observed and modeled CO₂ drawdowns or enhancements as the metric for explained variability, i.e., quantifying the degree to which model estimates can reproduce observed spatiotemporal variability across all observation locations and times. We use the transported signal based on spatiotemporal variability in APAR as a lower benchmark for model performance (see Sect. 2.1.3). If a model has an R² value that is lower than APAR's R² value, then that model is removed from the ensemble when this metric is applied (Fig. S2).

Similarly, we use the root mean squared error (RMSE) between observed and modeled CO₂ drawdowns or enhancements as the metric for the magnitude of CO₂ signals from modeled fluxes. We again use APAR as a lower benchmark, but in this case, we first rescale APAR by minimizing RMSE. This step comes down to performing a linear regression between the transported APAR signal and the observed CO₂ enhancements, which also implicitly embodies the necessary unit conversion. If a model has a higher RMSE than that of the rescaled APAR signal, then that model is removed from the ensemble when this metric is applied (Fig. S3).

We use four sub-metrics to assess seasonality. Seasonality is assessed at individual towers that have CO₂ concentration observations for at least 50 % of days in the study period and for which the maximum data gap is less than 31 consecutive days. Only four towers within our study region meet these criteria (red symbols in Fig. 1): LEF (Park Falls, Wisconsin, USA), AME (Mead, Nebraska, USA), WKT (Moody, Texas, USA), ETL (East Trout Lake, Saskatchewan, Canada). The towers included in this analysis fall within different biomes and have different average seasonal cycles, allowing for assessment of agreement between modeled and observed CO₂ seasonal cycles across various landscapes. We also conducted a sensitivity analysis to test if the number of towers used for the seasonality metric has a significant impact on model selection. To do so, we relaxed the criteria for tower selection to a maximum allowable data gap of 80 consecutive days yielding four additional towers: SGP (Southern Great Plains, OK, USA), AMT (Argyle, ME, USA), OFR (Fir, OR, USA), and EGB (Egbert, ON Canada) (see Fig. S4). We calculate the monthly average observed and modeled CO₂ seasonal cycle for each of these four towers across the four years (2007–2010).

Figure 1 Map of biomes in North America with the locations of continuous-monitoring towers used in this study. Symbols represent the locations of towers in the CO₂ observational network. Red filled symbols represent the subset of towers with high temporal coverage that were used to evaluate how well model-simulated NEE estimates reproduce the seasonality of atmospheric CO₂ observations, whereas all towers are used for the magnitude and variability metrics (see Sect. 2.3). Red triangles represent locations of towers with high temporal coverage where there are also eddy covariance flux towers nearby (see Fig. 2).

Figure 2 Example of consistency in atmospheric inversion and TBM ensembles across spatial scales centered at two grid cells located in the cropland (CRP) and deciduous broadleaf and mixed forest (DBMF) biomes. The AME tower, along with the US-Ne1, US-Ne2, and US-Ne3 flux tower sites, falls within the CRP biome. The LEF tower, along with the US-PFa and US-WCr flux towers, falls within the DBMF biome.

The first two seasonality sub-metrics are the R² and RMSE between the monthly averaged seasonal cycles of observed CO₂ drawdowns or enhancements and of CO₂ drawdowns or enhancements resulting from the transported carbon fluxes at each of the four tower locations, for each model. We again use transported signals resulting from spatiotemporal patterns of APAR as a lower benchmark for model performance. A model with an R² value greater than that of APAR and an RMSE value less than the RMSE value for rescaled APAR is considered to meet the sub-metrics of seasonal variability and magnitude, respectively.

The third sub-metric is the amplitude of the seasonal cycle, which is defined as the difference between the maximum and minimum monthly averaged CO₂ concentrations in the average seasonal cycle (Zhao et al., 2016). The model-estimated amplitude of the seasonal cycle is evaluated at each of the four tower locations. We use the amplitude of the seasonal cycle from rescaled APAR as a lower baseline; because the seasonal cycle of APAR is less peaked than that of NEE, the same will be true for the seasonal cycle of the transported signal based on APAR relative to observed CO₂ enhancements (Figs. S5, S6). If the amplitude estimated from a model is greater than the amplitude of the transported and rescaled APAR signal, then that model is considered to meet the minimum threshold for the amplitude sub-metric.

The fourth seasonality sub-metric is the timing of peak uptake for the monthly averaged CO₂ concentrations, defined as the month when peak uptake occurs. If the predicted peak uptake falls within one month of the observed peak uptake in atmospheric CO₂, then the model is considered to pass based on this sub-metric. When applying the overall seasonality metric, a model is kept in the ensemble if it can meet the lower benchmark for at least two of the seasonality sub-metrics at all four tower locations, but the sub-metrics that it meets can vary from tower to tower.

3 Results and discussion

3.1 Model consistency across scales

We find that the full ensemble of TBMs has more consistent carbon flux estimates across all examined spatial scales relative to inversions (Fig. 3a). This is evidenced by TBMs having a smaller standard deviation across the model ensemble for the majority of locations across North America, although there are a few regions where the ensemble of inversions has a smaller standard deviation or where the ensemble with the smaller standard deviation depends on the scale being examined. This result also holds at the biome and continental scale (Fig. 5 “all models”). This result is surprising, because atmospheric inversions are informed by large-scale atmospheric constraints, while TBMs are primarily constrained by process-based understanding derived at fine scales (e.g., plot scale). One would therefore expect that inversions would be more consistent at larger scales than TBMs.

Figure 3 Maps showing whether the TBM or the atmospheric inversion ensemble has the more consistent NEE estimates across spatial scales. Maps show where each ensemble has the most consistent estimate (smallest standard deviation) at each of the following scales: $\mathrm{1}{}^{\circ }\times \mathrm{1}{}^{\circ }$, $\mathrm{3}{}^{\circ }\times \mathrm{3}{}^{\circ }$, $\mathrm{5}{}^{\circ }\times \mathrm{5}{}^{\circ }$, $\mathrm{7}{}^{\circ }\times \mathrm{7}{}^{\circ }$, and $\mathrm{9}{}^{\circ }\times \mathrm{9}{}^{\circ }$. Panel (a) shows the most consistent ensemble when the statistical significance of the difference in consistency is not taken into account. Panel (b) shows the result when statistical significance is taken into account. Green regions represent where TBMs have the smaller standard deviation at every examined scale, while blue regions show where inversions are more consistent. Orange regions represent areas where there is no statistically significant difference in consistency at any spatial scale. Yellow regions represent areas where there are inconsistencies across spatial scales. More specifically, in panel (a), yellow regions represent areas where TBMs are more consistent at some scale while inversions are more consistent at other scales, whereas in panel (b) yellow regions represent areas where there is either a mix of statistically significant and not statistically significant differences across spatial scales or where there is a statistically significant difference across all scales, but neither inversion nor TBMs are more consistent across all scales.

However, the difference in the degree of consistency between inversions and TBMs is not statistically significant (p>0.05) for large portions of the continent (Fig. 3b). This is because both ensembles have very high inter-model spread (see examples in Fig. 2), reducing the statistical significance of their differences. These results underscore the importance of statistical significance testing in interpreting model differences.

Earlier studies comparing bottom-up and top-down models have shown varying results in terms of bottom-up versus top-down model consistency for North America. Hayes et al. (2012) reported an average annual NEE estimate for North America of −931 ± 670 TgC yr⁻¹ (mean ± 1 standard deviation) across seven inverse models and −511 ± 729 TgC yr⁻¹ across 12 TBMs for the period of 2000–2006, indicating a slightly higher consistency (lower standard deviation) across inversions at the continental scale. King et al. (2015), on the other hand, reported a mean ± 1 standard deviation annual net land–atmosphere exchange for North America of −890 ± 409 TgC yr⁻¹ across 11 inverse models and −364 ± 120 TgC yr⁻¹ across 10 TBMs for the period 1990–2009, indicating that TBMs were substantially more consistent. The synthesis presented in the State of the Carbon Cycle 2 (SOCCR-2) report further confirmed that the relative consistency among top-down versus bottom-up models varies across studies (Hayes et al., 2018). These earlier studies not only used older versions of model simulations than examined here, but also included far fewer TBMs. At the global scale, the most recent assessment shows that TBMs have a greater model spread (i.e., lower consistency) than do inversions (Friedlingstein et al., 2022). Though these studies are primarily focused on assessing agreement between bottom-up and top-down methods at a single large spatial scale, they demonstrate the difficulty in assessing consistency when ensemble size is limited, and statistical significance is not evaluated. In addition, comparing consistency of bottom-up and top-down models without assessing the statistical significance of observed differences may lead to misleading conclusions. Kondo et al. (2020) found TBMs to have a smaller inter-model spread in regional budget estimates, but the large spread in the seasonality of carbon uptake for TBMs made it difficult to deduce whether bottom-up models are more reliable than top-down models based on consistency alone.

Figure 4 Maps showing where the TBM or the atmospheric inversion ensemble has the more consistent NEE estimates across spatial scales when the ensembles are limited to those models that meet (a) variability, (b) seasonality, (c) magnitude, or (d) all three metrics. Colors are as defined in Fig. 3.

Figure 5 Estimated NEE for TBMs and atmospheric inversions for all models within their respective ensembles as well as subsets of the ensembles that meet the variability, seasonality, or magnitude metrics, or that meet all three. Boxplots represent the model-specific average NEE estimates across the models included in each ensemble; triangles represent the across-model mean. Panel (a) shows NEE for North America, while panels (b)–(d) show NEE for specific biomes. Stars represent cases for which there is a statistically significant difference (p<0.05) in the consistency of the TBM versus the inversion ensemble.

While the consistency of model ensembles can be used as one measure of uncertainty in modeling carbon uptake, both bottom-up and top-down methods carry uncertainties that may not be fully captured by ensemble spread alone. For example, the use of satellite data to augment data coverage for atmospheric inversions did not clearly improve consistency in inverse-model-based estimates in a recent intercomparison study (Crowell et al., 2019). Moreover, the differences in the processes incorporated in the various bottom-up models may impact the consistency across the TBM ensemble. We examined how using a simple definition of NEE ($R_{\mathrm{h}}+R_{\mathrm{a}}-$ GPP) impacted consistency within the MsTMIP-v2 ensemble and found that estimates from the MsTMIP-v2 ensemble are less consistent across models when using only GPP, R_h, and R_a in the calculation of NEE (Fig. S7). In other words, the models are more consistent when they include other components of NEE, even though those components differ from model to model. This seems to suggest that models may implicitly target a presumed net land sink irrespective of the processes included. This appears to be in line with earlier research that found some TBMs arrive at similar estimates of carbon uptake even though they show large disagreements on the primary driver of increased uptake in recent decades, while other TBMs arrive at dissimilar estimates despite having similar sensitivities (Huntzinger et al., 2017).

3.2 Impact of variability, magnitude, and seasonality on consistency

We evaluate whether greater consistency corresponds to greater accuracy by subsetting the model ensembles using metrics based on the variability, seasonality, and magnitude of atmospheric observations (see Sect. 2.3). Limiting the bottom-up ensemble to only include TBMs that reproduce the variability of atmospheric observations better than APAR reduced the ensemble size from 29 to 9 models. In other words, over two-thirds of TBMs represented the space–time variability of atmospheric CO₂ less well than did APAR (Table 1, Fig. S2). Conversely, all inversions performed better than this benchmark (Table 2), which is expected given that the inversions use all or a subset of the same observations in estimating fluxes. Sub-selecting models that could represent aspects of the seasonality of atmospheric observations reduced the ensemble of TBMs from 29 to 11 and the ensemble of inversions from 8 to 6. Limiting the ensembles to models that could represent the magnitude of observed atmospheric CO₂ signals reduced the ensemble of TBMs from 29 to 5, while the ensemble of inversions was reduced from 8 to 6. Only four TBMs (and six inversions) remained when all three metrics were applied. For the four TBMs that are in both the MsTMIP-v2 and TRENDY-v9 ensembles, we retained the TRENDY-v9 simulations as described in Sect. 2.1.1. Had we included the three MsTMIP-v2 simulations, however, none of them met any of the three metrics and they therefore would have been excluded from the subsets meeting the variability, seasonality, and magnitude metrics.

The large reduction in the ensemble size when even basic benchmarks related to observed atmospheric CO₂ signals are applied (i.e., filtering based on a minimum accuracy threshold) indicates that ensemble spread (i.e., consistency) is unlikely to be a good indicator of actual uncertainty in our understanding of North American carbon balance. An example of consistency not necessarily capturing uncertainty is shown in Fig. 2 at the $\mathrm{1}{}^{\circ }\times \mathrm{1}{}^{\circ }$ scale where TBMs and inversions show better agreement with eddy covariance flux tower observations in the deciduous broadleaf and mixed forest biome (Fig. 2b) than in the cropland biome (Fig. 2a) despite having similar model consistency in both biomes. It is unclear whether the models and observations disagree in the cropland biome due to sub-grid scale heterogeneity (Melton and Arora, 2014) versus inaccuracies in the models (Schuh et al., 2014; Guanter et al., 2014; Sun et al., 2021), but this comparison nevertheless illustrates how consistency may not capture the full extent of uncertainty in model simulations.

Fang et al. (2014) also noted that atmospheric observations can be used to evaluate flux patterns from TBMs in terms of models' ability to explain atmospheric observations. Using a similar approach here we see that differentiating between models that reproduce basic features of atmospheric CO₂ signals leads to a shift from TBMs having the smaller ensemble standard deviation across scales to there being no statistically significant difference between the ensemble standard deviations of TBMs and inversions across most of the continent. This is particularly true when models are selected based on consistency with all three metrics.

Applying the variability, magnitude, and seasonality metrics reduces the areas in North America for which TBMs have a statistically significantly greater consistency than do inversions (Fig. 4). Once all three metrics are applied, very few areas remain where TBMs have a higher consistency across all scales. This result is enlightening, because it indicates that once basic aspects of atmospheric observational constraints are taken into account, apparent differences in consistency between approaches disappear. The large reduction in the TBM ensemble size is part of the reason for this change, so this result must be interpreted with caution. In other words, this result is less a product of the consistency of remaining inversions increasing and more a product of the statistical significance of differences in consistency being reduced due to the reduction in the sample size. This filtering exercise also indicates that, although the ensemble of inversions examined here has lower consistency, it may actually exhibit higher accuracy as evidenced by the smaller reduction in ensemble size and higher number of models that meet all three criteria.

Table 2 Impact of filtering ensembles based on key metrics derived from observational constraints. The mean ± standard deviation (median) fluxes (in units of TgC yr⁻¹) for inversions and TBMs when metrics are applied to North America and biomes with the largest data constraint across 2007–2010.

3.3 Implications for understanding of North American carbon balance

3.3.1 Impact on model consistency

The impact of applying accuracy metrics to the ensembles offers a glimpse into the true North American carbon sink. While this approach did increase model consistency for TBMs, it did not impact model consistency for atmospheric inversions. The lack of increase in the consistency of inversions once accuracy metrics are applied indicates that there is a wide range of overall flux patterns and magnitudes that are consistent with large-scale atmospheric constraints. It is interesting, therefore, that model consistency among TBMs does increase substantially when accuracy metrics are applied (Fig. 5), although again this has to be interpreted with caution given the small number of remaining models. This contrast implies that, once accuracy metrics are applied, remaining TBMs all reproduce observed features of atmospheric observations using similar flux patterns, while remaining inversions reproduce observed features of atmospheric observations using a wider range of fluxes. The large spread of inversions that are also consistent with the same atmospheric constraints indicates that consistency across model ensembles is likely a poor indicator of overall uncertainty. This suggests that decreases in model spread do not necessarily indicate increases in model accuracy, as has also been suggested in previous studies (Annan and Hargreaves, 2010; Knutti et al., 2010; Kondo et al., 2020; Lovenduski and Bonan, 2017). Because TBMs were already informed by process-based understanding, subsetting the ensemble to those TBMs that are also consistent with broad features of atmospheric CO₂ observations may lead to the “best of both worlds”. The high consistency of the remaining (albeit small number of) TBMs is a sign that these models are more likely to capture the true North American carbon sink.

3.3.2 Impact on model agreement

Most promisingly, not only does the difference in consistency between top-down and bottom-up models decrease when ensembles are filtered based on key metrics derived from observational constraints, but the agreement between top-down and bottom-up models improves dramatically as well (Fig. 5a). Indeed, once all three metrics are applied, both the mean and the median across the remaining TBMs are very close to the mean and the median across the remaining inversions. In other words, filtering leads to better agreement between top-down and bottom-up estimates of North American carbon balance, a goal that has proven elusive up to now. The mean (median) North American fluxes across 2007–2010, when all three criteria are applied, is −465 TgC yr⁻¹ (−499 TgC yr⁻¹) for the TBMs and −597 TgC yr⁻¹ (−574 TgC yr⁻¹) for the inversions (Table 2).

Agreement for the three biomes that are best-constrained by atmospheric observations is also improved (Fig. 5b, c, d) although to a lesser extent. These biomes are evergreen needleleaf forests, croplands, and deciduous broadleaf and mixed forests biomes, based on an analysis by Sun et al. (2023) that showed these biomes to have the greatest sensitivity of atmospheric observations to fluxes (2023). When the seasonality metric is recalculated using eight towers selected from the sensitivity analysis (see Sect. 2.3), there are minimal changes to agreement between bottom-up and top-down estimates. The deciduous broadleaf and mixed forests, which is the biome with the lowest observational constraint of the three biomes, showed decreased agreement suggesting that using towers with high temporal data availability, rather than more towers, is important for capturing seasonality. The mean (median) uptake for evergreen needleleaf forests when all three metrics are applied is −118 TgC yr⁻¹ (−118 TgC yr⁻¹) for the TBMs and −146 TgC yr⁻¹ (−132 TgC yr⁻¹) for the inversions (Table 2). For croplands, the mean (median) is −111 TgC yr⁻¹ (−82 TgC yr⁻¹) for the TBMs and −164 TgC yr⁻¹ (−135 TgC yr⁻¹) for the inversions (Table 2). The mean (median) for deciduous broadleaf and mixed forests is −88 TgC yr⁻¹ (−80 TgC yr⁻¹) for the TBMs and −178 TgC yr⁻¹ (−194 TgC yr⁻¹) for the inversions (Table 2).

Lower improvement in the agreement between TBMs and inversions at the biome scale relative to the continental scale may have resulted from disagreement on where major sinks in North America lie. However, to understand why models disagree, it is necessary to understand what gives rise to differences. For example, deciduous broadleaf and mixed forests were found to account for the majority of the interannual variability in NEE for North America when using a top-down approach, but TBMs disagreed on whether forested or non-forested biomes contribute most strongly to interannual variability and what the primary environmental drivers of this variability are (Shiga et al., 2018b). TBMs have also been shown to have greater interannual variability in western temperate North America than in eastern temperate North America, albeit with substantial model spread (Byrne et al., 2020). A recent study found that TBMs that did well at reproducing observed CO₂ variability exhibited substantially stronger growing-season carbon uptake in croplands relative to other models (Sun et al., 2021). These studies highlight that there is still uncertainty about the geographic distribution, interannual variability, and climatic drivers of North American carbon uptake. This likely plays a role in the reduced agreement at the biome scale relative to the continental scale observed here.

3.3.3 Sensitivity analyses and study limitations

One potential reason for the lack of convergence between bottom-up and top-down approaches is the limited number of towers used for the seasonality metrics. These towers are primarily located in the midcontinent (Fig. 1), although observations from towers are also influenced by fluxes in other regions (see, e.g., Extended Data Fig. 2 in Sun et al., 2023). To test this this hypothesis, we conducted a sensitivity analysis to explore the impact of including four additional towers (see Sect. 2.3). We found that including more towers in other parts of North America did not change our primary findings. Specifically, we found that 11 models perform well based on the seasonality metric irrespective of which subset of towers is used (i.e., original four towers, additional four towers, or all eight towers), while six and four additional models also meet the metric when the original and additional sets of four towers are used, respectively. Though modifying the towers used impacts the consistency of ensembles based on the seasonality metric alone, the impact on the consistency and agreement of models that meet all three metrics is minimal (Fig. S8).

While increased data availability in biomes that are not well sampled would allow for a more robust evaluation of models' ability to capture key aspects of seasonality, several studies have found that even in areas with high data availability there is still noticeable disagreement between inversions (Bastos et al., 2020; Kondo et al., 2020). This indicates that model-specific issues may be more important than data availability in explaining discrepancies between bottom-up and top-down estimates. In the context of inverse modeling, the choice of fossil fuel inventory and transport model have been identified as significant sources of uncertainties, though there is not a clear consensus as to what the primary source of uncertainty is (Bastos et al., 2020; Gaubert et al., 2019; Peylin et al., 2011; Schuh et al., 2019). For TBMs, noteworthy uncertainties stem from how (and whether) models represent key processes such as land use change, wood and crop harvesting, fire, and vegetation dynamics (Ahlström et al., 2015; Bastos et al., 2020; Hardouin et al., 2022; Tharammal et al., 2019).

Several studies also highlight lateral fluxes as a significant source of CO₂ emissions, yet many TBMs do not account for them (Drake et al., 2018; Kondo et al., 2020; Raymond et al., 2013). To test the possible impact of lateral fluxes on our analysis, we examined the effect of including lateral fluxes in bottom-up estimates on the agreement between bottom-up and top-down methods. To do so, we added gridded lateral fluxes of river export, crop trade, and wood trade to TBMs to make them more comparable with fluxes seen by inversions. We used two different estimates of river export from Byrne et al. (2023). The first is the gridded product, which incorporates results from the Global NEWS model (Byrne et al., 2022). The second is the same gridded product rescaled, so total river export equals the country total river exports reported in Byrne et al. (2023), which is the mean of two models estimates (Global NEWS and DLEM). The differences in these two river export estimates highlight some of the uncertainties associated with estimating lateral fluxes (Byrne et al., 2023; Drake et al., 2018). We find that at the North American scale, incorporating lateral fluxes improved agreement between bottom-up and top-down models somewhat, but the change was not sufficient to explain discrepancies for the best-constrained biomes (Fig. S9). In the deciduous broadleaf and mixed forests and cropland biomes, lateral fluxes only partially explain discrepancies and applying the seasonality, variability, and magnitude metrics still improved agreement between bottom-up and top-down estimates (Fig. S9c–d). This contrasts with the evergreen and needle leaf forest biome where the inclusion of lateral fluxes led to better agreement between inversions and TBMs for the subset with all models included but ultimately exacerbated differences once models that meet all three criteria were selected (Fig. S9b). This aligns with the findings of Kondo et al. (2020) that accounting for lateral fluxes and using a standardized net CO₂ flux definition was not sufficient to explain discrepancies between inversions and TBMs at regional scales.

Next we explored whether agreement between top-down and bottom-up estimates was improved when a more consistent definition of NEE was applied to TBMs. To do so, we examined how using a simple definition of NEE ($R_{\mathrm{h}}+R_{\mathrm{a}}-$ GPP) impacted agreement between MsTMIP-v2 models and inversions; we found that agreement with inversions improved slightly (Fig. S7). Of the MsTMIP-v2 models we looked at, only four included both disturbance and product fluxes in their definition of NEE (CLM4, CLM4VIC, TEM6, VEGAS2.1), and how well these models agree with inversions varies by biome (Fig. S7). Bastos et al. (2020) found that models, irrespective of their inclusion of fire dynamics, exhibited similar performance, which could be attributed to heightened sensitivity of decomposition to temperature in models without fire. Moreover, Huntzinger et al. (2017) found TBMs can yield similar estimates despite diverging on their primary drivers. We therefore found that it is unlikely that varying NEE definitions provide a comprehensive explanation for the observed disparities.

There are additional sources of uncertainty that the sensitivity analyses described here cannot address. For example, the non-uniform spatial distribution of observation data throughout the domain limits the biomes that can be examined. Additional observations, particularly in regions that are not well sampled, would allow for additional analyses and would further increase the robustness of the continental-scale results. Leveraging larger ensembles of TBMs and inversions and incorporating an assessment of the impact of transport model choice would also allow for additional insight into the robustness of our findings. Though we showed that TBM estimates of North American carbon sink were actually more consistent when varying NEE definitions were used and that discrepancies between bottom-up and top-down models cannot be fully resolved by the incorporation of lateral fluxes, improved standardization across TBMs in terms of represented processes would allow for a more accurate comparison with inversions. Finally, it is important to recognize that better performance under present climate conditions does not necessarily translate to better model performance under future conditions.

4 Conclusions

Comparing estimates from bottom-up and top-down methods across spatial scales and evaluating estimates in light of atmospheric CO₂ observations is useful for exploring persistent differences between these approaches. We show that the difference in consistency between bottom-up and top-down ensembles is not statistically significant for large regions of North America because of large variability within both ensembles, highlighting the importance of significance testing in interpreting model differences. We also find that ensemble spread is unlikely to be a good indicator of overall uncertainty in the North American carbon balance. This is because when the same benchmarks based on observed atmospheric CO₂ signals are applied to both ensembles, inversions use a wider range of fluxes than TBMs to reproduce observed atmospheric observation features. Encouragingly, once models are sub-selected based on their ability to reproduce basic aspects of observed atmospheric CO₂ variability, seasonality, and magnitude, bottom-up and top-down estimates of North American carbon balance agree at the continental scale and for large biomes therein. Notably, these findings remained robust after several sensitivity analyses were performed. The convergence in flux estimates between top-down and bottom-up approaches demonstrates the usefulness of filtering models based on their agreement with even basic features of large-scale observational constraints for assessing our understanding of carbon budgets. This finding is encouraging because it presents a promising path towards both improving model consistency and reducing uncertainties. Thus, continued efforts to reduce uncertainties should focus on improving consistency at scales finer than large continental domains and leveraging top-down observational constraints to refine understanding of the North American carbon balance.