Bias-correcting carbon ﬂuxes derived from land-surface satellite data for retrospective and near real-time assimilation systems

. The ability to monitor and understand natural and anthropogenic variability in atmospheric carbon dioxide ( CO 2 ) is a growing need of many stakeholders across the world. Systems that assimilate satellite observations, given their short latency and dense spatial coverage, into high-resolution global models are valuable, if not essential, tools for addressing this need. A notable drawback of modern assimilation systems is the long latency of many vital input datasets, e.g., inventories, in situ measurements, and reprocessed remote-sensing data can trail the current date by months to years. This paper describes 5 techniques for bias-correcting surface ﬂuxes derived from satellite observations of the Earth’s surface to be consistent with constraints from inventories and in situ CO 2 datasets. The techniques are applicable in both short-term forecasts and retrospective simulations, thus taking advantage of the coverage and short latency of satellite data while reproducing the major features of long-term inventory and in situ records. Our approach begins with a standard collection of diagnostic ﬂuxes which incorporate a variety of remote-sensing driver data, viz. vegetation indices, ﬁre radiative power, and nighttime lights. We then apply an 10 empirical sink so that global budgets of the diagnostic ﬂuxes match given atmospheric and oceanic growth rates for each year. This step removes coherent, systematic ﬂux errors that produce biases in CO 2 which mask the signals an assimilation system hopes to capture. Depending on the simulation mode, the empirical sink uses different choices of atmospheric growth rates: estimates based on observations in retrospective mode and projections based on seasonal forecasts of sea surface temperature in forecasting mode. The retrospective ﬂuxes, when used in simulations with NASA’s Goddard Earth Observing System (GEOS), reproduce marine boundary layer measurements with comparable skill to those using ﬂuxes from a modern inversion system. The forecasted ﬂuxes show promising accuracy in their application to the analysis of changes in the carbon cycle as they occur. to and oceanic growth rates. We the ﬂux monthly ﬁner. empirical sink, we ﬂux to et al.

model simulations before assimilating satellite data, thus increasing the ability of the system to capture other signals of interest in the observations (Dee, 2005). For example, in any given year, estimates of global net biospheric exchange (NBE) from the TRENDY ensemble of terrestrial biosphere models (Sitch et al., 2015) have a range of roughly 4 petagrams of carbon (Pg C). 60 Transported through the atmosphere, this produces a 2 ppm spread in global CO 2 , which is twice the size of the atmospheric growth rate perturbation of the 2015-2016 El Niño. An under/over-estimation of NBE can lead to an under/over-prediction of the seasonal cycle amplitude of simulated CO 2 (Yang et al., 2007;Keppel-Aleks et al., 2012;Zhao et al., 2016), and has the potential for confusion with other, covarying error sources (Basu et al., 2011). Among flux inversion systems that ingest in situ data, e.g., those analyzed by Gaubert et al. (2019), the spread in NBE is usually much smaller, about 0.25 Pg C. Nevertheless, 65 flux inversions typically achieve this agreement through their reliance on long-term data records and transport simulations, and can trail the current date by a year or more. Our goal is to achieve a similar error reduction without the resulting latency, allowing the subsequent assimilation of satellite CO 2 data to focus on regional and seasonal errors rather than global budgeting errors.
This paper presents a collection of surface fluxes having both retrospective and forecasting modes that reproduce background 70 CO 2 measurements with comparable skill to a modern flux inversion system. The collection includes an additional, empirically derived land sink to our baseline flux collection to ensure that global flux totals are consistent with observed atmospheric CO 2 growth rates. Empirical adjustments of this kind date back at least to the work of Tans et al. (1990), who showed that their simulations recreated observed North-South gradients of CO 2 when they increased terrestrial uptake in the Northern an atmospheric growth rate based on in situ observations in the marine boundary layer (MBL; Dlugokencky and Tans, 2016b). 80 In forecasting mode, when many of the diagnostic flux products and observationally constrained growth rates are unavailable, the fluxes use extrapolation and a method for predicting the growth rate based on sea surface temperature forecasts (Jones and Cox, 2005;Betts et al., 2016). It differs from previous works in its functional form and its application to fluxes diagnosed from satellite measurements instead of those from a prognostic model. For now, we are willing to accept that the empirical sink adds a biophysical inconsistency into our system but hope to better address this in the future. 85 In the work that follows, the construction of the empirical sink is detailed in Section 2. The evaluation of the fluxes against other products and comparisons of transported mixing ratios to in situ measurements are presented in Section 3. In particular, we show that a GEOS simulation using the bias-corrected fluxes has similar skill reproducing MBL measurements to a simulation using fluxes from a modern inversion system. These findings and their potential scientific impact are summarized in Section 4.

90
Here we present the Low-order Flux Inversion (LoFI), a collection of carbon fluxes driven by remote-sensing land-surface data and bias-corrected to reproduce given atmospheric and oceanic growth rates. We use the term "low-order" to distinguish it from modern flux inversion systems which typically solve for fluxes at a regional, monthly scales or finer. To demonstrate the impact of the empirical sink, we also consider a "baseline" flux collection similar to that used in Ott et al. (2015) with the only difference being that the empirical land sink is removed. All products in both the LoFI and baseline collections are 95 conservatively regridded to the 0.5 • × 0.625 • dateline and pole centered grid used by GEOS for many applications including the MERRA-2 reanalysis of the meteorological inputs needed by the flux components (Section A1). All fluxes have a daily timestep, except the land biosphere which has a 3-hourly timestep to resolve the diurnal cycle. In retrospective mode, the LoFI collection trails the current date by up to a year. It forecasts fluxes for the current and next year using the modifications described in Section 2.3. As a rough metric of truth/plausibility, we compare both collections to widely-used ensembles of 100 terrestrial biospheric, ocean biogeochemical models, and flux inversions. These choices and the components of the LoFI fluxes are detailed below and summarized in Table 1.

Individual flux components
The LoFI flux collection consists of the following six components: Net Ecosystem Exchange (NEE) -An implementation of the Carnegie Ames Stanford Approach (CASA; Randerson et 105 al., 1996) and Global Fire Emissions Dataset version 3 (GFED 3; van der Werf et al., 2003;2010) referred to here as CASA-GFED 3. CASA-GFED 3 uses satellite-based measurements of land cover and vegetation changes along with meteorology from MERRA-2 to constrain carbon stocks and fluxes, viz. net primary productivity (NPP), which is determined from measurements of normalized difference vegetation index (NDVI; Pinzon and Tucker, 2014), and biomass burning, which is determined from Moderate Resolution Imaging Spectrometer (MODIS) burned area estimates Fossil fuel combustion -The Open-source Data Inventory for Anthropogenic CO 2 (ODIAC; Maksyutov, 2011, 2015;Oda et al., 2018). ODIAC is a global, monthly, high-resolution (1 km × 1 km) fossil fuel CO 2 gridded emission data product based on the disaggregation of country-level fossil fuel CO 2 emission estimates using a global power plant This approach reapplies the global mean climatological growth rate estimate of 1.5 µatm/yr for the partial pressure 130 of CO 2 in seawater (pCO sw 2 ) that Takahasi et al. (2009) use to derive their climatology. For the partial pressure in the atmosphere (pCO atm 2 ), we use weekly values of zonal-mean surface CO 2 from the NOAA MBL reference (Masarie and Tans, 1995;Dlugokencky and Tans, 2016a). Given the two partial pressures, estimates of the flux from the ocean surface to the atmosphere follow from the expression (Wanninkhof, 2014) F sw = kCU 2 10 (pCO sw 2 − pCO atm 2 ),

135
where k is a constant, U 10 is the 10-meter wind speed, and C is the fractional sea-ice coverage. To complete the flux calculation, we use daily, observationally-constrained estimates of U 10 and C from MERRA-2. This approach has been used by several previous studies and is derived from one of the ocean priors in the NOAA CarbonTracker System ( where ∆ + T m is the (non-dimensionalized) temperature increase from the previous month, and α is a constant scaling factor computed such that the global total fluxes for the year match a specified atmospheric growth rate ∆CO 2 . In other 150 words, where 2.124 · ∆CO 2 is the atmospheric growth rate in units of Pg C (Ballantyne et al., 2012), X is the area-weighted global, annual total in units of Pg C of a flux field variable X, and F is the sum of all baseline fluxes. In retrospective years (those preceding the current), we use growth rates derived from the NOAA MBL reference (Dlugokencky and 155 Tans, 2016b), and in NRT years (the current and following) we use projections based on seasonal forecasts of sea surface temperature described in Section 2.3.

Anthropogenic short-cycle burning and lateral fluxes
The above separation into component fluxes assumes that, added together, biomass burning and biofuel emissions account for emissions from both naturally occurring wildfires and anthropogenic burning due to gross land-use/land-cover change. Fur-160 thermore, emissions from ethanol, biodiesel, and other short-cycle fuels used for transportation are potentially underestimated since they are not included in the ODIAC fossil fuel and CASA-GFED 3 biofuel emissions, yet the removal of carbon due to the corn and soybean harvest in the Midwestern United States is included in CASA-GFED 3 (derived from USDA National Agricultural Statistics Service data for 2005). Excluding the lateral transport of significant amounts of carbon over continental scales, we do not expect these assumptions to have a noticeable effect on the comparisons to follow since they consider only 165 NBE, the sum of all land fluxes except fossil fuel emissions, and not the individual components.

Modifications needed for forecasting mode
Many of the products used in the above fluxes are unavailable until a few months to years following the end of a given year.
In particular, fossil fuel inventories and the NOAA MBL growth rate require data collection and analysis that must necessarily The primary motivation for using QFED biomass burning instead of GFED was QFED's ability to produce NRT estimates.

175
This ability is particularly important for biomass burning emissions whose inter-annual variability is comparable in magnitude to its annual mean and seasonal cycle amplitude and is much greater than its long-term trend. While the removal of terrestrial carbon by QFED fires does not match the corresponding fire loss in CASA carbon stocks, QFED emissions are calibrated to GFED (Darmenov and da Silva, 2015), resulting in a difference that is minor compared to that of the empirical sink (see the Supplementary Section A3).

180
For the empirical sink, we switch from using an atmospheric growth rate derived from in situ data in retrospective years to a linear functional fit to seasonal SST anomalies and anthropogenic emissions (Jones and Cox, 2005;Betts et al., 2016). This approach estimates the growth rate of CO 2 in ppm, ∆CO 2 as where N is the average of SST anomalies in the Niño 3.4 region (5 • N to 5 • S, 170 • W to 120 • W) from October of the previous 185 year to September of the current year in units of degrees Kelvin, and E is the global total anthropogenic emissions from fossil fuels and net land-use/land-cover change in units of Pg C. For the SST anomalies, we use the Reynolds analysis (Reynolds et al., 2007) until the last month it is available and fill in the remaining months with SST forecasts from the GEOS Subseasonal to Seasonal (S2S) forecast system (Molod et al., 2020), and for the anthropogenic emissions we use extrapolated totals from the Global Carbon Project (GCP) 2018 budget.

190
The coefficient 0.205 multiplying the anthropogenic term in the growth rate forecast (Equation 2) represents a constant airborne fraction of 0.44 when converted using a factor of 2.124 Pg/ppm, roughly in line with the findings of Raupach et al. (2014). Historical data records tend to show that the long-term trend in the airborne fraction is insignificant compared to its inter-annual variability (Knorr 2009;Ballantyne et al., 2012), but there is some indication of significant multi-decadal trends, including the possibility of a recent decrease (Keenan et al., 2016).

Flux and transport simulation analysis
Our primary method for evaluating the LoFI flux package, along with the baseline fluxes, are comparisons to the ensembles of bottom-up terrestrial biosphere and ocean biogeochemistry models and top-down flux inversions outlined in Table 1 , 1995;Ludwig et al., 1996) to estimate a riverine contribution of 0.45±0.18 Pg C (all uncertainties reported here are 1σ). Recently, Resplandy et al. (2018) showed that relationships between ocean heat and carbon transport to derive an estimate of 0.78 ± 0.20 Pg C. Since there is so much uncertainty about this discrepancy, in particular its distribution in time and space, we do not make any corrections There are several other flux evaluation metrics, which we do not apply here, each with their own limitations. For example, fluxes can be measured directly from towers with eddy covariance techniques (Dabberdt et al., 1993), but the spatial footprint 210 of a flux tower is typically much smaller than ten kilometers (Raczka et al., 2013), and upscaling flux tower data to a global, gridded product has thus far been unable to produce reliable global NBE budgets (Jung et al., 2011). Eddy covariance measurements from aircraft (Desjardins et al., 1989), e.g., NASA's Carbon Airborne Flux Experiment (CARAFE; Wolfe et al., 2018), are representative of much longer horizontal length scales than towers, but the currently available campaign data are sparse in space and time. As shown in Section 3.3, another choice is to transport our fluxes through the atmosphere with the GEOS 215 general circulation model (GCM) and compare the simulated CO 2 mixing ratios with atmospheric observations. While this approach can test the ability of the fluxes to reproduce large-scale, long-term signals, it is susceptible to the misinterpretation of a transport error as a flux error.

Ocean exchange
Over the past few decades, global pCO sw 2 has increased by roughly 1.5 µatm/yr, yet with considerable regional differences By imposing observed trends in pCO sw 2 and pCO atm 2 , our ocean exchange fluxes produce a sink that is generally consistent 225 with the inversion ensemble and GCP 2018 ocean biogeochemical model ensemble (Figures 1 and A1), again with the provision that the biogeochemical model ensemble does not include outgassing due to riverine input. Averaged over 2000-2010, our mean ocean sink is 1.55 Pg C. This is consistent with annual budgets based on atmospheric measurements of a combination of of O 2 and CO 2 called atmospheric potential oxygen (APO; Stephens et al., 1998): Keeling and Manning (2014) estimate an ocean sink of 2.50 ± 0.60 Pg C for 2000-2010, which reduces to 1.72 ± 0.60 Pg C after removing 0.78 Pg C to account for riverine 230 input (see discussion above). The pCO sw 2 based products of Landschützer et al. (2016) and Rödenbeck et al. (2014), which do not require a riverine adjustment, produce sinks of 1.37 and 1.74 Pg C (averages from GCP 2018; Le Quéré et al., 2018), also in line with our budgets (several more comparable estimates are available in Jacobson et al., 2007). Even during the strong ENSO of 2015-2016, the constant growth in pCO sw 2 that we impose produces a global sink within the ranges of the top-down and bottom-up ensembles (Figure 1). While linear growth is not appropriate for simulations of several decades, over which 235 pCO sw 2 increases exponentially (Raupach et al., 2014), it is sufficient for our 15-year study period and something we hope to address in future developments. and biomass burning-by over 3 Pg C (see Figure 2a and Table 2). This is due to the assumption of a neutral biosphere in CASA-GFED, which thus has no long-term net sink, and is the primary motivation for the derivation of the empirical land sink. In the Northern Hemisphere, the neutral biosphere assumption causes CASA to systematically underpredict seasonal cycle amplitudes in comparison to measurements from flux towers (Carvalhais et al., 2008) and a combination of aircraft and ground-based remote-sensing retrievals (Yang et al., 2007). This effect is seen in the diagnostic fluxes as an underprediction of 245 the global sink strength in March through July (Figure 2) limited primarily to the Northern Extratropics (Figure 3, first column).
In contrast, including the empirical sink in the LoFI fluxes moves its annual totals and seasonal cycle significantly closer to the ranges of the comparison ensembles. Since the empirical sink increases the seasonal cycle and its magnitude increases as net terrestrial uptake grows, it produces a seasonal cycle amplitude that increases in time, consistent with observations of in high northern latitudes (Graven et al., 2013) and their attribution to increased extratropical terrestrial uptake (Barnes et al., 2016).  Table 2 for comparable ranges from the TRENDY and inversion ensembles).  The net terrestrial sink found by inversions has been shown to be consistent with growth in temperate and boreal terrestrial ecosystems driven primarily by carbon fertilization and forest regrowth (Schimel et al., 2015;Fernández-Martinez et al., 2019) with the possibility of an additional contribution from agriculture (Zeng et al., 2014). Notably, the forest inventory analysis Taken together, the findings described above suggest an empirical sink proportional to the CO 2 growth rate (e.g., driven by 275 carbon fertilization) and monthly temperature increase (e.g., focused to the extratropical growing season). The choice of the last remaining factor in Equation 1, heterotrophic respiration (HR) over, for example, Net Primary Production (NPP), requires further investigation. We chose HR rather than NPP for several reasons: 1) boreal forests allocate a much greater percentage of biomass below ground than tropical forests do (Pan et al., 2011), 2) the NDVI driver data of CASA tends to put a strong, observationally-driven constraint on NPP, and 3) the uncertainty in the HR response to temperature changes through its so-280 called Q 10 function (Lloyd and Taylor, 1994;Tjoelker et al., 2001;Huntzinger et al., 2020). In any case, the temperature increase factor of the empirical sink makes it relatively insensitive to the choice of the flux factor (e.g., HR, NPP, or ecosystem respiration).
Even after adding the empirical sink, the LoFI fluxes have some discrepancies with the comparison ensembles worth noting.
In particular, the LoFI fluxes predict a sink from 0 • S to 15 • S during JJA, while the comparison ensembles predict neutral 285 fluxes ( Figure 3, first column, fourth row difference between blue line and grey shading). The opposite difference, although less noticeable, is present during DJF (Figure 3, first, column, second row). It is unclear, however, how accurate either ensemble is in this case: the inversions are hindered by the limitation of in situ data and uncertainties in transport over this latitude band, while the biospheric models are potentially hindered by inaccuracies in meteorological driver data and deficiencies in the representation of disturbance (Molina et al., 2015;van der Laan-Luijkx et al., 2015). This approach may overestimate the net 290 sink over the Midwestern United States, where our version of CASA-GFED 3 already includes a corn and soybean harvest (other harvests are not currently represented). Understanding the interaction of these two adjustments is the subject of ongoing work.

Transport simulations
Ultimately, a major goal of this effort is to develop a realistic collection of fluxes that improves the ability of transport sim- constrained by analysis fields every six hours, while physical processes such as convection, turbulence, and radiative transfer are recalculated at a high temporal resolution. This computationally efficient framework provides the ability to simulate realistic meteorology with a tight coupling between fine scale atmospheric transport processes and trace gas emissions. For the analyzed meteorology, we use MERRA-2, resulting in a transport simulation sharing many of the properties of that used by Ott et al. (2015).

310
The evaluation of the two model runs (one with LoFI fluxes, the other with CT2016 fluxes) against the NOAA MBL surface sites is shown in Figure 4. Neither of the runs is clearly superior. This suggests that, at least in terms of aggregate statistics over multiple years, the empirical sink produces a correction to the baseline, diagnostic fluxes with a similar skill as running a formal inversion system based on MBL data. While our approach may degrade the skill of CT2016 by using a different transport model, it may also improve it by running at a higher resolution, and our calculated differences are consistent with

Growth rate forecasts
At the start of 2016, the LoFI fluxes switch from retrospective mode to NRT. While it would be possible to extend retrospective magnitude. This suggests that misrepresentations of the spatial and temporal variability of El Niño may have a greater impact on our ability to represent the Mauna Loa data record than errors in the growth rate forecast.
A major factor in the ability of NRT LoFI to reproduce atmospheric CO 2 observations is that, even in forecasting mode, it continues to use the retrospective meteorological analysis from MERRA-2. If we were to use meteorological forecasts instead, 330 we would expect the skill to degrade within a few months and last about two years or less (Ilyina et al., 2020).

Conclusions
This paper presented an adjustment to a diagnostic collection of surface fluxes designed to bias-correct its global budgets to match inventory data and in situ observations. For the adjustment, we developed an empirical sink that was the product of the monthly increase in temperature with heterotrophic respiration, thus focusing the correction to the Northern Extratropics  Barnes et al., 2016).
This suggests that diagnostic models of the terrestrial biosphere have the potential to reproduce these patterns with the addition of a transient carbon pool representing the net sink in the NE.

340
Using the fluxes bias-corrected with the empirical sink in transport simulations reproduced atmospheric measurements in the marine boundary layer with the same skill as transport simulations using fluxes from a flux inversion system. In particular, the annual total errors were consistently less than a ppm, with the true value almost always falling within a quartile of the differences. Globally, the errors are on the order of a few tenths of a ppm. The empirical sink thus enables the study of carbon cycle anomalies, like the 2015-2016 El Niño, whose effect on the atmospheric growth rate is a few ppm or less. By removing 345 the dominant error contribution to our baseline fluxes, the lack of a net sink, the bias-correction also opens the possibility of correcting for errors in other flux components, e.g., fossil fuel emissions, a major, long-term scientific goal needed in the implementation and verification of international climate accords.
There are several benefits to using this approach to bias-correct a system's surface fluxes. In comparison to a flux inversion system, it is exceedingly simple to implement. It is also less susceptible to errors due to the particular transport model, data 350 selection, and error covariance models. For example, constraining the global growth rate alone requires few, if any, assumptions about atmospheric transport or decorrelation times and lengths. When used in simulations or as a prior in an assimilation system, this approach significantly reduces biases due to misspecification of the growth rate. Failing to remove such a bias before assimilating data limits the ability of the assimilation system to account for other signals of interest in the observations (Dee, 2005), e.g., synoptic-scale variations due to passing weather systems, regional and seasonal anomalies due to drought, 355 and changes in anthropogenic emissions. Finally, our approach produces high-resolution fluxes which are often precluded by the computational demands of flux inversion systems.
Since observational estimates of the global growth rate are currently only available at the end of each year, using the empirical sink developed here in an NRT atmospheric monitoring system requires a prediction of the global growth rate. We projected growth rates in 2016 and 2017 based on forecasts and analyses of SST (Jones and Cox, 2005;Betts et al., 2016). The values of 360 2.93 and 2.63 ppm were reasonable estimates of the values measured in the MBL of 2.85±0.09 and 2.15±0.09 ppm which are unavailable until a few months after the year's end. The predicted CO 2 mixing ratios showed comparable skill in reproducing in situ observations. Combined with the future ability to assimilate satellite retrievals of CO 2 lagging real time by just a few days, we expect to be able to monitor and predict growth rates in NRT.
Data availability. The LoFI fluxes are available upon request from the corresponding author. NOAA CarbonTracker, CarbonTracker Europe, and the CAMS flux inversions are all publicly available from their institution's websites. All other data products used in this work (viz., the TRENDY ensemble and the Jena CarboScope flux inversion) must be obtained from the respective project leads.  (Denning et al., 1995). To resolve the diurnal cycle of the fluxes, we first downscale the monthly LoFI fluxes to daily using the the algorithm described below. We stop at daily for all fluxes other than terrestrial NEE, which is the difference of ecosystem respiration (ER) and GPP. Daily terrestrial ER and GPP are downscaled to 3-hourly following the approach of Olsen and Randerson (2004) with the slight modification of starting from daily instead of monthly fluxes. This approach is quite similar to the downscaling approach used by NOAA's CarbonTracker system. It has the advantage 375 of avoiding the noticeable discontinuities at monthly boundaries that are present in the monthly to 3-hourly downscaling, but the disadvantage of possibly missing synoptic scale disturbances that occur over multiple days since the downscaling from monthly to daily uses interpolation.
Fluxes are downscaled to a higher resolution either spatially or temporally by finding the smoothest interpolant that preserves the averages at the original, coarser resolution. This interpolant is found by minimizing a quadratic cost function subject to a 380 linear constraint. The quadratic cost function is the square of the discrete approximation to the Laplacian of the downscaled field. When interpolating in space, we use the spherical Laplacian, which takes the shrinking distance between grid boxes near the poles into account. The linear constraint is that the averages of the downscaled field at the coarser resolution be the original values. It is imposed using Lagrange multiplier, which transforms the problem into an unconstrained quadratic optimization problem on a higher-dimensional space. This approach is used to downscale monthly NEE to daily as described above and is 385 applied and to downscale ocean pCO sw 2 from its native 4 • × 5 • spatial grid to the 0.5 • × 0.625 • grid of MERRA-2 (as shown in Figure A1). vegetation models except LPJ-GUESS, which did not submit monthly results, and the ocean model ensemble is the collection of global annual ocean exchange totals reported in GCP 2018. While comparison to these ensembles is not a true validation of our fluxes and is susceptible to uncertainties in lateral exchanges between the land and ocean, we do expect it to indicate when and where our surface flux product is an outlier compared to other estimates. In particular, we use the ensembles to identify coherent, systematic surface flux errors over wide zonal bands and multiple months. 400 Figure A1. Climatologies (2003Climatologies ( -2015 of ocean exchange from the LoFI fluxes. The first row represents annual averages and each subsequent row represents averages over different seasons. The left column depicts zonal means from the LoFI fluxes (dash-dot red) and the ranges of zonal means from the inversion ensemble (dark grey); and the right column shows seasonal gridded maps of the diagnostic fluxes.
The ensemble comparisons only consider net biospheric exchange (NBE), the sum of net ecosystem exchange (NEE), biomass burning emissions, and all other emissions not due to the combustion of fossil fuels. For the inversion ensemble, we compute NBE by subtracting a common fossil fuel product, which we use in our diagnostic fluxes, from the total surface flux. No distinction is made here between changes in NEE and biomass burning due to natural effects or to anthropogenic activities such as land use change. To avoid issues associated with statistics on small sample sizes or bias due to over representation 405 of certain model configurations, we use only the range of the minimum and maximum values over each ensemble.

A3 QFED versus GFED in LoFI
The main motivation behind replacing GFED in CASA-GFED with QFED is to avoid having to extrapolate biomass burning emissions in the forward processing product. Unlike other flux components, we do not expect the construction of a month-bymonth linear climatology to have much skill for biomass burning. Rather than switching from GFED to QFED in the switch 410 from reanalysis to forward processing fluxes, we chose to always use QFED. This prevents the introduction of a jump in the fluxes and allows a more direct comparison of averages but has the downside of introducing a biomass burning component that is not in balance with the carbon stocks of CASA. In any case, as is shown in Figures A2 and A3, the difference between QFED and GFED biomass burning emissions is quite small at almost every time and place. The only noticeable difference is that the baseline fluxes are slightly lower in the Amazon and Congo rain forests during JJA. Still, this difference is minor in 415 comparison to the empirical sink.  The left column depicts zonal means of the baseline fluxes using GFED and QFED, the range of the inversion ensemble (dark grey), and range of the TRENDY V7 ensemble (light grey); the middle column depicts the distance of the baseline fluxes using GFED to the range of the inversion ensemble; and the right column the difference of the baseline fluxes using QFED with those using GFED, i.e., the difference of QFED with GFED.