Reviews and syntheses: An empirical spatiotemporal description of the global surface–atmosphere carbon ﬂuxes: 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 charac-terization of C cycling requires detailed information on spatiotemporal patterns of surface–atmosphere ﬂuxes. 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 syn-thesize a wide range of observation-based spatially explicit surface–atmosphere CO 2 ﬂuxes from 2001 to 2010, to identify the state of today’s observational opportunities and data limitations. The considered ﬂuxes include net exchange of open oceans, continental shelves, estuaries, rivers, and lakes, as well as CO 2 ﬂuxes related to net ecosystem productivity, ﬁre emissions, loss of tropical aboveground C, harvested wood and crops, as well as fossil fuel and cement emissions. Spatially explicit CO 2 ﬂuxes are obtained through geostatis-tical and/or remote-sensing-based upscaling, thereby mini-mizing 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 this study is to identify key uncertainties and observational shortcomings that need to be prioritized in the expansion of in situ observatories. Uncertainties for NCE and its components are derived using resampling. In many regions, our NCE estimates agree well with independent estimates from other sources such as process-based models and atmospheric inversions. This holds for Europe (mean ± 1 SD: 0.8 ± 0.1 PgC yr − 1 , positive numbers are sources to the atmosphere), Russia (0.1 ± 0.4


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 humancaused perturbation of the C cycle is largely driven by emissions of CO 2 from burning fossil fuel C previously stored 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, during the years 2000-2009 about 45 % of the anthropogenic CO 2 emissions each year stay in the atmosphere, the rest being taken up by the oceans (27 %) and land (27 %) (Le Quéré et al., 2015).
Recently, a case has been made for a globally policyrelevant 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 • C (UNFCCC, 2015).In such a policy-relevant C observing system, an uncertainty assessment for each data stream of 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 outgassing 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 surface-atmosphere 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 net 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 surfaceatmosphere 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.Section 4 addresses the benefits and limits of the current observational system for constraining global net CO 2 fluxes.Section 5 summarizes the main conclusions drawn from this synthesis.

Data and methods
We collected ensembles of data-driven estimates of the NCE for the major subsystems of the Earth from 2001 to 2010 (Table 1).Each data set was aggregated to 1 × 1 • spatial resolution.All data sets have an original spatial resolution of at least 1 × 1 • and were aggregated to the lower-resolution common grid.The temporal resolution is monthly.For data sets that were only available at yearly timescale or once over the complete time period (Table 1), we distributed fluxes evenly across all months.In this synthesis, we include NCE 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 (NCE) derived from the combination of all the above products with the growth rate of atmospheric CO 2 (CGR).Combining all fluxes, the overall NCE between the Earth's surface and the atmosphere is given as (1) www.biogeosciences.net/14/3685/2017/Biogeosciences, 14, 3685-3703, 2017

Uncertainty estimation and propagation
For five of the nine variables contributing to the NCE of Eq. (1), we obtained multiple spatiotemporally explicit estimates directly from the raw data products.These estimates are expected to sample the uncertainty in each variable.We could obtain 10 different estimates for the terms Marine and Crops respectively, 50 estimates for Rivers, eight estimates for NEP, two estimates for E LUC , and only one estimate for each of the remaining four variables (Table 1).To estimate NCE including spatiotemporally explicit uncertainties, we randomly draw ensemble members from each of these nine terms to create an integrated ensemble for NCE (Eq.1).
With the available estimates, we could in principle create 10×50×8×10×2 = 80 000 spatiotemporal explicit estimates of NCE.From these 80 000 possible NCE estimates we randomly select 200 (to reduce computational expense) to construct the NCE ensemble, which is used in the remainder of the paper.This resampling approach is illustrated in Fig. 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 1 standard deviation (SD) of the annual mean across all realizations), interannual variability (IAV, defined here as 1 SD of annual means across all available years) and the coefficient of variation (CV = IAV/mean) for each of the nine 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 method's degrees of freedom for internal validation.Repeating this relationship building process and withholding different sets of validation data has created the five ensemble members used for this study.For the Jena CarboScope mixed-layer scheme, we used the five sensitivity cases with varied parameters such as changes in correlation length 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 airsea pCO 2 difference ( pCO 2 ) (Jähne et al., 1987) and a quadratic dependence of the wind speed at a height of 10 m (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 m 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 CO 2 product.

Shelves
For continental shelf seas we derived the pCO 2 from 3 × 10 6 surface pCO2 measurements extracted from the SO-CATv2 database (Bakker et al., 2014) and observational atmospheric pCO 2 data (GLOBALVIEW-CO2, 2012).The local CO 2 air-sea flux values were then obtained using a winddependent quadratic formulation parameterized as in Wanninkhof et al. (2013) and wind speeds extracted from a crosscalibrated 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 timescale depending on the data coverage.These temporally and regionally averaged air-sea CO 2 fluxes were then disaggregated using a 1 • 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 four 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 airwater 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 • resolution from the length of the coastline and a type-specific length to estuarine surface ratio (Dürr et al., 2011).

Marine
We created 10 marine estimates by combining the 10 estimates for the open oceans with 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 a 3 × 3 mean window.This same filter is also applied to the rest of the merged data set 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
Fifty 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 • (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., 2012(Raymond et al., , 2013) ) 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 simula-tion 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 Database (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 • 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 eight 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 ("FLUXCOM-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 offsite (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 to 2011 at 0.05 • 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 • , added all uncertainty estimates within one 1 • pixel and used them as estimates for 1 standard deviation on the new 1 • grid.Assuming Gaussian-distributed errors we sampled 1000 values at each pixel and used 10 maps of the 5th, 15th, . . ., 95th quantiles as different ensemble members.Data were 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 one estimate of globally gridded forest harvesting data around year 2000 as described in the Supplement 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 Supplement (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 1 × 1 • resolution.
We used E LUC only in those pixels where the average of the two estimates (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 simulation indicated an uncertainty of about 20 % (1SD) for continentalscale 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 area 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 × 0.1 • ) 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_Fuel.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 UN-FCCC.

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 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 of NCE, we used annual means of 10 atmospheric CO 2 inversions collected in Peylin et al. (2013) 2001-2010, Le Quéré et al., 2015).Thus, there is a large mismatch with our NCE of 9.7 ± 2.0 PgC yr −1 .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 (Schwalm et al., 2015), 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 (Fig. 3a).According to these estimates, tropical land areas are a larger CO 2 sink than Biogeosciences, 14, 3685-3703, 2017 www.biogeosciences.net/14/3685/2017/ 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 the 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; Fig. 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 (Fig. 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 undersampled, 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 homogeneous 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 (Fig. 3c).Only in the Southern Ocean is the relative uncertainty 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 (Fig. 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 (Fig. 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 (Fig. 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 (black line in Fig. 4d shows NCE-FF).FF contribute little 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 yr −1 , which is unrealistic, if compared to, for example, forest biomass inventories (Pan et al., 2011).Without this large sink, global NCE would be of −0.4 ± 1.8 PgC yr −1 .This corrected estimate (assuming neutral C exchange in tropical forests) is still a sink more than 4 PgC yr −1 larger than the global NCE accurately constrained by CO 2 growth rate observations (4.3 ± 0.1 PgC yr −1 ).Including missing fluxes (e.g.biogenic fluxes and emissions from wetlands; see Sect.4.2) for which we do not have spatially explicit estimates (see Sect. 4.4) could close this gap.These considerations suggest that the CO 2 sink of tropical forests from FLUXCOM is probably strongly overestimated and responsible for at least half of the global mismatch with the observed CO 2 growth rate (see Sect. 4.1 for further discussion).2013) and Lauerwald et al. (2015).Thus, these estimates are not fully independent of those presented in this study because they use the same river fluxes, fire emissions, and FF emission.All other fluxes are independent.The RECCAP budgets were based on inventories, and in some instances on process models results.For E LUC , RECCAP publications used regional data sets or bookkeeping models (a bookkeeping model combines local or regional observationbased estimates of C stocks before land use change and trajectories of C stocks after land use change including slash, biomass, soil carbon, and harvested wood products with changing areas of different land use types, Houghton and Nassikas, 2017), that are independent from estimates gathered in Sect.2.3.4.The RECCAP regions include North America (NA; King et al., 2015), South America (SA; 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, the Middle East, Ukraine, Kazakhstan and New Zealand are omitted in regional RECCAP studies because of the difficulty of obtaining local ground-based observations.Ciais et al. ( 2017) collected the regional estimates and combined them with estimates of lateral transport to estimate C budgets for each region.NEE in Ciais et al. (2017) minus C export by rivers should in principle be equal to our NCE estimates without FF over the same regions (Fig. 5).In regions without tropical forest except NA (that is, EU, RU, EA, SAs, and AU) the estimates by Ciais et al. (2017) are within the uncertainty range of our assessment.For NA and regions containing the tropics, our approach shows a much stronger C sink.
Given that Ciais et al. (2017) rely on an independent method, this demonstrates that a relatively good understanding and observational coverage of net C fluxes exists for EU, RU, AU, SAs, 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 (see Sect. 3.2, but also Gloor et al., 2012;Valentini et al., 2014) and underlines the importance of long-term ground-based measurement campaigns in those regions (e.g.RAINFOR, http://www.rainfor.org/:Malhi et al., 2002;ATTO: Andreae et al., 2015;Zhou et al., 2014).

RECCAP over ocean
We compare our estimates of marine 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).While there is large overlap between the pCO 2 data used in the RECCAP estimates and our NCE estimates over oceans, different independent methods have been used to obtain flux estimates of ocean CO 2 exchange (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 process-based 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 (Fig. 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 CO 2 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 (Fig. 7).
The comparison shows 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 −1 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 (Fig. 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 www.biogeosciences.net/14/3685/2017/Biogeosciences, 14, 3685-3703, 2017 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.
4 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 surface-atmosphere 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 observationdriven 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 data set 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 data sets (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 res-piration 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 (Fig. 4).
Given the difference between NCE and inversions in the tropics (Fig. 4), we can assume that a bias of FLUXCOM NEP towards a too high C 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-thanaverage NEP are possibly over-represented in FLUXNET.However, such an age dependency (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 the land surface and the 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 if known): emissions from biogenic volatile organic compounds (VOCs) amount to approximately 0.76 PgC yr All known missing fluxes add up to an additional C release of about 4 PgC yr −1 .Although substantial, they do not cover the mismatch of more than 9 PgC yr −1 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 yr −1 , Sect.3.2).This significant amount of missing fluxes prohibits constraining FLUX-COM 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 and 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 data streams used (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 model-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 data sets (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 providing not only 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, of 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 www.biogeosciences.net/14/3685/2017/Biogeosciences, 14, 3685-3703, 2017 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 were 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 yr −1 .
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.
Data availability.All variables used in Eq. ( 1) (including all ensembles runs listed in Table 1) are available for download (2001-2010, monthly, 1 • ).For NCE we provide the 200 ensemble runs which are used in this study.The landing page of the DOI is a data portal.After registration any user can download the data without restriction (https://dx.doi.org/10.17871/GEOCARBON_synth_obs_v2).
Competing interests.The authors declare that they have no conflict of interest.

Figure 1 .
Figure 1.Schematic explanation of the uncertainty propagation.Each spatiotemporal estimate of NCE is computed as the sum of randomly selected estimates of the nine 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.

Figure 2 .
Figure 2. Different components of observation-driven C exchange between the Earth's surface and the atmosphere.Red arrows denote a flux from the surface to the atmosphere (net source), green arrows denote a flux from the atmosphere to the surface (net sink).Units are in PgC yr −1 .

Table 1 .
Data sets used in this study including reference, time period and number of ensemble runs.If not specified, temporal resolution is monthly.
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.Data scarcity precludes including all 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.
, available at the same spatial and temporal resolution.Atmospheric CO 2 inversions estimate the spatiotemporal explicit CO 2 exchange between the Earth surface and atmosphere using atmospheric CO 2 measurements and a transport model.

Table 2 .
Net carbon exchange for different subsystems and variables that contribute to NCE (Eq.1; negative numbers are surface uptake).Uncertainty (Unc.) is SD over ensemble runs.IAV is SD over annual values(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)and CV is coefficient of variation, computed as IAV divided by the mean.PgC yr −1 .That means that the observation-based data sets suggest 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 to being a net CO 2 source to the atmosphere of 4.3 ± 0.1 PgC yr −1 (