An empirical spatiotemporal description of the global surface - atmosphere carbon fluxes: opportunities and data limitations

. Understanding the global carbon (C) cycle is of crucial importance to map current and future climate dynamics relative to global environmental change. A full characterization of C cycling requires detailed information on spatiotemporal 10 patterns of surface-atmosphere fluxes. However, relevant C cycle observations are highly variable in their coverage and reporting standards. Especially problematic is the lack of integration of the carbon dioxide (CO 2 ) exchange of the ocean, inland freshwaters and the land surface with the atmosphere. Here we adopt a data-driven approach to synthesize a wide range of observation-based spatially explicit surface-atmosphere CO 2 fluxes from 2001 and 2010, to identify the state of today’s observational opportunities and data limitations. The considered fluxes include net exchange of open oceans, 15 continental shelves, estuaries, rivers, and lakes, as well as CO 2 fluxes related to net ecosystem productivity, fire emissions, loss of tropical aboveground C, harvested wood and crops, as well as fossil fuel and cement emissions. Spatially explicit CO 2 fluxes are obtained through geostatistical and/or remote sensing-based upscaling; thereby minimizing biophysical or biogeochemical assumptions encoded in process-based models. We estimate a bottom-up net C exchange (NCE) between the surface (land, ocean, and coastal areas) and the atmosphere. Though we provide also global estimates, the primary goal of


Introduction
The global carbon (C) cycle is crucial for sustaining life on Earth (Vernadsky, 1926).Humans have largely modified the C cycle over centuries if not millennia (Ruddiman, 2003;Pongratz et al., 2009).In the Industrial Era, the human-caused perturbation of the C cycle is largely driven by emissions of CO 2 from burning fossil fuel C previously in geological deposits, and changes in land use, which transfer CO 2 from C stocks in the land biosphere to the atmosphere, but can also result into CO 2 removal and increase of land stocks.As those anthropogenic C emissions are partly taken up by oceans and terrestrial ecosystems not affected by land-use change, the different reservoirs of the global C cycle and the fluxes between them change over time (Houghton, 2007).A precise knowledge of the various stocks and fluxes in the C cycle is a prerequisite to monitor these changes and make well-informed predictions under future climate change.
The Global Carbon Project (GCP) has made major efforts in this direction and its annual updates of the global C budget have become a key source of information for the scientific community and policy makers (Le Quéré et al., 2015).The GCP annual C budget quantifies the partitioning of anthropogenic C emissions among the atmosphere, land, and ocean components of the global C cycle, and separates the net land flux into land use change emissions and a so called 'residual land C sink' obtained by difference with other terms of the budget and thus corresponding to the net land-atmosphere CO 2 flux over non land-use affected ecosystems.The budget of the GCP focuses on annual values integrated at the global scale.An important point is that the GCP budget quantifies solely the anthropogenic perturbation of CO 2 fluxes, i.e., it provides information about the fate of anthropogenic CO 2 emissions in natural reservoirs (Ciais et al., 2013).According to the GCP, about 44% of the anthropogenic CO 2 emissions each year stay in the atmosphere, the rest being taken up by the oceans (26%) and land (30%) (Le Quéré et al., 2015).
Recently, a case has been made for a globally policy-relevant integrated C observation and analysis system (Ciais et al., 2014).This system would go beyond the update of global budgets, for which the CO 2 growth rate accurately measured at a single station (e.g.Mauna Loa) is sufficient to constrain the global annual time-space integral of all CO 2 sources and sinks.It proposes to quantify regional CO 2 fluxes with sufficient spatial details to monitor the effectiveness of CO 2 mitigation and to detect and monitor trends of CO 2 losses and gains by land and terrestrial systems.This is partly relevant for monitoring country-level Intended Nationally Determined Contributions (INDCs) to incept a CO 2 emission trajectory consistent with global warming below 2 degrees Celsius (UNFCCC, 2015).In such a policy-relevant C observing system, an uncertainty assessment for each data stream and CO 2 fluxes at different spatial and temporal scales is important to, for instance, identify significant regional emission hotspots and trends in emissions and sinks (Ciais et al., 2014).
The steadily increasing number of Earth observations, in particular since the start of the satellite era, has improved our knowledge of the Earth system (Berger et al., 2012;Tatem et al., 2008).Especially C cycle science has benefited from globally available satellite observations and community efforts to unify in-situ observational networks such as FLUXNET on land (Baldocchi, 2014), the Surface Ocean CO 2 Atlas (SOCAT) (Bakker et al., 2014) and more recently CO 2 outgasing from lakes and rivers (Raymond et al., 2013).Combining these available point measurements of either CO 2 fluxes (e.g., from eddy-covariance towers on land), or variables that can be directly related to CO 2 fluxes (e.g., pCO 2 over aquatic surfaces) with climate fields and remotely sensed variables (e.g., vegetation greenness), provides a basis to robustly upscale surfaceatmosphere CO 2 exchange to larger areas using statistical models (Jung et al., 2011;Rödenbeck et al., 2015).
In this study we aim at characterizing the ability of current C cycle observations on ground for quantifying a spatiotemporally explicit picture of the net CO 2 exchange between the Earth's surface (terrestrial and aquatic) and the atmosphere (NCE).Unlike the GCP global budget of anthropogenic CO 2 , we consider here the full contemporary exchange of surface-atmosphere CO 2 fluxes.We focus our analysis on fluxes that can be directly derived from observations.That is, we use data-driven empirical models instead of process-based models that are only indirectly constrained by observations.Further, we only consider 'bottom-up' estimates derived from measurements at the Earth's surface or from satellites.
Inversions, which largely rely on atmospheric measurements in combination with a transport model, are not directly included but used for comparison.The goal of this analysis is to test the up-scaling of local flux-related observations to regional and global budgets, and point out the limitations of the current observational networks and data-driven models used to interpolate point-scale CO 2 fluxes across larger scales, for quantifying the most important CO 2 fluxes exchanged between the Earth's surface and the atmosphere.
One of the major innovations of this study is combining data-driven estimates of oceanic, inland waters and terrestrial ecosystems CO 2 exchange and providing spatially explicit maps of the CO 2 exchange between the surface and the atmosphere at a monthly scale for the decade 2001-2010.At the same time, by adding emissions from fossil fuels and cement production and comparing with the annual growth rate of CO 2 , we identify the limits of a C budget purely driven by surface data.We characterize regions in which surface-atmosphere CO 2 fluxes are most uncertain based on the currently available data and the models used for upscaling, and thus point out regions where either more observations or a better understanding of the processes are necessary.It is not the primary goal of this study to provide the best global CO 2 flux inventory, but rather to identify the key uncertainties and observational shortcomings that need to be prioritized in the expansion of in-situ observatories.
The paper is structured as follows.In Sect. 2 we introduce the different data streams used in the analysis, including spatially explicit estimates of aquatic and terrestrial CO 2 exchange.In Sect. 3 we present the resulting combined synthesis as global maps, regionally aggregated fluxes, absolute and relative uncertainties, latitudinal averages and seasonal cycles.Sect. 4 addresses the benefits and limits of the current observational system for constraining global net CO 2 fluxes.Sect. 5 provides an outlook on future requirements to achieve better observation-based net CO 2 flux estimates and discusses the necessity for more consistent uncertainty estimates.

Data and Methods
We collected ensembles of data-driven estimates of the net CO 2 exchange between the Earth's surface and the atmosphere (NCE) for the major subsystems of the Earth from 2001-2010 (Table 1).Each dataset was aggregated to 1 x 1 degree spatial resolution.All datasets have an original spatial resolution of at least 1 x 1 degree such that no information was lost through re-gridding.The temporal resolution is monthly.For datasets that were only available at yearly time scale or once over the complete time period (Table 1), we distributed fluxes evenly across all months.In this synthesis, we include net CO 2 exchange from open oceans, continental shelves, estuaries, rivers, lakes, and terrestrial ecosystems, which we combine with estimates of fossil fuel and cement emissions (FF).The terrestrial ecosystem component accounts for fire emissions (Fire), loss of tropical above-ground biomass assumed to be released as CO 2 to the atmosphere (E LUC ), emissions of the CO 2 contained in harvested wood (Wood) and crops (Crops), and Net Ecosystem Productivity (NEP).We combine fluxes from oceans, shelves, and estuaries into a homogeneous marine flux product in order to account for overlapping or missing regions from the different aquatic products (Marine, Sect.2.2.6).We further compare the net CO 2 exchange derived from the combination of all the above products with the growth rate of atmospheric CO 2 (CGR).Data scarcity precludes including all

Uncertainty estimation and propagation
To estimate NCE including spatiotemporally explicit uncertainties, we combine randomly drawn ensemble members from all of the 9 datasets contributing to NCE (Eq.1).With the available realizations, we could in principle create 10x10x50x8x10x2 = 800000 spatiotemporal explicit estimates of NCE (see Table 2 for the available number of realizations per dataset).From these 800000 possible NCE estimates we randomly select 200 to construct the NCE ensemble, which is used in the rest of the paper.This resampling approach is illustrated in Figure 1 and ensures a consistent propagation of spatiotemporally correlated uncertainties.For instance, by aggregating each member of NCE to the desired region and estimating uncertainty through the 200 members, we can compute regional uncertainties.In addition, we computed mean fluxes, uncertainty (defined as one standard deviation (SD) of the annual mean across all realizations), interannual variability (IAV, defined here as one SD of annual means across all available years) and the coefficient of variation (CV = IAV/mean) for each of the 9 flux terms in Eq. 1.

Oceans
For the global open ocean flux estimate we used two complementary data-driven estimates (Table 1).Both approaches computed maps of the sea surface partial pressure of CO 2 (pCO 2 ).They relied on the surface ocean CO 2 observations from the SOCATv2 database (Bakker et al., 2014) and filled data gaps by either establishing relationships between auxiliary driver data and observations, which can then be applied to extrapolate pCO 2 in regions without data coverage (SOM-FFN, Landschützer et al., 2014), or by assimilating the available observations in a mass-balance model of the mixed layer and directly interpolating data gaps (Jena CarboScope mixed-layer scheme oc_v1.2,Rödenbeck et al., 2014).To test the established predictor-target relationship, the SOM-FFN method holds back a certain fraction of the observations proportional to the methods degrees of freedom for internal validation.Repeating this relationship building process and withholding different sets of validation data has created the 5 ensemble members used for this study.For the Jena CarboScope mixedlayer scheme, we used the 5 sensitivity cases with changes in correlation length etc. as described in Rödenbeck et al (2014).
The pCO 2 fields of both methods have been validated against independent observations (Landschützer et al., 2014;Landschützer et al., 2015;Rödenbeck et al., 2014) and were compared with other complementary data based interpolation methods (Rödenbeck et al., 2015), illustrating their good performance in reconstructing interannual variation.
Both methods calculate the air-sea flux using a bulk formulation of the air-sea CO 2 transfer, driven by the air-sea pCO 2 difference (∆pCO 2 ) (Jähne et al., 1987) and a quadratic dependence of the wind speed at a height of 10 meters (Wanninkhof, 1992) updating the gas transfer coefficient to fit a mean transfer velocity of 16 cm per hour following Wanninkhof (2013).
High-resolution wind speeds at 10 meters are calculated from the u and v wind components of the ERA-interim wind speed analysis (Dee et al., 2011) and atmospheric pCO2 fields, required to calculate the ∆pCO 2, are calculated are estimated from the GLOBALVIEW-CO2, 2012 Marine Boundary Layer CO2 product.

Shelves
For continental shelf seas we derived the ∆pCO 2 from 3x10 6 surface pCO2 measurements extracted from the SOCATv2 database (Bakker et al., 2014) and observational atmospheric pCO 2 data (GLOBALVIEW-CO2, 2012).The local CO 2 airsea flux values were then obtained using a wind-dependent quadratic formulation parameterized as in Wanninkhof et al. (2013) and wind speeds extracted from a cross-calibrated multiplatform (CCMP) high-resolution data product for ocean surface winds (Atlas et al., 2011).The resulting local fluxes were then integrated spatially over 150 coastal regions (COSCATs -COastal Segmentation and related CATchments; Laruelle et al. (2013); Meybeck et al. (2006)) using distinct integration methods depending on the data density (Laruelle et al., 2014).In addition, a temporal integration was also performed at the monthly, seasonal or yearly time scale depending on the data coverage.These temporally and regionally averaged air-sea CO 2 fluxes were then disaggregated using a 1-degree resolution map excluding land areas and open ocean waters using the shelf break as outer limit (Laruelle et al., 2014).

Estuaries
The CO 2 emissions from estuaries were derived from 161 annually averaged local CO 2 air-water exchange rates reported in the literature (Laruelle et al., 2013).The data were allocated to one of the 45 coastal MARCATS regions (MARgins and CATchments Segmentation) defined in Laruelle et al. (2013) and further categorized among the 4 dominant estuarine types (i.e., small deltas, tidal systems, lagoons, fjords, see (Dürr et al., 2011)) to calculate regionally-averaged, type specific CO 2 emission rates.In MARCATS regions devoid of estuarine data, the global average type-dependent air-water CO 2 flux was used from Laruelle et al. (2013).These flux densities were then multiplied by the estuarine surface areas for each type, estimated at 1-degree resolution from the length of the coastline and a type-specific length to estuarine surface ratio (Dürr et al., 2011).

Marine
We combined open oceans, shelves, and estuaries to a consistent marine product.For pixels with observations from multiple products (e.g., estuaries and oceans) we follow a "priority rule" whereby the shelves, estuaries, or oceans observation value only (in that order) is retained.Empty pixels are gap-filled with 3 x 3 mean window.This same filter is also applied to the rest of the merged dataset to smooth out hard borders between the different estimates.This application does not significantly change the overall flux estimates, but arguably results in a more realistic interface.Note that in the merged Marine product, uncertainty and IAV could only be assessed for the ocean flux.

Rivers
Estimates of CO 2 evasion from streams and rivers were derived from a spatially explicit, empirical model of river water pCO 2 and global maps of stream surface areas and gas exchange velocities at a resolution of 0.5 degree (Lauerwald et al., 2015).The empirical pCO 2 model was trained on 1182 river catchments from the GLORICH database (Hartmann et al., 2014) for which averages of pCO 2 could be calculated.Steepness of terrain, terrestrial net primary production, average air temperature as well as population density were identified as predictors (R 2 =0.47).The global maps of stream surface area and gas exchange velocities were obtained by a GIS-based application of published empirical scaling laws (Raymond et al., 2013;Raymond et al., 2012) using topography (Lehner et al., 2008) and runoff (Fekete et al., 2002).The CO 2 evasion was calculated as product of water-air pCO 2 gradient (assuming an atmospheric pCO 2 of 390 µatm), river surface areas, and gas exchange velocities.A Monte-Carlo simulation based on standard errors of the predictors in the pCO 2 model and uncertainty ranges for estimates of stream surface area and gas exchange velocity was run to produce 50 CO 2 evasion estimates.

Lakes
Estimates of CO 2 evasion from lakes and reservoirs were taken from Raymond et al. (2013), which reports average lake pCO 2 , total lake/reservoir surface area, and total CO 2 evasion for 231 COSCAT regions (including endorheic regions).For the total lake/reservoirs surface area, data from the Global Lakes and Wetland Data base (GLWD, Lehner and Döll, 2004) were combined with an estimate for small lakes and reservoirs not represented in the GLWD using a scaling law.Here, we used the GLWD data to downscale the estimates of Raymond et al. (2013) to a continuous 1-degree resolution.For this purpose, we combined a uniform air-water CO 2 flux (per unit surface area) within each COSCAT region with a spatially explicit estimate of the lakes/reservoirs surface at this resolution.The small lakes/reservoirs not represented in the GLWD were assumed evenly distributed over the COSCAT area.

NEP
We used empirical, machine learning based products from FLUXCOM (www.fluxcom.org)for net ecosystem productivity (NEP), derived from more than 200 FLUXNET sites and exclusively remote-sensing based predictor variables ("FLUCOM-RS", see Tramontana et al., 2016).The eight machine learning methods used here include artificial neural networks, four variants of model or regression tree ensembles, kernel methods (support vector machines, kernel ridge regression), and multivariate adaptive regression splines (Tramontana et al., 2016).All methods were trained on 8-daily tower based NEP estimates.

Crops
About 42% of global crop biomass is harvested, transported, and respired off site (Wolf et al., 2015a).The impact of this lateral C transport on fluxes can be seen at the country scale in the form of import and exports, but even more so at subregional scales where the movement of crop biomass to feed livestock and humans is evident (Hayes et al., 2012;West et al., 2011).To capture the spatial distribution of CO 2 fluxes from agricultural harvest, we used livestock and human CO 2 emissions estimates (Wolf et al., 2015b) that are available from 2005-2011 at 0.05 degree spatial resolution.CO 2 that has previously been taken up from the atmosphere by the harvested biomass of crops is included in the NEP estimates from FLUXCOM.We aggregated best estimates of the data to 1 degree, added all uncertainty estimates within one 1 degree pixel and used them as estimates for one standard deviation on the new 1 degree grid.Assuming Gaussian distributed errors we sampled 1000 values at each pixel and used 10 maps of the 5 th , 15 th ,…, 95 th quantiles as different ensemble members.Data was then linearly extrapolated back to 2001-2004.In a final step, and because it is not known in which months the emissions occur, we further distributed the annual estimates equally across all 12 months.

Wood
We used globally gridded forest harvesting data around year 2000 as described in the Supplementary Information S1.These data include fuelwood and roundwood harvested volumes in m 3 .We translated wood volumes into units of C using a value of 0.275 MgC/m 3 from FAO (http://www.fao.org/docrep/w4095e/w4095e06.htm),assuming wood density of 0.55 t/m 3 .To avoid double counting wood harvest with aboveground biomass loss in tropical areas exposed to land use change, we use wood harvesting data only in locations where the amount of harvested wood (in C) exceeds E LUC (Sect.2.3.4).We assume that 100% of the harvested wood is respired back to the atmosphere within a year, thus assuming no change in C stock of wood products and constant harvesting rates across years.However, C contained in harvested wood is usually emitted at a different location than where the harvest took place.We thus incorporated lateral shifts of harvested wood by redistributing wood harvest according to the consumption of wood as explained in the Supplementary Information S1 (see also Fig. S2).

E LUC
We used two estimates for CO 2 fluxes due to tropical deforestation and degradation.It is assumed here that 100% of biomass loss is converted to a CO 2 flux being released instantly (within a year) to the atmosphere.In reality, a fraction of lost tropical biomass decays in ecosystems (belowground biomass and slash) and a fraction is used in wood products of various lifetime.
However, slash is decomposed fast and biomass from deforested areas is transformed on average to short-lived products (≈ 5 years after Earles et al. ( 2012)).
1) Gross tropical deforestation emissions were taken from Harris et al. (2012).They represent total (above-and belowground) C loss from gross forest cover loss in the tropical regions due to human or natural causes (e.g.disturbances without forest recovery) for the period of 2000-2005.
2) More recent estimates of aboveground C loss in the tropics from stand-replacement disturbance of forest cover due to human or natural causes were provided by Tyukavina et al. (2015).Sample-based estimates of mean 2000-2012 aboveground C loss for each 30-m resolution forest C stratum were attributed to all pixels of the corresponding stratum and averaged to the 1x1 degree resolution.
We used E LUC only in those pixels where the average of 1) and 2) exceeds wood harvesting (Sect.2.3.3).

Fire
We used fire emissions from the Global Fire Emissions Database version 4 with small fires (GFED4s, http://www.globalfiredata.org)based on burned area from Giglio et al. (2013) and Randerson et al. (2012) and an updated version of the biogeochemical modelling framework of van der Werf et al. ( 2010) to convert burned area to C emissions.We included all fire types except tropical deforestation and degradation fires, which are included in E LUC and should thus not be counted twice (Sect.2.3.4).For an earlier version of fire emissions (GFED3) a Monte Carlo simulations indicated an uncertainty of about 20% (1σ) for continental-scale estimates but these estimates turned out to be not very reliable (van der Werf et al., 2017).For example, the inclusion of small fire burned led to an increase in burned area exceeding the previously assumed uncertainty and the current version therefore has no uncertainty assessment at pixel level.Note that GFED fire emissions depend on estimates of net primary production, and combustion factors as computed by the CASA model.

FF
We use the IER-EDGARv4.2product for fossil fuel and cement emissions, which was derived within the CARBONES project by the Institute für Energiewirtschaft und Rationelle Energieanwendung (IER).It is based on the Edgar v4.2 fossil fuel spatial distribution (with the highest spatial resolution of 0.1 x 0.1 degree) and uses national consumption and global production statistics.Based on the sectorial distinguished EDGARv4.2 emissions, sector-specific and country specific temporal profiles were included.A detailed description of the construction of the product is given at http://www.carbones.eu/wcmqs/project/ccdas/#Fossil%20Fuel.It is important to note that FF emissions here are not observation based as the IER-EDGARv4.2product is partly based on national estimates from official inventories reported by countries to the UNFCCC.

Atmospheric growth rate
We used the atmospheric rate of change of CO 2, which is equal to the space and time integral of all emissions and sinks at the surface, using the calculations made by the GCP (Le Quéré et al., 2015).These calculations and are based on the global growth rate of atmospheric CO 2 (CGR) provided by the US National Oceanic and Atmospheric Administration Earth System Research Laboratory (NOAA/ESRL) and were derived from multiple stations selected from the marine boundary layer sites with well mixed background air (Ballantyne et al., 2012;Masarie and Tans, 1995).They applied conversion from concentrations to carbon mass is 1 ppm = 2.12 PgC (Prather et al., 2012).

Inversions
For a comparison of yearly variability, spatial patterns and latitudinal bands, we used annual means of 10 inversions collected in Peylin et al. (2013), available at the same spatial and temporal resolution.The mean and uncertainty for each year is taken over all available inversions for that year, as not all inversions were available until 2010.Atmospheric CO 2 inversions estimate surface CO 2 fluxes such that they best fit observed atmospheric constraints.They usually rely on prior information provided by terrestrial and oceanic biogeochemical models but are mostly independent from the bottom-up datasets included in the present synthesis.They further use FF as an input and then provide the surface-atmosphere flux excluding FF.

Global net carbon exchange
Mean fluxes, their uncertainties, interannual variability (IAV), and CV (the mean-normalized IAV) for all individual fluxes contributing to NCE are presented in  Quéré et al., 2015).Thus, there is a large mismatch with our NCE, which over-estimates the CO 2 sink at the surface by 9.7±2.0PgC / year.This highlights that our observation-based NCE is biased towards a too large sink.Potential reasons for this mismatch are discussed in Sect. 4. For most fluxes, uncertainty estimates strongly exceed IAV (Table 2).Interestingly, process-based models, which are only indirectly constrained by observations, provide an NCE that matches roughly the CO 2 growth rate (Le Quéré et al., 2015).Developers of process-based models have access to CO 2 growth rate data and may be in the position to tune their models so that they give realistic NCE values, whereas in our bottom-up approach, we conducted a blind up-scaling of ground measurements without trying to match the CO 2 growth rate.

Spatial patterns of net carbon exchange
The 200-member NCE ensemble and the uncertainty distribution of each flux component enables us to provide a best estimate for a gridded average surface-atmosphere CO 2 flux map for the time period 2001-2010 (Figure 3a).According to these estimates, tropical land areas are a larger CO 2 sink than the mid-latitudes despite the visible forest bands in North America and Russia that function as sinks.In contrast, the high latitudes indicate a relatively small source.In the ocean, these patterns are reversed, with sources in the tropics and a sink in the mid-latitudes.Clearly, there is a strong land-sea contrast and land NCE is much higher in magnitude compared to ocean NCE.In areas with high human population densities and active industry (Europe, Eastern China, US, South Africa), emissions from fossil fuels and cement production clearly dominate over land CO 2 fluxes.
Absolute uncertainty of NCE generally scales with the mean flux and is highest in the most productive areas over land (Amazon basin, Congo basin, Indonesia; Figure 3b).Due to the small contribution of the oceans, absolute uncertainties are barely discernible there.Although gross air-sea exchange fluxes have typical uncertainties of more than 20%, their differences are determined from independent measurements with a much higher accuracy (Ciais et al., 2013).
Relative uncertainties however show very distinct patterns (Figure 3c).These are high on land in semi-arid and arid, and in mountainous regions (i.e., rather unproductive areas with near-zero mean) such as Australia, the Middle East, the Midwest US, the Sahel, South Africa, the Andes, and around the Tibetan Plateau.Marine-atmosphere CO 2 exchange is most uncertain in relative terms in the Bay of Bengal and in the Southern Ocean, which is known to be under-sampled, and where the two data-driven NCE fluxes show substantial regional patterns (Landschützer et al., 2014;Rödenbeck et al., 2014).In addition, linear features with high relative uncertainty are visible, especially in the Southern Hemisphere.These are related to the borders of the clusters used for deriving homogenous regions of sea-air exchange in one of the ocean-exchange products, which result in this product in strong spatial gradients in the sea surface pCO 2 (Landschützer et al., 2014).Relative uncertainties are mostly below 100% for the median across latitudinal bands (Figure 3c).Only in the Southern Ocean the relative uncertainty is substantially higher, reflecting difficulties in reconstructing seasonal to interannual variabilities with sparse observational constraints (Landschützer et al., 2014;Rödenbeck et al., 2014).Nevertheless, Landschützer et al. (2015) have shown that there is a better agreement between the estimates of Landschützer et al. (2014) and Rödenbeck et al. (2014) when low frequency variability, such as decadal variability, is analysed.
Averaged over latitudinal bands, the tropics are clearly a CO 2 sink (Figure 4a), a feature of the FLUXCOM models used for NEP, whereas mid-latitudes form a net CO 2 source, mostly due to fossil fuel and cement emissions surpassing natural CO 2 sinks.This latitudinal pattern is strongly driven by the terrestrial fluxes (Figure 4b).Marine and land aquatic CO 2 exchange in turn is about 4 times smaller in magnitude and shows CO 2 sources in the tropics and CO 2 sinks in the extratropics (Figure 4c).The aquatic CO 2 source in the tropics is not only the result of the ocean air-sea exchange, but also of the very intense river outgassing in low latitude regions (Lauerwald et al., 2015).NCE in the mid-latitudes is dominated by fossil fuel emissions (blue line in Figure 4d shows NCE-FF).FF have little contribution in the tropics and the high-latitudes but offset land and ocean CO 2 sinks in the northern mid-latitudes so that the net CO 2 balance of this latitude band is a net CO 2 source.
We use the land cover map of FLUXCOM to identify tropical forests (all pixels where broadleaved evergreen trees dominate).Tropical forest, which covers about 3.5% of the Earth's surface, are allocated a CO 2 sink of -5.0±0.6 PgC / year, which is unrealistic, if compared to e.g.forest biomass inventories (Pan et al., 2011)

3.3
Net carbon exchange over the RECCAP regions

RECCAP over land
Here we compare our NCE estimates over land with largely independent estimates of net ecosystem exchange (NEE) over continental-scale regions collected in RECCAP (REgional Carbon Cycle Assessment and Processes).The RECCAP budgets were based on inventories, and in some instances on process models results.NCE components between RECCAP and this study that are not independent from each other are fire emissions, and FF emissions.For E LUC , RECCAP publications used regional datasets or bookkeeping models, that are independent from estimates gathered in Secti.2.3.4.These regions include North America (NA, King et al., 2015) 2014)).

RECCAP over ocean
We compare our estimates of mean annual C exchange over the ocean with estimates from the RECCAP initiative.The Pacific Ocean is divided into North Pacific extratropics (NP), Tropical Pacific (TP) and South Pacific extratropics (SP) (Ishii et al., 2014).The Atlantic Ocean is divided into Arctic Ocean (AR), Northern Subtropics (NS), Equatorial (EQ) and Southern Subtropics (SS) (Schuster et al., 2013).Further, there are estimates for the Northern (NI) and Southern Indian Ocean (SI) (Sarma et al., 2013) as well as for the Southern Ocean (SO) (Lenton et al., 2013).The RECCAP estimates of NCE over oceans are independent from the two estimates that we use in this study (Sect.2.2.1).Overall, the estimates from both sources agree very well (Fig. 6) and show ocean net C release in tropical regions (TP, EQ and NI) and net C uptake in all other regions.In SO our estimates predict a smaller sink compared to the RECCAP estimates, a difference probably owing to the weak observational constraint (Landschützer et al., 2014;Rödenbeck et al., 2014).Our estimates generally have smaller uncertainty ranges, which is because i) the RECCAP studies include many more approaches (including processbased models, atmospheric and ocean inversions) in their estimates and ii) in our analysis we include the uncertainty from the ocean pCO2 products and their realizations but do not account for the uncertainty in the kinetic gas transfer.

Comparison with inversions
We compare NCE without FF (NCE-FF) with annual values from 10 inversions estimating the surface-atmosphere CO 2 flux without FF (Peylin et al., 2013).While both estimates agree well in the mid-latitudes, they show opposite patterns in the tropics and northern high latitudes (Figure 4d).The estimates of NEP in our NCE-FF probably have a substantial bias towards too much uptake over tropical land (Sect.4.1).The comparison suggests that C fluxes are comparably well constrained in the mid-latitudes where bottom-up and top-down approaches agree.Similar results have been obtained in a comparison of a bottom-up upscaling approach with a more recent inversion based on CO 2 concentration data from the Greenhouse gases Observing SATellite (GOSAT, Kondo et al., 2015).The temporal evolution between both estimates show little agreement except the trend towards more net C uptake by the Earth's surface (Figure 7).The comparison suggests that NCE-FF estimated from this study has lower interannual variability compared to inversion estimates.Uncertainties are very high for our NCE-FF.In addition, the mean annual C uptake in our estimates is nearly 10 PgC/yr higher than for inversions.

Monthly variability and mean seasonal cycle
NCE in the Northern hemisphere (NH) exhibits a much stronger mean seasonal cycle compared to the Southern hemisphere (SH), ranging from a net C uptake of nearly 2 PgC (per month) in July to a net C release of about 0.9 PgC in December and January (Figure 8).The SH is always a net C sink, ranging between slightly under 0.8 PgC uptake in January to roughly 0.1 PgC in August and September.This illustrates the 'breathing of the Earth', that is, vegetation activity largely follows the annual cycle of the sun.NH NCE is strongly offset by fossil fuel emissions.The uncertainties for the SH seasonal cycle are generally much lower than for the NH fluxes due to the larger contribution of the latter to the overall flux pattern, which is related to the distribution of land areas.If compared to inversions, we find that both estimates only match in the summer of the NH.In all other months and in the SH, our NCE estimates show a consistently much larger surface C sink.In addition, the NCE estimates from this synthesis show a smaller amplitude of the mean seasonal cycle compared to the inversions.The difference in amplitude of the mean seasonal cycle is on average 0.7 PgC for the NH and 0.4 PgC for the SH.

Current limitations of a bottom-up spatiotemporal assessment of net carbon exchange
Our study shows that today's spatiotemporally explicit and independent bottom-up observation-driven estimates of surfaceatmosphere CO 2 exchange suffer from large bias, such that they do not match the global NCE well constrained from the CO 2 growth rate.This statement is not downgrading the advances in the area, but rather a systematic reflection of the state of current research and monitoring.In fact, at the regional scale, those estimates are often well constrained and may be used for model-data integration studies and validation purposes.The regions where observation-driven CO 2 exchange is constrained the best include Europe, Russia, South Asia, East Asia, Australia and all oceanic regions except the Southern Ocean.The most likely candidate for inducing the mismatch between data-driven estimates and the atmospheric CO 2 growth rate is terrestrial NEP.In particular, tropical NEP estimates suggest a too large tropical sink.In the following sections, we discuss (i) the possible reasons for the large bias in NEP (Sect.4.1), (ii) which fluxes are missing in our synthesis (Sect.4.2), (iii) how this synthesis dataset could be used for model-data fusion (Sect.4.3), (iv) uncertainties in fire emissions (Sect.4.4), and (v) the impact of missing seasonal cycles in some of the datasets (Sect.4.5).

Difficulties in estimating NEP over land
Correctly predicting NEP from remote sensing requires establishing universal relationships between those predictors and respiratory processes (Jägermeyr et al., 2014;Tramontana et al., 2016).However, predicting such processes still poses major challenges to researchers (Trumbore, 2006).The CO 2 flux related to heterotrophic decomposition processes, for instance, relates to factors controlling biological activity via temperature, moisture availability, and the decomposable substrate material.The question how soil respiration or total ecosystem respiration depends on these variables is not yet entirely understood.Advancing our knowledge on these processes is challenging due to both a lack of theory of respiration and the difficulty of obtaining relevant data to test models (Trumbore, 2006).
In addition to a good theory for respiration, information on disturbance history (e.g., time since last fire) and forest age would improve the upscaling of NEP from sites to regions (Ciais et al., 2014).Disturbances that cause physical damage to vegetation properties tend to temporarily increase respiration and reduce photosynthesis and thus alter the balance between gross C uptake and release.Disturbed ecosystems are thus initially assumed to be strong C sources until plant production recovers.However, how these regrowth processes compensate a given disturbance regime cannot yet be quantified at global scales, as the area covered by disturbed ecosystems is variable and unknown (Ciais et al., 2014).For example, regrowth of vegetation after fires and other disturbances is not well sampled neither in the FLUXNET stations nor in the set of predictors used by the FLUXCOM models and is assumed to be implicit in our NEP estimate.Furthermore, management can have strong effects on annual NEP of croplands, which form large parts of the land surface (Jung et al., 2011).However, not all of the relevant predictors (i.e.disturbance maps, management practices, soil moisture) are currently available to be included in empirical upscaling exercises (Tramontana et al., 2016).
In addition to the above difficulties, some regions are undersampled by eddy-covariance towers and thus NEP is not well constrained.This is the case for tropical forests and the northern high latitudes.In the tropics, undersampling leads to a large overestimation of net CO 2 uptake in comparison to inversion and forest inventories (Peylin et al., 2013;Pan et al., 2011) whereas in the high latitudes it leads to a comparably large CO 2 release (Figure 4).
Given the difference between NCE and inversions in the tropics (Figure 4), we can assume that a bias of FLUXCOM NEP towards a too high CO 2 sink is the main reason why the C budget is not closed in our approach.This raises the question why upscaled NEP has such a strong systematic bias towards a sink, particularly in the tropics (see also Jung et al., 2011).We suspect that the eddy-covariance towers collected in FLUXNET, which provide the empirical basis for the global data driven estimates (see Sect. 2.3.1)do not represent the different age classes of forests very well.For instance, young and regrowing forests with a generally higher-than average NEP are possibly overrepresented in FLUXNET.However, such an agedependency (Amiro et al., 2010;Coursolle et al., 2012;Hyvönen et al., 2007;Magnani et al., 2007) has not yet been included in global upscaling of NEP.This hypothesis should be tested in future upscaling exercises.

Missing fluxes
Due to a focus on spatially explicit maps, not all known fluxes between land surface and atmosphere are considered in our analysis.We assume that including the following fluxes may have an influence on the regional and global flux estimates (estimates of the flux magnitude are given in brackets if known): • Emissions from biogenic volatile organic compounds (VOCs) amount to approximately 0.76PgC / year globally (Sindelarova et al., 2014) • CO 2 emissions from wetlands, estimated globally at around 2.1 PgC /year (Aufdenkampe et al., 2011) • CH 4 emissions from biogenic sources and animals • Crop residues burning in households All known missing fluxes add up to an additional C release of about 4 PgC / yr.Although substantial, they do not cover the mismatch of more than 9 PgC / yr by far (Sect.3.1).However, they would suffice to close the budget if tropical forests are assumed to be C neutral (tropical forests are responsible for a net C sink of about 5 PgC / year, Sect.3.2).This significant amount of missing fluxes prohibits constraining FLUXCOM runs with all the remaining fluxes.In other words, we cannot be certain of the bias in upscaled NEP as long as the major fluxes are not quantified in a spatially explicit manner.Emissions from VOCs and wetlands should thus receive particular attention if a consistent spatiotemporal picture of vertical CO 2 exchange is to be obtained.

Uncertainty estimates and model-data fusion
Our uncertainty estimates of ocean and land C exchange likely underestimate the true uncertainty.In particular, Landschützer et al. (2014) estimated that the choice of sea-air gas transfer formulation (also including other relationships than quadratic) and the pCO 2 mapping mismatch lead to an uncertainty of 37% for the global average over sea-air exchange between 1998-2011, with the majority of this uncertainty stemming from the gas transfer formulation.Furthermore, the uncertainty of NEP is likely underestimated because all upscaling methods in FLUXCOM use the same set of predictors (Tramontana et al., 2016).Hence, the uncertainty estimates only cover the uncertainty related to the upscaling method but do not contain uncertainties related to the selection of predictors or observational uncertainty of the predictors themselves.
A comprehensive spatiotemporally explicit bottom-up estimate of NCE can be a powerful ingredient for model-data integration exercises (Rayner et al., 2005).Yet, model-data integration requires uncertainty characteristics of all used data streams (Raupach et al., 2005).Furthermore, it is important that uncertainties can be described in terms of random errors (Ciais et al., 2014).Error estimates at the local or regional level are difficult to use if no spatial error covariance matrix is available.The uncertainty analysis presented in this study obtained through Monte Carlo sampling aims to be of use for modal-data-integration studies.Errors are automatically propagated through different spatial resolutions by aggregating the individual ensembles of NCE.Naturally, efforts should be made to obtain error estimates for all integrated datasets (i.e., Wood, Fires, Shelves, Estuaries, and Lakes).Nevertheless, this first integrated NCE estimate offers new possibilities for approaches such as the Carbon Cycle Data Assimilation System (CCDAS, Rayner et al., 2005), by not only providing a full spatiotemporal grid of fluxes, but also a transparent and consistent error propagation scheme.This can have also practical applications, for instance for designing new measurement campaigns in regions with high uncertainties to reduce knowledge gaps in the global CO 2 fluxes.

Uncertainties in fire emissions
Fire emission estimates combine satellite-based fire data with ecosystems models.Uncertainties in global fire emission estimates are substantial and different fire products vary largely by location, vegetation type and fire weather (Ciais et al., 2014;French et al., 2011).
While GFED4 burned area estimates come with regional uncertainty estimates (Giglio et al., 2013), the actual uncertainty of C emissions from fires are probably much larger, on the order of 50% (van der Werf et al., 2017).The uncertainties of fire emission estimates depend regionally and temporally on the various input data sets such as burned area, small fire burned area, biomass loadings, and combustion completeness.Better quantifying this uncertainty requires an assessment of the burned area estimates as well as new field data on fuel consumption and emission factors.In this study we cannot propagate this uncertainty into the NCE estimates as this would require spatiotemporal error covariance matrices.

Seasonality for coastal and inland waters, wood and crop harvest emissions
Recently, major steps have been undertaken to resolve the spatial variability of coastal and inland water CO 2 fluxes (Laruelle et al., 2013;Laruelle et al., 2014;Lauerwald et al., 2015;Raymond et al., 2013).Estimates of the seasonal variation in these fluxes are necessary to better constrain seasonal variations in NCE.For inland waters, seasonality has so far only been investigated at regional scale (Laruelle et al., 2015;Richey et al., 2002).For shelves some seasonal estimates are currently available in temperate and high latitudes, indicating that net C uptake is highest in spring whereas C release is highest in summer (Laruelle et al., 2014(Laruelle et al., , 2017)).These estimates indicate that seasonal differences in shelf net C exchange are as high as the annually integrated latitudinal gradient.An analysis performed over Atlantic shelves suggests that the seasonal variability in the air-sea CO 2 exchange is most pronounced over temperate latitudes.In these regions, shelves generally behave as strong CO 2 sinks in winter and spring, partly sustained by CO 2 fixation during the spring phytoplankton bloom, but can become mild CO 2 sources to the atmosphere in summer due to the effect of temperature-driven decrease CO 2 solubility in water (Laruelle et al., 2014).Such behaviour is consistent with that of the open ocean at similar latitudes (Takahashi et al., 2009).In the continental shelves surrounding other oceanic basins, however, a recent study suggests more complex seasonal patterns involving the contributions of processes other than temperature to the seasonality of coastal pCO 2 (Laruelle et al., 2017).
Biogenic C emissions related to tropical aboveground biomass loss as well as crop and wood harvest are equally distributed across months in this study.When exactly C emissions from humans and livestock occur is difficult to predict and would require more detailed consumption data (Wolf et al., 2015a).

Conclusions
From the presented synthesis, we draw the following main conclusions: i) Current estimates of surface-atmosphere CO 2 exchanges that are spatiotemporally explicit and entirely driven by surface observation are not sufficiently well constrained to close the C budget at the global scale.The data-driven estimates show a large bias towards too much C uptake by the Earth surface of nearly 10 PgC / year. ii) The most likely candidate for inducing the mismatch between data-driven surface-atmosphere CO 2 exchange and the atmospheric CO 2 growth rate is land NEP.In particular tropical NEP estimates appear to be strongly overestimated (too large land sink).Understanding this bias will help designing better upscaling approaches (e.g., by including currently missing relevant predictors) and pinpointing variables that need to be (better) monitored in the future.
iii) Regionally, the estimates of NCE are partly well constrained and may be used for model-data integration studies, validation of models, and evaluating claims and potentials of net C uptake within the framework of the Paris agreement (UNFCCC, 2015).These regions include Europe, Russia, South Asia, East Asia, Australia and most oceanic regions.Better constraining C fluxes in regions with currently high uncertainties should be a priority of future research.
known vertical CO 2 fluxes in this study.Missing fluxes include geological CO 2 fluxes, erosion related fluxes, non-CO 2 fluxes, wood product pools decay, and biofuel burning.Combining all fluxes, the overall net CO 2 exchange (NCE) between the Earth's surface and the atmosphere is given as: NCE = Marine + Lakes + Rivers -NEP + Crops + Wood + E LUC + Fire + FF. (1) All units were transformed into fluxes of C per unit time.If all CO 2 fluxes were included, NCE would translate into the CGR.By convention negative fluxes indicate an uptake by the Earth surface.

Figure 1 .Figure 2 .
Figure 1.Schematic explanation of the uncertainty propagation.Each spatiotemporal estimate of NCE is computed as the sum of randomly selected estimates of the 9 fluxes contributing to NCE (see Eq. 1, here denoted by F i ).For this study we compute 200 estimates of NCE.Uncertainties can be assessed at different spatial scales by first aggregating all NCE estimates to the desired scale and then using the 200 members for uncertainty estimation.5

Figure 3 .
Figure 3. Gridded spatial patterns of NCE.a) 2001-2010 decadal mean.b) Uncertainty; 1SD across the NCE ensemble.c) Relative uncertainty; uncertainty normalized by absolute mean.Latitudinal plots in b) and c) denote median across latitudes.

Figure 5 .
Figure 5. NCE (2001-2010 decadal mean) in RECCAP regions over land, including (red) and without fossil fuels (blue).Shown are median, interquartile range (box) and 1.5 x interquartile range (whiskers).The regional estimates including uncertainties of NCE collected in Ciais et al (in revision) are underlain in grey.NA: North America, SA: South America, EU: Europe, AF: Africa, RU: Russia, EA: East Asia, SAs: South Asia, SEA: South East Asia, AU: 5 Australia, Rest: remaining land areas.

Figure 6 .
Figure 6.NCE in RECCAP regions over the ocean.Shown are median, interquartile range (box) and 1.5 x interquartile range (whiskers) of the Marine fluxes.The RECCAP estimates including uncertainties are underlain in grey.NP: North Pacific extratropics, TP: Tropical Pacific, SP: South Pacific extratropics, AR: Arctic Ocean, NS: Northern Subtropics, EQ: Equatorial, SS: Southern Subtropics, NI: Northern Indian Ocean, SI: Southern Indian 5 Ocean, SO: Southern Ocean.

Figure 7 .
Figure 7.Comparison of NCE-FF with NCE from inversions (by construction without FF) on interannual time scales.Both time series were zero-centered by adding an offset of 13.23 PgC / year for NCE-FF and 3.74 PgC / year for NCE from inversions.Lines show mean, shading is 1 SD.

Table 2 .
(2001fluxes are also summarized graphically in Figure2(mean over2001-  2010).Our best surface-data driven bottom-up global estimate of NCE is -5.4±2.0PgC/ year.That means, that the observation-based datasets suggests a large net sink, even if FF and E LUC are included in NCE.By contrast, the accurately measured CO 2 growth rate constrains NCE being a net CO 2 source to the atmosphere 4.3±0.1 PgC / year(2001( -2010, Le    , Le Andreae et al. (2015)4)imates without FF over the same regions (Figure5).In regions without tropical forest except NA (that is, EU, RU, EA, SAs, and AU) the estimates byCiais et al (in revision)are within the interquartile range of our assessment.For NA and regions containing the tropics, our approach shows a much stronger C sink.Given thatCiais et al. (revision)rely on an independent method, this demonstrates that a relatively good understanding and observational coverage of net C fluxes exists for EU, RU, AU, and EA to some extent.It is somewhat surprising that both approaches largely differ over North America, where good observational coverage for instance through eddy-covariance towers exist.The comparison also reveals the high uncertainties and biases in bottom-up estimates of NCE over tropical forests (seeSect.3.2, but also Gloor et al. (2012)andValentini et al. (2014)) and underlines the importance of long-term ground based measurement campaigns (e.g.RAINFOR, http://www.rainfor.org/,Malhietal.(2002),andATTO,Andreae et al. (2015),Zhou et al. ( Haverd et al., 2013)Gloor et al., 2012), Europe(EU, Luyssaert et al., 2012), Africa(AF, Valentini et al., 2014), Russia (RU,Dolman et al., 2012), East Asia (EA,Piao et al., 2012), South Asia (SAs,Patra et    al., 2013), and Australia (AU,Haverd et al., 2013).No regional study is yet available for Southeast Asia (SEA).Greenland, Middle East, Ukraine, Kazakhstan and New Zealand are omitted in regional RECCAP studies because of the difficulty to obtain local ground-based observations.Ciais et al. (in revision) collected the regional estimates and combined them with estimates of lateral transport to estimate C budgets for each region.NEE inCiais et al. (in revision)minus C export by rivers should in principal