Patterns and controls of inter-annual variability in the terrestrial carbon budget

The terrestrial carbon fluxes show the largest variability among the components of the global carbon cycle and drive most of the temporal variations in the growth rate of atmospheric CO2. Understanding the environmental controls and trends of the terrestrial carbon budget is therefore essential to predict the future trajectories of the CO2 airborne fraction and atmospheric concentrations. In the present work, patterns and controls of the inter-annual variability (IAV) of carbon net ecosystem exchange (NEE) have been analysed using three different data streams: ecosystem-level observations from the FLUXNET database (La Thuile and 2015 releases), the MPI-MTE (model tree ensemble) bottom–up product resulting from the global upscaling of site-level fluxes, and the Jena CarboScope Inversion, a top–down estimate of surface fluxes obtained from observed CO2 concentrations and an atmospheric transport model. Consistencies and discrepancies in the temporal and spatial patterns and in the climatic and physiological controls of IAV were investigated between the three data sources. Results show that the global average of IAV at FLUXNET sites, quantified as the standard deviation of annual NEE, peaks in arid ecosystems and amounts to∼ 120 gCm−2 y−1, almost 6 times more than the values calculated from the two global products (15 and 20 gCm−2 y−1 for MPI-MTE and the Jena Inversion, respectively). Most of the temporal variability observed in the last three decades of the MPI-MTE and Jena Inversion products is due to yearly anomalies, whereas the temporal trends explain only about 15 and 20 % of the variability, respectively. Both at the site level and on a global scale, the IAV of NEE is driven by the gross primary productivity and in particular by the cumulative carbon flux during the months when land acts as a sink. Altogether these results offer a broad view on the magnitude, spatial patterns and environmental drivers of IAV from a variety of data sources that can be instrumental to improve our understanding of the terrestrial carbon budget and to validate the predictions of land surface models.

either of the two fluxes can cause large variations in their difference.
It has been long debated whether it is GPP or TER that controls the spatial and temporal variability in NEE. Several studies have ascribed inter-annual variability in NEE to variability in either GPP (Ahlström et al., 2015;Janssens et al., 2001;Jung et al., 2011Jung et al., , 2017Stoy et al., 2009;Urbanski et al., 2007) or TER (Morgenstern et al., 2004;Valentini et al., 2000) or both (Ma et al., 2007;Wohlfahrt et al., 2008b). GPP and TER show comparable ranges of IAV, typically larger in absolute terms than that observed for NEE due to the temporal correlation between the two gross fluxes (Richardson et al., 2007). Given that photosynthesis and respiration may respond differently to environmental drivers (Luyssaert et al., 2007;Polley et al., 2008), the interpretation of climate impacts on the variability in NEE requires an understanding of the relation between the IAV of NEE and that of GPP and TER (Polley et al., 2010).
The environmental factors driving the IAV of NEE (IAV NEE ) include climate, physiology, phenology, and natural and anthropogenic disturbances (Marcolla et al., 2011;Richardson et al., 2007;Shao et al., 2015). Understanding the spatio-temporal variability in NEE and its controlling mechanisms is essential to assess the vulnerability of the terrestrial carbon budget, to evaluate the land mitigation potentials and to quantify the ecosystem capacity to store carbon under future climatic conditions (Heimann and Reichstein, 2008). Besides, quantifying inter-annual variability in NEE is a prerequisite for detecting longer-term trends or step changes in flux magnitude in response to climatic or anthropogenic influences and identifying its drivers (Cox et al., 2000;Lombardozzi et al., 2014).
Despite the broad literature on the subject, very few examples of IAV analysis based on multiple data streams are available in the literature (Desai et al., 2010;Pacala, 2001;Poulter et al., 2014). In the present study, patterns and controls of the inter-annual variability in NEE have been analysed using three different data streams: ecosystem-level data from the FLUXNET database, the MPI-MTE (model tree ensemble) bottom-up product resulting from the statistical upscaling of in situ flux data (le Maire et al., 2010) and the Jena CarboScope Inversion top-down product, which estimates land (and ocean) fluxes from atmospheric CO 2 concentration measurements and atmospheric transport modelling (Rödenbeck et al., 2003). In particular, this analysis aims to (i) assess the magnitude and the spatial pattern of the IAV of NEE (IAV NEE ), (ii) investigate the role of key climatic variables, like temperature and precipitation, in driving the spatial pattern of IAV, and (iii) identify the role of photosynthesis and respiration as sources of IAV NEE . Finally, the consistencies and discrepancies among the different data products are analysed and critically evaluated.
2 Materials and methods

Datasets
Data on an ecosystem scale were retrieved from two releases of the FLUXNET dataset, namely La Thuile (http://fluxnet. fluxdata.org/data/la-thuile-dataset/) and 2015 (http://fluxnet. fluxdata.org/data/fluxnet2015-dataset/). These datasets contain half-hourly data of carbon dioxide, water vapour and energy fluxes that are harmonized, standardized and gap-filled. Time series of NEE and of the component fluxes GPP and TER, together with air temperature and precipitation, were used in the present analysis. Flux data have the advantage of representing direct observations of in situ IAV; however, at most sites the time series are still too short for a proper analysis of the temporal variability in NEE (Shao et al., 2015). For this reason only sites with a minimum of 5 years of observations and an open data distribution policy were selected. A subset of 89 sites satisfied the two criteria, among which there were 27 evergreen needle-leaf forests, 5 evergreen broadleaf forests, 12 deciduous broadleaf forests, 6 mixed forests, 12 grasslands, 8 croplands, 6 sites with closed and open shrublands, 7 wetlands, and 6 sites with savannahs and woody savannahs.
On a global scale, two sources of gridded data were used: a bottom-up data product, namely the MPI-MTE product (Jung et al., 2009), and, as a top-down product, the Jena CarboScope CO 2 Inversion (Rödenbeck et al., 2003). The MPI-MTE dataset is built with a machine learning technique to upscale in space and time the flux observations from the global network of eddy covariance sites (FLUXNET) integrated with climate and remote-sensing data for the time period 1982-2011 (Jung et al., 2009). Global maps for GPP and TER at 0.5 • spatial resolution and monthly temporal resolution were used, while NEE fields were calculated as the difference between the gross fluxes. This product has become a reference dataset to evaluate process-oriented land models and remote-sensing estimates of primary productivity, despite the uneven distribution of eddy covariance sites on which it is trained. It integrates a large amount of in situ measurements and remote-sensing and meteorological observations using a machine learning technique and has been proved to reproduce spatial patterns and seasonal variability in fluxes well (Jung et al., 2009). On the other hand, the prod-uct has some shortcomings: for instance, the effects of land management, land use change and CO 2 fertilization are not represented. The MPI-MTE has been recognized as underestimating the inter-annual variability in carbon fluxes. These limits may be due to the missing representation of some key determinants of IAV like changes in soil and biomass pools, disturbances (e.g. fires), ecosystem age, management activities and land use history. Finally, the lag between external forcing and ecosystem response is not represented in the product (Jung et al., 2011).
To derive surface fluxes, the Jena CarboScope Inversion combines modelled atmospheric transport with highprecision measurements of atmospheric CO 2 concentrations. Atmospheric transport is simulated by a global threedimensional transport model driven by meteorological data. For consistency with the MPI-MTE product, monthly averaged NEE land fluxes from the s81_v3.6 version of the product were used here, at a spatial resolution of 5 • × 3.75 • . The Jena Inversion is particularly suited for the analysis of temporal trends and variability since it is based on a temporally constant observation network (14 atmospheric stations for the version s81_v3.6). Weaknesses of the Jena Inversion product are linked (i) to the sparse density and biased spatial distribution of the sampling network, whose geometry affects the flux estimates in a systematic way, (ii) to data gaps, (iii) to incompleteness of the accounted fluxes, since the calculation is based on CO 2 data only, while atmospheric carbon comes also from CO and volatile organic compounds, and (iv) to potential systematic errors in the transport model (Baker et al., 2006).
As the inversion estimates the total land flux, calculated as the difference between the total surface flux and prescribed anthropogenic emissions, it also includes CO 2 emissions from fires in addition to NEE. For improving the consistency with the other two datasets (MPI-MTE and FLUXNET) that do not account for fire emissions, we subtracted fire emissions from the inversion estimates using a harmonized combination of the products RETRO (Schultz et al., 2008) for the period 1982-1996and GFED4 (van der Werf et al., 2010 for the period 1997-2013. RETRO is a global gridded dataset (at 0.5 • spatial resolution) for anthropogenic and vegetation fire emissions of several trace gases, covering the period from 1960 to 2000 with monthly time resolution. GFED4 combines satellite information on fire activity and vegetation productivity to estimate gridded monthly fire emissions at a spatial resolution of 0.25 • since 1997. RETRO and GFED4 were harmonized using the overlapping years 1997-2000 to calculate calibration coefficients as the ratio of GFED4 to RETRO for latitudinal bands of 30 • . The RETRO time series was then multiplied by these coefficients and the resulting time series of fire emissions was finally subtracted from the land flux of the Jena Inversion. It is worth noting that the remaining flux from the inversion is the sum of land use change emissions and NEE while the MPI-MTE does not account for the land use change flux.
In order to analyse the role of climatic drivers in the interannual variability, global maps of temperature and precipitation were used. Gridded air temperatures were obtained from the Climatic Research Unit (CRU) at the University of East Anglia on a monthly timescale and 0.5 • × 0.5 • spatial resolution, based on an archive of monthly mean temperatures provided by more than 4000 weather stations (Jones et al., 2012). Precipitation fields were obtained from the Global Precipitation Climatology Centre (GPCC) product at 0.5 • and monthly time steps (Schneider et al., 2014). This product is based on a large dataset of monthly precipitation from more than 85 000 stations and is provided by NOAA ESRL Physical Sciences Division (PSD; Boulder, Colorado, USA). The MODIS MCD12C1 land cover product (Friedl and Brodley, 1997) was used to classify the land pixels and to calculate statistics by plant functional type. MCD12C1 provides the dominant land cover types at a spatial resolution of 0.05 • using a supervised classification algorithm that is calibrated using a database of land cover training sites. Product resolutions were harmonized using the aggregate function of the raster R package (Hijmans and van Etten, 2014).

IAV analysis
The inter-annual variability in NEE was estimated as the standard deviation of annual NEE values generated by trend and residuals, computed on time windows of 12 months shifted with a monthly time step (Luyssaert et al., 2007;Shao et al., 2015;Yuan et al., 2009) and calculated with the same methodology for the three data streams used in the analysis. Average values of IAV for plant functional type (PFT) were determined using the PFT classification of FLUXNET sites and the MCD12C1 product (aggregated at the appropriate spatial resolution using the dominant PFT) for the MPI-MTE and the Jena Inversion. Map grid cells were also classified according to mean annual temperature and precipitation, and the mean value of IAV NEE and normalized IAV NEE were calculated for each climate bin. The IAV of both MPI-MTE and the Jena Inversion was normalized by the average GPP of the specific climate bin from the MPI-MTE.
For the two gridded products, which provide a 30-yearlong time series , the IAV was partitioned into two components, namely the variance explained by the temporal trend and that due to annual anomalies (Ahlström et al., 2015). For this purpose a linear model was fitted to the time series at each pixel, and the determination coefficient of the regression was used to measure the fraction of variance explained by the trend, whereas its complement to 1 was the fraction of variance due to anomalies.
The spatial correlation between IAV and climatic drivers (air temperature and precipitation) was analysed on a global scale for the MPI-MTE by calculating the spatial correlation coefficient between the temporal standard deviation (IAV amplitude) of NEE and the average annual temperature or precipitation in moving spatial windows of 15 • × 11.5 • (which means 31 × 21 pixels for MPI-MTE). The latitudinal averages of these correlation coefficients were calculated for latitudinal bands of 30 • . This analysis was not replicated on the Jena Inversion because on a fine scale the spatial variability in the fluxes in this product is mainly controlled by the prior estimates. In fact, the optimization algorithm of the inversion spatially allocates the fluxes proportionally to the prior; hence, grid cells with higher productivity will change more if compared to cells with a lower prior. Finally, in order to identify which process between photosynthesis and respiration drives IAV NEE , for FLUXNET and MPI-MTE linear regressions between NEE and GPP and NEE and TER were fitted for each site/pixel and the difference between the determination coefficients of the two linear regressions was computed. Since GPP and TER cannot be derived from inversion products, we performed a similar analysis using NEE of the carbon uptake period (CUP, sum of negative monthly NEE) and of the carbon release period (CRP, sum of positive monthly NEE) as proxies of GPP and TER for all the three data streams. Also in this case NEE was linearly regressed with NEE CUP , and NEE CRP to detect which of the two processes drives the variability in NEE. IAV NEE and IAV controls were also analysed in a climatic space defined by mean annual temperature and precipitation. Finally, annual anomalies of the two global products used in the present analysis were compared with the estimates derived from the Global Carbon Project (GCP) (Le Quéré et al., 2016).

IAV patterns
The spatial pattern of inter-annual variability for the three datasets is shown in Fig. 1. The IAV of NEE at individual FLUXNET sites ranges from 15 to 400 gC m −2 y −1 with an average of 130 gC m −2 y −1 . On average the most northern sites show a lower temporal variability both in Europe and in North America (Fig. 1a). A global map of IAV NEE is shown also for MPI-MTE (Fig. 1b) and the Jena Inversion ( Fig. 1c) at the original spatial resolutions of the two products. The observed range of IAV is similar for the two gridded products and substantially lower than that observed at the site level, probably due to the spatial averaging of the land fluxes that dampens the temporal variability. The mean global value of IAV is in fact 15 and 20 gC m −2 y −1 for MPI-MTE and the Jena Inversion, respectively, and hence about 1/6 of the sitelevel IAV. The two gridded products confirm the decreasing trend of IAV toward northern latitudes observed at flux sites. A general decrease in IAV NEE at higher latitude for both evergreen needleleaf forests (ENF) and deciduous broadleaf forests (DBF) was also observed by Yuan et al. (2009) although for none of the two PFTs were these trends significant. In terms of IAV, the two global products show a reasonable qualitative correspondence for North America and Eurasia, whereas they disagree for South America, with MPI-MTE showing a minimum of IAV in the humid tropics, where the inversion product shows large variability. MPI-MTE in particular shows maximum values along the eastern coast of South America while the Jena Inversion shows an almost opposite pattern. A similar behaviour is also observed in Africa, where the top-down product shows a maximum in central Africa, while MPI-MTE shows a minimum in the Congo Basin and higher values in arid zones like the Sahel and southern Africa. These discrepancies could, on the one hand, be ascribed to the limits of the bottom-up approach in dealing with the low seasonality of the fraction of absorbed radiation (FaPAR) in evergreen broadleaf forests, given the relevance of this predictor in the MPI-MTE estimates. A second reason for the discrepancy could be the CO 2 emissions from land use change that are particular relevant in some tropical areas but are not accounted in the MPI-MTE estimates. On the other hand, the fine-scale estimates of the inversion are largely determined by the a priori weighting pattern, which has been chosen proportional to time-mean net primary production (NPP; from the LPJ model) as a vegetation proxy (Rödenbeck et al., 2003). As the atmospheric data can only constrain larger-scale patterns comparable to the distances between the stations, this means that IAV will be locally higher where mean NPP is high and vice versa.
As far as the Northern Hemisphere is concerned, a good correspondence is observed in western Eurasia, while some discrepancies are observed in other zones; for example MPI-MTE shows a large IAV in India, probably driven by the changes in FaPAR related to agricultural intensification, which does not emerge from the inversion product that has little observational constraint in this area. To summarize, the spatial pattern of IAV in the two products better agrees in the Northern Hemisphere for temperate and cold temperate zones, whereas for the Southern Hemisphere, and in particular for the humid evergreen forests, the two products show a poor match. In general, it has to be considered that both the MPI-MTE product and the Jena Inversion are driven by data from surface networks that are very limited in the tropics and the Southern Hemisphere, and, therefore, these observationdriven estimates are under-constrained in those areas. These results highlight that for achieving more robust and consistent estimates of the terrestrial carbon fluxes, it is of key importance to increase the availability of atmospheric and ecosystem flux observations in the tropical region, either by establishing new sites where the network is sparse or improving the sharing of data where the monitoring stations are available but not connected to global networks (e.g. flux stations in Amazonia).
The results presented in the maps of Fig. 1 are summarized in the climate space in Fig. 2. The left panels show that peak values of IAV are located in different climate regions for the two gridded products (temperate humid for MPI-MTE and tropical humid for the Jena Inversion). These results highlight that top-down and bottom-up estimates do not agree on the main sources of temporal variability in the terrestrial carbon budget and call for more investigation to pin down the reasons for these large discrepancies. Given that the IAV of NEE increases with the primary productivity at the FLUXNET sites (Fig. 3), in Fig. 2 (right column) we normalized the IAV of both MPI-MTE and the Jena Inversion by the average GPP of the specific climate bin from the MPI-MTE. Normalization using GPP (which is always positive) offers a more robust metric of relative IAV if compared to normalization with NEE (that spans 0). Figure 2 shows either the mean IAV (left column) or the ratio of the mean IAV and GPP (right column) in each climate bin, since this metric is less sensitive to outliers than the mean of ratios and gives more weight to points with larger fluxes. The normalized IAV shows a consistent pattern between the three differ- q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0 −100 −300 −500 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q ent data products, with a clear decreasing trend with increasing temperature and precipitation (i.e. increasing productivity). Ultimately arid regions seems to have the higher relative variation in land carbon fluxes, in accordance with previous findings (Ahlström et al., 2015). Interestingly, the two gridded products show slightly different climatic location for the peak in relative IAV, with MPI-MTE pointing to warm arid regions, whereas the Jena Inversion points to cold arid systems.
The dependency of IAV NEE on GPP and on NEE CUP is shown in Fig. 3 for the three datasets. Both for FLUXNET We think that this mismatch is due to the prominent role that FaPAR has in the MTE approach. In fact, canopy greenness is particularly stable in the tropical humid forests, generating this unusual pattern of low relative IAV in regions of high productivity. These contrasting results for key regions like the Amazon and the Congo Basin confirm the large uncertainty of the IAV estimates in areas with limited observational constraints. In these regions, climate sensitivities derived from estimates of the inter-annual variability in the terrestrial carbon budget therefore have to be carefully interpreted (Fang et al., 2017). The importance of the spatial scale of analysis for the IAV NEE has been explored both for FLUXNET sites and the global products (i.e. MPI-MTE and the Jena Inversion) (Fig. 4). The two global products show good agreement at the native Inversion resolution (5 • × 3.75 • ) and at a global level, when only one global value is retrieved by spatially averaging all the pixels of the original maps. For the MPI-MTE product, the observed IAV decreases regularly at decreasing map resolution. By contrast, the Jena Inversion shows a rapid descent followed by a stabilization. FLUXNET sites and their aggregation at increasing distance also show a decreasing IAV with higher values compared to the global prod-ucts. The slope of the lines in Fig. 4 reveals the degree of spatial compensation between anomalies (steeper slopes are generated by stronger compensation and therefore lower spatial coherence), which leads to a decrease in IAV NEE at the increase in the spatial extent of the observations. Of the three products, MPI-MTE shows the gentler slope and therefore the larger spatial coherence of the anomalies. This is possibly due to the missing representation of land disturbances in the MTE methodology, which may ultimately lead to an overestimation of the spatial coherence in the land CO 2 flux anomalies.
The fractions of IAV NEE generated either by temporal trends or by annual residuals are summarized in Fig. 5 for the two global gridded products. For MPI-MTE, more than 80 % of the IAV is explained by residuals at all latitudes. Only in limited zones like the Congo Basin and western Amazonia does MPI-MTE show a relative minimum in the importance of residuals, but this global product might underestimate the total variability in these zones (see Fig. 1b). Residuals explain the largest share (between 62 and 90 %, average 77 %) of the temporal variability also in the Jena Inversion, with a higher relevance of trends in the Southern Hemisphere. The inversion product shows several hotspots of trend-driven variability, like southern Africa, South America and northern Eurasia, which is indeed reported as an area of increasing productivity in the last decades (Forkel et al., 2016). In the interpretation of these results it is important to consider that MPI-MTE is generated by the statistical upscaling of FLUXNET data, using climate and FaPAR as predictors. This methodology relies on the assumption of a constant ecosystem response to climate drivers, and for this reason the product cannot reproduce the influence of some environmental factors (e.g. increasing CO 2 concentration or nitrogen deposition) that may alter these responses and that are not reflected in input variables like climate or FaPAR. By contrast, inversion products do not make any assumption on the climate dependence of ecosystem functioning but also includes emission from land management and land use change that may hide or emphasize the NEE trends. In summary, it is important to note that, despite the important climate trends, in the last 30 years the temporal variability in the land carbon balance has been driven by annual residuals, confirming the dominant role of climate variability in the terrestrial C budget (Le Quéré et al., 2014).
For the two gridded products the analysis of IAV (either in terms of absolute IAV NEE or normalized with NEE CUP ) was disaggregated by plant functional type (Fig. 6). The analysis in terms of absolute IAV NEE shows that savannahs and woody savannahs (WSAV-SAV) are the PFTs characterized by the larger IAV and variability within the PFT. This was found both for the MPI-MTE and the Jena Inversion product and confirms the results of a recent study (Ahlström et al., 2015) in which semi-arid ecosystems were found to account for the largest fraction (39 %) of the global IAV in net biome productivity. This variability was found to be signif- icantly related to the length of the growing season (Ma et al., 2007) and is driven by the uncertainty in water supply in arid systems. In terms of normalized IAV the two gridded products show different behaviours, shrublands (CSH-OSH) being the most variable PFT for MPI-MTE, while the inversion data show a higher variability for evergreen broadleaf forests (EBF) and WSAV-SAV. As observed on a pixel scale in Fig. 1, even at PFT level the results obtained from FLUXNET sites show a higher variability than the gridded products. In general, at FLUXNET sites IAV is proportional to ecosystem productivity ( Fig. 4) with the maximum values observed in EBF, deciduous broadleaf and mixed forests (DBF-MF), and grassland-croplands (GRA-CRO) and the minimum in wetlands (WET). The large value of IAV observed in GRA-CRO is presumably also affected by the potential large impact of management in these ecosystems that can either reduce (e.g. by irrigation) or increase the climate-induced variability (e.g. by changing crops or fertilization schemes). In general, the disaggregation of IAV NEE by PFTs shows rather similar results between the two gridded products, both in terms of magnitude and distribution. The largest difference is observed in the evergreen broadleaf forests whose absolute and relative variability is much larger in the inversion, possibly as a result of the intensive disturbances that have occurred in these ecosystems over the last decades and that are not captured with the MTE methodology.

Climate dependence of IAV
The climatic dependence of the spatial variability in IAV NEE on a global scale for the MPI-MTE product (Fig. 7) shows a clear pattern, with positive correlations in temperaturelimited areas at northern latitudes and a negative temperature dependence in water-limited zones (Braswell et al., 1997). These observations agree with Reichstein et al. (2007), who report that GPP shifts from soil water content to air temperature dependency at around 52 • N. These opposite temperature dependences will probably lead to future contrasting changes in IAV. In fact, under a global warming scenario, the northern latitudes will be characterized by a larger sink (Zhao and Running, 2010) but also by a larger temporal variability, while arid zones like the Mediterranean Basin, central eastern Australia and sub-Saharan Africa will probably experience a reduction in IAV linked to large-scale droughts and consequent reduction in primary productivity (Ciais et al., 2005). Concerning precipitation, the MPI-MTE product shows more complex spatial patterns, with negative correlation in the humid tropics, temperate Europe and the southeast USA and positive correlation elsewhere. It is worth noting that this analysis does not account for the potential lag between climate events and ecosystem responses, in line with the assumption adopted in the formulation of the MPI-MTE statistical model (Jung et al., 2011). Potential delays in the ecosystem responses to precipitation anomalies as observed from field studies (Doughty et al., 2015) may eventually explain the contrasting spatial pattern shown by this data product.
The climate dependencies of IAV are further separated between the variability due to trends and anomalies (Fig. 7b, d). The two components of IAV NEE mostly show an agreement in the sign of the climatic controls, meaning that the environmental drivers have the same effects on trends and anomalies and therefore on the long and short timescales. This is a relevant finding because it supports the use of IAV to investigate medium-term climatic responses. In general, anomalies show a higher correlation than trends, probably due to the larger magnitude of the variance attributed to this component. In conclusion, the spatial patterns shown in the maps of Fig. 7 and the agreement between the two components of IAV shown in the bar plots indicate that the temperature controls of the IAV of NEE are in general the same as for the primary productivity (i.e. positive in colder biomes and negative in warmer regions), while the contrasting results observed for precipitation suggest that the role in the spatial and temporal variability played by water availability is unclear, proba-bly because of the temporal correlation between precipitation and temperature anomalies, as shows by Jung et al. (2017). The analysis of the climate drivers of IAV was not performed for the Jena Inversion because for this product local variation in IAV is heavily driven by the prior estimates of NPP and therefore results have limited sensitivity to the atmospheric constraints.

Physiological drivers of IAV
An improvement in the mechanistic understanding of IAV NEE can be achieved by partitioning the net flux into its two components: GPP and TER. Partitioned fluxes are available for FLUXNET sites and for derived products like MPI-MTE, while they cannot be derived from atmospheric inversions. For this latter product the fluxes during the CUP (NEE < 0) and during the CRP (NEE > 0) were used in this analysis as proxies of GPP and TER, respectively.
To investigate how good these proxies are, the ratios TER / GPP during CUP and GPP / TER during CRP were analysed at FLUXNET sites and for each pixel of the MPI-MTE product and averaged by PFT (Fig. 8a). As far as the MPI-MTE product is concerned, TER ranges from 55 to 78 % of GPP during the CUP, while GPP is 56 to 80 % of TER during the CRP; hence, on average about two-thirds of the signal comes from GPP (TER) in the CUP (CRP). These ratios show a certain variability among PFTs, with ENF having the larger imbalance between the two fluxes and the lowest TER / GPP ratio during CUP (due to the strong seasonality of GPP in this PFT), while the two fluxes are not so well partitioned for EBF (ratio ∼ 0.8), which are characterized by a long growing season with consistently large fluxes of GPP and TER. The other PFTs show an average ratio value of ∼ 0.65 both in CUP and CRP. In summary, it can be inferred that NEE during CUP is dominated by the GPP signal, while NEE during CRP is dominated by TER albeit to a smaller extent as is shown by the frequency distributions in Fig. 8bc calculated from the MPI-MTE product. The distribution of the TER / GPP ratio during the CUP is in fact narrower and peaks at a value of 0.7, while a broader distribution is observed for the GPP / TER ratio during the CRP. As expected, there is a larger spread in the composition of NEE during CRP across the world, and this is linked to the larger variability in the seasonality of GPP that may actually go to 0 in the dormancy season, while TER is always positive.
In order to identify which of the gross fluxes controls the variability in the net land flux, we assessed the fraction of variance (R 2 ) of NEE explained by GPP or TER (for MPI-MTE and FLUXNET) and CUP or CRP (for all products). Results given in Fig. 9 show the difference of the determination coefficients between the two regressions (NEE versus GPP or TER; NEE versus NEE CUP or NEE CRP ) and are used to determine which component dominates the inter-annual variation in NEE. Blue zones in Fig. 9 are regions where IAV NEE is driven by photosynthesis (GPP or CUP), the dif- ference R 2 GPP −R 2 TER (or R 2 CUP −R 2 CRP ) being positive, while in the red zones IAV NEE is mainly controlled by respiration (TER or CRP). Figure 9a shows that, in most of the land area, the IAV NEE is driven by GPP both at FLUXNET sites and for the MPI-MTE product. The same data products show an even clearer dominance of NEE CUP over IAV (Fig. 9b). The Jena Inversion product shows that, although most of the globe is NEE CUP driven, there are quite a few areas that are weakly CRP driven, like the eastern US, arid regions in Africa and the Amazon Basin, probably because these areas are estimated to be CO 2 sources in this inversion and therefore NEE is dominated by NEE CRP (data not shown). When latitudinal profiles are considered, all the products show that GPP and NEE CUP control the temporal variability in yearly NEE more than TER or NEE CRP (le Maire et al., 2010). Results shown in the global maps of Fig. 9 are represented in the climatic space in Fig. 10. Map pixels were classified according to mean annual temperature and precipitation. For each climate bin the difference between the determination coefficients for NEE versus GPP and TER is shown. Across the whole climate space, the IAV retrieved from the MPI-MTE product is mostly controlled by CUP and GPP, although the difference in R 2 in the case of GPP and TER is low. The Jena Inversion by contrast shows climate areas where IAV is CRP driven, especially in intermediate to high temperature classes. Simi- lar results have been reported across several PFTs by Yuan et al. (2009) and Ahlström et al. (2015) using FLUXNET site data and MTE products. A higher correlation of IAV with GPP rather than with TER in deciduous forests has also been reported by Barr et al. (2002) and Wu et al. (2012). These results suggest that ecosystem fluxes during the CUP, and in particular photosynthesis more than respiration, consistently control the inter-annual variability in NEE on all the spatial scales for both bottom-up and top-down data products (Janssens et al., 2001;Luyssaert et al., 2007;le Maire et al., 2010;Urbanski et al., 2007;Wohlfahrt et al., 2008a;Wu et al., 2012;Yuan et al., 2009). These results highlight that temporal variations in photosynthesis and in ecosystem CO 2 exchange during the carbon uptake period therefore drive the short-term climate sensitivity of the global carbon cycle consistently across different regions and climates. The possibility of interpreting these short-term responses as long-term potential impacts of climate change is therefore to be disputed, given the limited role that respiration appears to play in modulating the rapid reactions of the terrestrial biosphere to environmental drivers.
Finally, in order to place our analysis in a broader context, global annual values of the gridded products used in the present analysis have been compared with the estimates of the Global Carbon Project (Fig. 11). On an annual timescale, the Jena Inversion shows excellent agreement with the GCP, and this is not surprising since GCP land fluxes are estimated as a residual term from the atmospheric CO 2 budget and are  therefore not completely independent of the Jena Inversion product. On the other hand, this analysis highlights how the MTE bottom-up approach is barely correlated (R 2 = 0.015, p = 0.52) with the top-down estimates, both in terms of trend and of anomalies. These discrepancies may be partially explained by the missing representation of land disturbances (land use change, land management) in the MTE product.

Conclusions
In conclusion, this study assessed the temporal variability in the terrestrial C budget with three different datasets to diagnose common patterns and emerging features. Some discrepancies between data products have emerged, in particular in the tropics, where a chronic deficiency of atmospheric and ecosystem observations severely limits the accuracy of large-scale assessments. On the other hand, several important global features have been identified and confirmed by the different products, like (i) the dominant role played by photosynthesis in the short-term variability in the land carbon budget, (ii) the high relative IAV in water-limited ecosystems, and (iii) the dependence of IAV on spatial scales and ecosystem productivity. Ultimately, the variability in the land fluxes observed in the recent decades proved to be extremely valuable to investigate the controlling mechanisms and the sensitivity and vulnerability of the terrestrial C balance to climate drivers.
Edited by: Alexey V. Eliseev Reviewed by: Nir Krakauer and two anonymous referees