Global covariation of carbon turnover times with climate in terrestrial ecosystems

A global, observation-based assessment of whole-ecosystem carbon turnover times shows that the overall mean global carbon turnover time is about 23 years and that locally its spatial variability depends on precipitation at least as strongly as on temperature. Feedback between the terrestrial carbon cycle and climate is partly determined by changes in the residence time of carbon in land ecosystems. Using observation-based gross primary production estimates, remote-sensing based estimates of vegetation biomass and new estimates of total soil organic carbon to full depth, Nuno Carvalhais et al. calculate a spatially explicit estimate of global mean residence times of carbon in land ecosystems. They arrive at an overall mean global carbon turnover time of about 23 years with carbon residing, on average, in the vegetation and soil near the Equator for a shorter time period than at northern latitudes. The paper reports expected dependencies of carbon turnover time with temperature, but also identifies a strong association with precipitation, implying that future carbon cycle climate feedbacks may depend more strongly on changes in the water cycle than currently expected based on Earth system model studies. The response of the terrestrial carbon cycle to climate change is among the largest uncertainties affecting future climate change projections1,2. The feedback between the terrestrial carbon cycle and climate is partly determined by changes in the turnover time of carbon in land ecosystems, which in turn is an ecosystem property that emerges from the interplay between climate, soil and vegetation type3,4,5,6. Here we present a global, spatially explicit and observation-based assessment of whole-ecosystem carbon turnover times that combines new estimates of vegetation and soil organic carbon stocks and fluxes. We find that the overall mean global carbon turnover time is years (95 per cent confidence interval). On average, carbon resides in the vegetation and soil near the Equator for a shorter time than at latitudes north of 75° north (mean turnover times of 15 and 255 years, respectively). We identify a clear dependence of the turnover time on temperature, as expected from our present understanding of temperature controls on ecosystem dynamics. Surprisingly, our analysis also reveals a similarly strong association between turnover time and precipitation. Moreover, we find that the ecosystem carbon turnover times simulated by state-of-the-art coupled climate/carbon-cycle models vary widely and that numerical simulations, on average, tend to underestimate the global carbon turnover time by 36 per cent. The models show stronger spatial relationships with temperature than do observation-based estimates, but generally do not reproduce the strong relationships with precipitation and predict faster carbon turnover in many semi-arid regions. Our findings suggest that future climate/carbon-cycle feedbacks may depend more strongly on changes in the hydrological cycle than is expected at present and is considered in Earth system models.

The response of the terrestrial carbon cycle to climate change is among the largest uncertainties affecting future climate change projections 1,2 . The feedback between the terrestrial carbon cycle and climate is partly determined by changes in the turnover time of carbon in land ecosystems, which in turn is an ecosystem property that emerges from the interplay between climate, soil and vegetation type [3][4][5][6] . Here we present a global, spatially explicit and observation-based assessment of whole-ecosystem carbon turnover times that combines new estimates of vegetation and soil organic carbon stocks and fluxes. We find that the overall mean global carbon turnover time is 23 z7 {4 years (95 per cent confidence interval). On average, carbon resides in the vegetation and soil near the Equator for a shorter time than at latitudes north of 756 north (mean turnover times of 15 and 255 years, respectively). We identify a clear dependence of the turnover time on temperature, as expected from our present understanding of temperature controls on ecosystem dynamics. Surprisingly, our analysis also reveals a similarly strong association between turnover time and precipitation. Moreover, we find that the ecosystem carbon turnover times simulated by state-of-the-art coupled climate/carbon-cycle models vary widely and that numerical simulations, on average, tend to underestimate the global carbon turnover time by 36 per cent. The models show stronger spatial relationships with temperature than do observation-based estimates, but generally do not reproduce the strong relationships with precipitation and predict faster carbon turnover in many semiarid regions. Our findings suggest that future climate/carbon-cycle feedbacks may depend more strongly on changes in the hydrological cycle than is expected at present and is considered in Earth system models.
The largest global gross exchanges of carbon occur at the interface between the atmosphere and the terrestrial biosphere 7 . Changes in the net exchange of CO 2 between the land and the atmosphere may provide positive or negative feedbacks to increasing atmospheric CO 2 concentrations and, thus, changes in climate 1,8 . The response of the net exchange of CO 2 to climate depends on the response of carbon uptake (gross primary production (GPP)) and on how carbon residence times change simultaneously. It is thus of importance to quantify the time that carbon resides in terrestrial ecosystems and its spatial covariation with climate. Furthermore, global modelling studies show a stronger convergence of GPP estimates in comparison with wider ranges in whole ecosystem carbon stocks 9 , reflecting the strong spread in the modelled residence times of carbon.
In steady state, that is, when the exchange of carbon between two reservoirs is balanced, the turnover time (t, yr) equals the mean residence time 10 . Assuming a balance between assimilation and losses of ecosystem carbon, t can be calculated from the total reservoir size (C total , kgC m 22 ) and the influx or the outflux (kgC m 22 yr 21 ) as 10 t 5 C total /flux (1) For terrestrial ecosystems, the total reservoir size equals the carbon stocks in vegetation and soils. The influx is the carbon uptake through gross primary production GPP and the outflux includes all carbon losses (terrestrial ecosystem respiration, fire emissions, lateral export and so on). We relax the strict steady-state assumption, calling t the apparent whole ecosystem turnover time, computed as the ratio of C total to GPP, and interpret the quantity as an emergent diagnostic at the ecosystem level (Methods section on turnover times). We note that the turnover time, or mean residence time, of carbon in ecosystems emerges from the turnover of compartments varying greatly in their individual turnover times 5,11 (for example leaves, wood, different soil organic carbon fractions). Hence, carbon allocation, leaf, root and wood turnover, plant mortality and soil carbon decomposition are key processes that regulate the terrestrial turnover times, and which are controlled by climate variability 12 .
Here we combine and enhance recently derived estimates of the carbon stocks in vegetation and soil to obtain global total terrestrial ecosystem carbon stocks and their observation-based uncertainties at 0.5u resolution (Methods section on global carbon estimates). We merge remotesensing-based carbon stock estimates for tropical and Northern Hemisphere vegetation (Methods section on vegetation carbon) with enhanced soil organic carbon estimates based on the Harmonized World Soil Database (HWSD) and a dedicated circumpolar soil organic carbon map (which better accounts for carbon in permafrost-affected high-latitude soils; Methods section on soil organic carbon). Total soil organic carbon stocks are estimated down to full depth (that is, beyond the 100 cm depths often reported, but see Methods section on soil organic carbon). Total carbon stocks (soils and vegetation) amount to 2,807 z855 {555 Pg of carbon C (PgC) and are predominantly within tropical forests, which contain almost 25% of the total global stocks, followed by boreal forests (18%) (Methods section on global carbon estimates). Per unit area, the total carbon stocks vary largely between and within biomes (Fig. 1a), where tropical forests and northern high latitudes exhibit the highest stocks. Substantial uncertainties are located in tundra regions (interquartile range over the mean, ,38%) and in tropical savannahs and grasslands (,30%) (Extended Data Fig. 1 and Extended Data Table 1).
Using these total ecosystem stocks and the observation-based GPP estimates in equation (1), we derive a global t of 23 yr (ranging between 18 yr (percentile 2.5) and 29 yr (percentile 97.5)). We note that this duration is an estimation of the mean residence time of a carbon atom in terrestrial ecosystems from its initial fixation by photosynthesis until its respiratory (including autotrophic respiration) or non-respiratory loss. A previous collection of global estimates of net primary production, carbon in soils and carbon in vegetation allows an approximate mean estimate of t of 21 6 7 yr, assuming 50% of autotrophic respiration costs 13 . Still, the spatial variation of t observed in our analysis spans almost two orders of magnitude, that is, between ,7 yr (first percentile) and 439 yr (99th percentile) ( Fig. 1b and Methods section on mean carbon turnover times).
The longest turnover times are found in cold biomes (tundra, 65 z13 {20 yr; boreal forests, 53 z20 {8 yr), in temperate grassland and shrubland (41 z13 {9 yr) and in desert regions (36 z14 {9 yr), whereas tropical forests and savannah exhibit the shortest turnover times (14 z4 {3 and 16 z6 {4 yr, respectively) ( Table 1). Using localized regression analysis (Methods section on correlation analysis), we find that, spatially, t covaries significantly (P , 0.05) with mean annual temperature or with total annual precipitation in 86% of the globe ( Fig. 1c and Supplementary Information Section 2). There is, however, a strong variability in the spatial correlations between t and temperature and t and precipitation ( Fig. 1c and Extended Data Fig. 2). Negative correlations between temperature and t are widely observed, and can be linked to the expected decomposition responses to temperature [14][15][16] . However, significant positive correlations emerge in regions of forest/herbaceous cover transitions (or patchiness) and in warm arid environments ( Supplementary Information Section 4), where precipitation shows the strongest correlations with t, indicating that moisture effects may dominate and override temperature effects. No clear dominant patterns are observed in tropical forests, suggesting that climate has a limited effect on the spatial variability of t there and that nutrient availability 17 or natural and human disturbances 18,19 , or both, have greater effects. Globally, we observe a higher frequency of stronger spatial correlations between t and precipitation (in ,55% of land grid cells) than between t and temperature (,45%) (Extended Data Figs 3 and 4).
Turnover times vary considerably with latitude, ranging from 255 yr (mean t above 75u N) in the high northern latitudes to 15 yr in the equatorial tropics (Fig. 2a). We find that the most rapid latitudinal changes exist between the sub-Arctic zones and the temperate zones, and near the tropical circles (between 20 and 40u N). Within the tropics (between 20u N and 20u S), the variations in t are comparatively minor. In the Northern Hemisphere, in the transition zone between 50 and 65u N, the spatial covariation of t is strongest with temperature, but south of 50u N precipitation is the dominant associated variable (until the Equator and below 40u S (Fig. 2c, d, in blue)). Across all latitudes, higher precipitation is associated with shorter residence times (negative partial correlations), whereas correlations between t and temperature are more variable across latitudes, being low in the northern tropical zone and carbon, turnover times of carbon in terrestrial ecosystems and its spatial covariation with climate variables. Global distribution of estimated total ecosystem carbon (C total ) density in each grid cell (kgC m 22 ), and total ecosystem carbon mass per biome (PgC)) (a), median of turnover times of carbon (t) in terrestrial ecosystems (b), and the spatial covariation (r, Pearson correlation coefficient) of t with climate variables (c). The insets in a and b show the C total and t estimates per biome according to previous classifications (Extended Data Table 1); the uncertainty bars per biome report the 95% confidence intervals (CI, ranging between percentiles 2.5 and 97.5). The ranges in C total span 4 kgC m 22 (approximately the 5th percentile) to 254 kgC m 22 (maximum estimate of mean carbon stocks), and the colour code in the t map is binned between t # 4 yr and 439 yr (the 99th percentile). Spatial correlations (moving window of 5.5u by 5.5u) with temperature (tas) and precipitation (pr) are shown only for confidence levels above 95%. The missing regions in southern Australia and New Zealand are due to missing data in the vegetation carbon data set (Methods section on total vegetation carbon). Total ecosystem turnover times of carbon per biome, estimated using equation (1) and stocks and fluxes aggregated per biome. Data estimates of t were aggregated by biomes defined previously (Extended Data Table 1), and the ranges reported are the 2.5th (P 2.5) and 97.5th (P 97.5) percentiles from the ensemble of t estimates, which can be interpreted as the confidence intervals in these estimates. The total represents the global t including all biomes.

RESEARCH LETTER
also north of 60u N. Overall, across the latitudinal range, the spatial correlations of t with temperature and precipitation are significant, but only low to moderately so. We note that carbon turnover in ecosystems will depend on the time-integrated effect of climate variables and summary statistics, and that mean total precipitation and mean temperature can serve only as simple proxies. The significance of soil organic carbon stocks in explaining the spatial variability of t is pervasive (Extended Data Fig. 5b). For approximately 80% of the land surface, t covaries more strongly with soil carbon stocks than with vegetation carbon stocks. However, the residence times of carbon in terrestrial ecosystems should tend to increase with vegetation longevity, and with allocation towards woody biomass. Hence, we expected extensive positive correlations of t with tree cover. Nevertheless, we also find negative correlations (Extended Data Fig. 5a). Precipitation, which is associated with tree cover, could overshadow the tree effect by increasing turnover times disproportionally, but the negative correlations persist even when we control for precipitation (Extended Data Fig. 6 and Supplementary Information Section 3). An increasing probability of fire related to increasing fuel loads with above-ground biomass 20,21 , thereby reducing turnover times, is a possible explanation for this apparent paradox. Others are the contribution of trees to wetter microclimates 22,23 and increasing nutrient availability 23,24 in regions of low tree density and transition regions. Additionally, other factors like natural and anthropogenic disturbances 25 or management activities 26 can accelerate rates of turnover and, consequently, reduce mean residence times.
We calculated the turnover times of carbon (equation (1) in models from the Coupled Model Intercomparison Project Phase 5 (CMIP5)) (see Methods section on CMIP5 and Supplementary Information Section 5). The broad latitudinal patterns of t in the CMIP5 ensemble ( Fig. 2b) are consistent with the observations (Fig. 2a) (Pearson correlation coefficient, r 5 0.88; P , 0.0001), but with a mean underestimation bias in the latitudinal profile of 47% (normalized average error). However, the zonal mean carbon turnover times vary by a factor of 2 to 40 across the analysed CMIP5 models (Fig. 2b). The models differ strongly with respect to correlations of t with their (modelled) climate variables 27 . Across almost all latitudes, we find a range from positive to negative correlations for both temperature and precipitation. In the Northern Hemisphere, modelderived correlations with temperature tend to be more negative than the observation-derived correlations 27 (Fig. 2c). For latitudes below 50u N, the model-derived correlations with precipitation are more positive when compared with the observation-derived correlations (Fig. 2d). In the tropical zones, the CMIP5 ensemble predicts increasing turnover times associated with higher precipitation, which contrasts with the observationderived estimates.
Overall, the CMIP5 models correlate with the observation-derived estimates of t (r 2 5 0.38, P , 0.001), but exhibit shorter turnover times (Fig. 3); this is reflected also in the global turnover times (,36% lower; Extended Data Table 2). The bias is particularly pronounced in the high northern latitudes and in the seasonally dry biomes of northern tropical Africa, North America, central southern Asia and east Australia. An underestimation of turnover times in the high latitudes can potentially be ecosystem turnover times of carbon and associations to temperature and precipitation from data and models. Evaluation of latitudinal patterns and climate association of t for data and models. The latitudinal gradients in t from data (a) and from models (b) show distinctive associations with temperature (tas, c) and precipitation (pr, d). For consistency, the temperature and precipitation data sets considered for the model analysis are also model outputs (Methods section on CMIP5). The comparisons are based on partial correlations, controlling for precipitation when evaluating the association of t with temperature (and vice versa), and are performed at the spatial scale of the NorESM1-M model output, to minimize artefacts in the correlations caused by differences in spatial resolution.

LETTER RESEARCH
explained by either neglecting, or having an incomplete representation of, permafrost processes 28,29 , although this can only partly explain differences in biases between North America and Siberia. Globally, in 31% of the land grid (35% of the global land area analysed), fewer than onequarter of the models are within the confidence intervals of the data. Even assuming a 50% error in the observed carbon stocks and, consequently, in turnover times, this would not explain modelled turnover times differing by more than a factor of two from observations. Furthermore, despite representing soil organic carbon pools with long residence times, CMIP5 models do not provide an explicit representation of soil organic carbon vertical profiles. This could partly explain the observed differences but not the systematic underestimation of soil organic carbon up to 1 m depth in northern latitudes (Extended Data Fig. 7), which warrants attention to the representation of soil organic pools and vertical profiles in models. Biases in simulated climate may also lend significant biases to model estimates of t, although a comparison of differences in t and spatial covariations with climate reveals that models fall short in describing the climate responses seen in the observation-derived data (Supplementary Information Section 6). Other possible reasons for this pronounced model bias may include responses to, or biases in, modelled soil moisture 27 or insufficient sensitivity of decomposition to drought. Moreover, adaptation of vegetation to dry conditions includes leaf sclerophylly, long leaf lifespans and higher wood densities in shrubs, which together lead to increased turnover times. In addition, interactions with nutrient cycles (for example that of nitrogen) may slow the turnover of carbon in ways which are not represented in models 30 . The spatial analysis shown here does not imply that the relationships between t and climate factors are the same in the temporal dimension: these relationships emerge from the effects of climate-and other factors-through time, to which models should be comparable. In this regard, the emergence of appropriate model-data integration frameworks is essential for a consistent transfer of information from observation-based estimates of t to modelling approaches 31 .
We have presented an observation-based estimate of the total terrestrial carbon pool size and whole-ecosystem carbon turnover times and its spatial variation at 0.5u resolution with associated uncertainties. Our findings suggest significant hydrological control of carbon turnover, probably as relevant as temperature, adding to the well-known coupling between carbon and water cycling for photosynthesis, and calling for a better understanding of changes with the hydrological cycle. Although the ensemble mean of state-of-the-art coupled climate/carbon-cycle models reproduces the temperature-driven latitudinal patterns of carbon turnover times, we note an important underestimation bias and differences between models of more than one order of magnitude. The pronounced underestimation of whole-ecosystem carbon turnover times in semi-arid regions calls for a more accurate description of hydrological processes and water-carbon interactions. We expect that improved representations of the adaptation of vegetation to water availability, fire dynamics, and physicochemical and microbial soil organic carbon stabilization mechanisms 32 , in addition to permafrost dynamics, will probably help to address the aforementioned biases. Overall, these results emphasize the role of water on the carbon dynamics in the terrestrial biosphere and suggest that future climate/carbon-cycle feedbacks will be more sensitive to changes in the water cycle than expected and represented in state-of-the-art models.  turnover times between models (t mod ) and observation-derived (t obs ) data ensembles. Stippling indicates locations where fewer than onequarter of the models are within the 5th and 95th percentiles of the data. Here the deserts are filtered out according to a Köppen-Geiger classification (Extended Data Table 2) and to a minimum GPP of 10 gC m 22 yr 21 . The missing regions in southern Australia and New Zealand are due to missing data in the vegetation carbon data set (see Methods). Author Information Reprints and permissions information is available at www.nature.com/reprints. The authors declare no competing financial interests.

RESEARCH LETTER
Readers are welcome to comment on the online version of the paper.

METHODS
Estimates of total soil organic carbon based on global databases. The Harmonized World Soil Database 33 condenses a comprehensive collection of geographic information on soil physical and chemical properties from regional and national inventories all over the world. The HWSD is organized in mapping units, each consisting of particular combinations of different soils referred to as 'soil ID' from here on. For every soil ID, among other variables, the database reports texture, bulk density and concentration of organic carbon for the top (0 to 10 or 30 cm) and subsoil (from 30 cm to 1 m depth) layers. Estimates of total organic carbon for each soil ID can then be computed per layer as follows: Here soil organic carbon stocks (SOC, kg C m 22 ) is estimated from organic carbon content (OC, wt%), layer thickness (D, m), gravel content (G, vol%) and bulk density (BD, kg m 23 ). Such an approach allows for estimates of SOC in the top layer (0-30 cm) and in the subsoil layer (30-100 cm). We used these two estimates to fit two empirical models of cumulative SOC (equations (2) and (3)), which were then integrated until the full soil depth (D f ) per soil ID was reached, as follows: SOC~a log Dbz1 ð Þ ð3Þ Here K, I, a and b are empirical parameters estimated per soil ID. Estimates using equation (2)  we included an alternative model formulation that has shown a strong fit to data, as well as a faster saturation of cumulative SOC with depth ( Supplementary Fig. 1). Here the latter model represents a more conservative estimation of full-depth SOC. The full depth of the soil (D f ) was extracted from the Global Soil Texture And Derived Water-Holding Capacities database 36 . This database contains standardized values of soil depth and textures for the globe, which were selected for the same soil types in the same continents, according to continents defined in ref. 37.
The HWSD was still a work in progress at the time of our study, such that data from certain regions in the world still needed updating and were therefore considered less reliable. Two such regions were North America and northern Eurasia 33 . Therefore, as an alternative to the HWSD, we also considered the Northern Circumpolar Soil Carbon Database 38,39 (NCSCD) in our SOC estimates for northern latitudes 38,39 . We generated a set of global SOC estimates, which included factorial combinations of SOC from the HWSD, extrapolated to full soil depth using both empirical models, and also used the NCSCD data set for northern latitudes. These data sets were aggregated from ,1 km 2 (0.01u by 0.01u) to ,55 km 2 (0.5u by 0.5u) resolution. Our global estimate of total soil organic carbon was 2,397 z860 {561 PgC (mean; upper limit, percentile 97.5; lower limit, percentile 2.5) (Supplementary Fig. 2). The global soil carbon stocks are comparable to a previous estimate of total soil organic carbon of 2,344 PgC in the top 3 m of soil 34 (Supplementary Table 1). The range of estimated SOC values varied significantly between the different biome types across the world (Supplementary Table 1). According to our estimates, tropical biomes (forest, savannahs, grasslands) together account for 32% of the global soil carbon stock, and the areas of largest integrated stock are found in tropical (20%) and boreal (19%) forests.
A comprehensive assessment of the uncertainties in the SOC estimates should integrate the uncertainties from several sources: from (1) the uncertainties in the soil profiles and depth information to (2) uncertainties in the spatial extrapolation to a global extent, including (3) uncertainties in the extrapolation to full depth and (4) those emerging from the different data sources considered here (the HWSD and the NCSCD). The ability to quantify the uncertainties stemming from all these sources, and propagate them to the final SOC estimates, is limited by the available information. But here we are able to explicitly propagate the uncertainties that stem from the methods used to extrapolate SOC to full depth (from both the aforementioned empirical models; equations (2) and (3)) and from the different data sources (HWSD and NCSDC), by creating individual SOC estimates, which are combined individually with the estimates of vegetation stocks and GPP for explicit propagation of the uncertainties in t. The quantification of the uncertainties in the total SOC stemming from depth are propagated by exploring the variability present in the soil depth from the WISE data set for the same soil types as in the HWSD. We do so by contrasting soil depth standard deviations against soil depth (s Df 5 0.19D f ; r 5 0.28, P , 10 210 , N 5 4,790), for which we draw additionally 50 random samples of depth (with a s Df of 19%) to estimate the uncertainties in total SOC that may stem from the soil depth considered.
Overall, this approach is based on the best available information, and addresses uncertainties by evaluating the spread in the generated ensembles. We acknowledge that the integration of uncertainties stemming from the soil profiles and depth, and from the regionalization to the global scale, could alter the uncertainties presented here. The provision of information on these sources of uncertainty, and the ability to tackle them in future estimates is essential to a more comprehensive assessment of data uncertainty. Deriving total vegetation carbon. Our global estimates of total vegetation carbon were derived from a collection of estimates for pan-tropical regions 40 and for northern and temperate forests 41 based on radar remote-sensing retrievals 42 . Above-and below-ground biomass uncertainty for the tropical regions was propagated from errors in measurements, allometric relations, sampling and predictions 40 . In the Northern Hemisphere, estimates accounted for uncertainties in the BIOMASAR GSV data, wood density data and biomass compartment data 41 . On regional scales, the Northern Hemisphere biomass map in comparison with inventory-based data showed strong agreement (Russia: r 2 5 0.78 and NRMSE 5 0.35; United States: r 2 5 0.90 and NRMSE 5 0.32; Europe: r 2 5 0.70 and NRMSE 5 0.40; NRMSE denotes the root mean squared error divided by the mean of the observations) and can thus be considered a very suitable product at 0.5u resolution. Evaluation results for the United States and Europe have shown that this data set might slightly underestimate high carbon densities due to the use of C-band radar data, but there was no systematic error detected in the intercomparison in Russia 41 .
One shortcoming of the above two products is the sole consideration of tree forms in their estimates. Therefore, to account for the herbaceous biomass in our estimates, we assumed a mean turnover time of one year in the live vegetation fraction per grid cell, and, given that the costs of autotrophic respiration vary significantly 43 , we took the respiratory costs (a) to lie in the range 25%-75% (uniformly distributed): Here C H is the herbaceous component of carbon in vegetation; GPP is the gross primary production, based on the newest data driven estimates 44 ; a is the respiration cost; and f H is the fraction of each 0.5u grid cell considered as herbaceous in the SYNMAP 45 . The correlation (r 2 ) between the data set accounting for non-woody stocks and the original vegetation stock estimates was 0.98 6 0.015, and the normalized mean absolute error was 0.07. By accounting for herbaceous cover in our global C stock estimates, we obtained differences of less than 1% in the mean global t, with the highest differences observed for croplands (1.3%), temperate grasslands and shrublands (0.6%), and wetlands (0.6%). Therefore, accounting for herbaceous plants in our global carbon stocks did not make much of a difference to the final estimates of t. We note that these differences may not fully reflect the dynamics in natural vegetation types, which may include below-ground perennial roots or rhizomes. However, even a threefold increase in vegetation mean residence times of carbon would result in a difference of less than 2% in total ecosystem turnover times globally. These results reflect the large contribution of woody vegetation and, mostly, of soil organic carbon to the global carbon stock estimates. Overall, terrestrial vegetation holds about 442 6 146 PgC (mean 6 standard deviation), which is ,16% of the global organic carbon estimated on land (Supplementary Table 2). Excluding herbaceous vegetation, our estimate of 429 6 144 PgC in forests encompasses the latest rounded estimate of 300 PgC derived from global inventory data 46 . The most significant part of vegetation carbon is found in tropical forests and tropical savannahs and grasslands (62%), and temperate and boreal forests, and temperate grasslands and shrublands incorporate circa 25%. Global carbon estimates. The ensemble of vegetation carbon pools was composed of 200 members, assuming normally distributed uncertainties in the satellite-derived C stocks and in GPP estimates, and a uniform distribution of autotrophic respiration costs (see previous Methods section). Each member of the ensemble of total soil carbon was individually added to each of the members of vegetation carbon ensemble. We randomly sampled 200 members of this ensemble to achieve the final data ensemble of total ecosystem carbon stocks (Fig. 1a) with uncertainties (Extended Data Fig. 1). The uncertainty bands in the latitudinal profiles in Extended Data Fig. 7 report the 5th and 95th percentiles of the data ensemble per latitudinal window. Global mean turnover times of carbon. Turnover time is commonly defined as the ratio between the total size of a reservoir and its outflux 47 . For terrestrial ecosystems, the total reservoir size is equal to the carbon stock in vegetation and soils, and the outflux comprises all carbon losses (respiration of autotrophic plants, respiration of heterotrophic organisms, losses by fire and harvest). Under the assumption that the ecosystem is neither gaining nor loosing carbon (steady state), the turnover time can equivalently be calculated as the ratio between the carbon stock in vegetation (C veg ) and soils (C soil ), and the flux into this reservoir, GPP: In steady state, t is the average time that newly assimilated carbon spends in terrestrial ecosystems before it is respired, burnt or harvested. Acknowledging that our definition of t (equation (4)) hinges on the steady-state assumption, we call t the RESEARCH LETTER apparent whole-ecosystem turnover time and interpret the quantity as an emergent diagnostic at ecosystem level ( Supplementary Information Section 1). An ensemble of apparent whole-ecosystem turnover times (t) was obtained by applying equation (1) to a random permutation of the mean annual GPP  and total ecosystem carbon data sets (N 5 200). The resulting uncertainties had wide ranges in space (Supplementary Fig. 3) and between biomes ( Table 1). The uncertainty bands in the latitudinal profiles in Fig. 2a report the 5th and 95th percentiles of data ensemble per latitudinal window. Benchmarking our current results against other observation-based estimates of global carbon turnover time is hampered by the fact that previous studies have mostly focused on soils. Global turnover times of carbon in soils were reported to range from 27 yr (ref. 48) to 32 yr (ref. 49), and were generally considered to lie between 30 and 40 yr, assuming strong variations for different ecosystems 50 . The spatial variations in the residence times of carbon in soils also reflect climate controls, exhibiting longer residence times in cold biomes at high latitudes (as shown in ref. 51 using 14 C for forest soils). Such spatial variations are also seen in t (Fig. 1a), which shows a latitudinal range from 255 yr in high northern latitudes (mean t north of 75u N) to 15 yr in the equatorial tropics (Fig. 2a). In ref. 13, several estimates of global NPP (table 13), soil organic carbon (table 14) and vegetation stocks (table 12) were reported. Assuming respiration costs of around 50%, whole-ecosystem turnover time can be estimated as t 5 C total /(NPP 1 RA) (RA, autotrophic respiration). By combining the different quantities of stocks and fluxes, a global t of 21 6 7 yr (mean 6 standard deviation) can be estimated.
For model evaluation, the separation between the soil and vegetation components of t is useful for a more detailed diagnostic of model performance. We define t soil as the ratio between C soil and net primary production (NPP) and t veg as the ratio between C veg and GPP. However, the spatial representation of observationbased empirical estimates of NPP is still hampered by the difficulty of accounting for autotrophic respiration fluxes. Hence, the confidence level for NPP spatial estimates (for example, from ref. 52) still falls short to the levels of GPP estimates 44 . An approach is to consider an ensemble of NPP fields supported by observationbased empirical estimates based on climate patterns. A set of state-of-the-art global NPP fields from refs 52-57 is used to build a data-driven empirical ensemble. Each of these NPP fields is combined individually with the C soil ensemble members to build an ensemble of t soil for comparison to the CMIP5 model results. For the construction of t veg ensemble, we relied on the GPP and vegetation carbon stocks described above.
The data used here can be obtained from http://www.bgc-jena.mpg.de/geodb/ BGI/tau.php. Climate data. Climate data are based on the European Centre for Medium-Range Weather Forecasts (ECMWF) ERA-Interim reanalysis product 58 and have been bias-corrected as described in ref. 59. We obtained daily data with a grid cell size of 0.5u during 1979-2010 from ECMWF. This reanalysis data was bias-corrected against the WATCH forcing data 60 using an overlapping period of 1979-2001 and following the standard procedure 61 . Ideally, this approach conserves statistical moments of the distribution, for example the mean and variance. In addition, the number of rainy days remains unchanged. The WATCH forcing data serves as the reference data set because it has already been bias-corrected against other climatic data sets 60 . For the analysis, we computed the mean fields of each of the climate variables by averaging the whole data sets per grid cell between 1982 and 2005. Correlation analysis. The association between t and climate is assessed locally for each grid cell of the global data sets by computing local correlations using a 5.5uby-5.5u moving window (11 by 11 grid cells). This approach enables the assessment of the local importance of climate factors (Supplementary Information Section 2) and yields a global estimate of regional covariation of t with climate variablestemperature and precipitation. The association between t and climate is determined by analysing the Pearson correlation, partial correlation coefficients and the nonparametric Spearman's rank correlation coefficient. To disentangle the relative importance of temperature and precipitation in determining the spatial patterns of t, we (1) used the Lindeman-Merenda-Gold (LMG) method 62 , implemented in ref. 63, which allows to quantify the contribution of different correlated regressors (here temperature and precipitation) to a multiple linear regression model; and (2) quantified the changes in the residual sums of squares by removing each independent variable from a bivariate regression between temperature and precipitation: RI v( RSS v {RSS tas,pr )=TSS, where RI v is the relative importance of variable v, RSS v is the residual sum of squares of the regression of t with variable v, RSS tas,pr is the residual sum of squares of the regression of t against temperature and precipitation, and TSS is the total sum of squares. Finally, the metrics (RI tas and RI pr ) are normalized (divided by r 2 ) to sum to 1 (ref. 63). In addition, we performed a conditional independence test on rejecting the null hypothesis that t is independent of precipitation or temperature given the dependence on temperature or, respectively, precipitation 64 . Latitudinal correlations between t and climate variables are based on partial correlation coefficients between t and temperature or precipitation, controlling for precipitation or, respectively, temperature. The analysis is conducted on a common grid size for the CMIP5 models (see next section) and the observation-derived t and climate for a latitudinal window of ,9.5u (5 grid cells). Partial correlations are computed individually for each ensemble member of the data and the models. The uncertainty bands in the data ensemble (Fig. 2c, d) represent the 5th and 95th percentiles of the partial correlations per latitudinal window. Processing Earth system model outputs from CMIP5. We analysed historical simulations outputs from ten Earth system models from CMIP5 65 (Supplementary  Table 3). The historical scenario simulations (also known as the 20th-century simulations) for CMIP5 were carried out for the period from the start of the industrial revolution to near present: 1850-2005. The Earth system model here is the atmosphereocean coupled global climate model coupled to a carbon-cycle model, and was forced in diagnostic mode by observed changes in atmospheric composition from natural and anthropogenic sources, volcanoes, greenhouse gases and aerosols, as well as changes in solar output and land cover. The model outputs evaluated relate to climate (temperature (tas), precipitation (pr), net and shortwave downward radiation (Rn and rsds, respectively)); carbon fluxes (net ecosystem exchange (nee), gross primary production (gpp), net primary production (npp) and autotrophic and heterotrophic respiration fluxes (ra and rh, respectively) to determine ecosystem respiration (reco 5 ra 1 rh); and carbon pools (accounting for leaf (cLeaf), wood (cWood) and roots (cRoot) in vegetation; and accounting for soil (cSoil), litter (cLitter) and woody debris (cCwd) in soil).
The spatial fields of the variables were obtained by computing mean annual values between 1982 and 2005. The ranges stand for the common period between data availability for GPP fluxes 44 and the historical runs from CMIP5. Like for the data, the modelled values of t are estimated from equation (1) using these simulation outputs.
Model outputs were always processed at the native spatial resolutions. To perform comparisons between models and between models and data, we constructed a common grid model ensemble. The common model ensemble was built by aggregating all model outputs to a common spatial grid, corresponding to the resolution of the NorESM1-M model (,1.89u by 2.5u, latitude by longitude). The aggregation consisted of computing an area-weighted mean per grid cell. Because each model shows a different number of model realizations (number of model ensembles in Supplementary Table 3), we averaged all ensembles per model to avoid overweighting models with a higher number of realizations. The original grid outputs were used for within-model evaluations (partial correlations, latitudinal gradients, global and biome statistics).