Apparent ecosystem carbon turnover time: uncertainties and robust features

The turnover time of terrestrial carbon (τ) controls the global carbon cycle – climate feedback and, yet, is poorly simulated by the current Earth System Models (ESMs). In this study, by assessing apparent carbon turnover time as the ratio 15 between carbon stocks and fluxes, we provide a new, updated ensemble of diagnostic terrestrial carbon turnover times and associated uncertainties on a global scale using multiple, state-of-the-art, observation-based datasets of soil organic carbon stock (Csoil), vegetation biomass (Cveg) and gross primary productivity (GPP). Using this new ensemble, we estimated the global average τ to be 42$% &' years when the full soil depth is considered, longer than the previous estimates of 23$) &* years. Only considering the top 1 m (assuming maximum active layer depth is up to 1 meter) of soil carbon in circumpolar regions 20 yields a global τ of 35$) &' years. Csoil in circumpolar regions account for two thirds of the total uncertainty in global τ estimates, whereas Csoil in non-circumpolar contributes merely 9.38%. GPP (2.25%) and Cveg (0.05%) contribute even less to the total uncertainty. Therefore, the high uncertainty in Csoil is the main factor behind the uncertainty in global τ, as reflected in the larger range of full-depth Csoil (3152-4372 PgC). The uncertainty is especially high in circumpolar regions with a uncertainty of 50% and the spatial correlations among different datasets are also low compared to other regions. Overall, we 25 argue that current global datasets do not support robust estimates of τ globally, for which we need clarification on variations of Csoil with soil depth and stronger estimates of Csoil in circumpolar regions. Despite the large variation in both magnitude and spatial patterns of τ, we identified robust features in the spatial patterns of τ that emerge regardless of soil depth and differences in data sources of Csoil, Cveg and GPP. Our findings show that the latitudinal gradients of τ are consistent across different datasets and soil depth. Furthermore, there is a strong consensus on the negative correlation between τ and 30 temperature along latitude that is stronger in temperate zones (30oN-60oN) than in subtropical and tropical zones (30oS30oN). The identified robust patterns can be used to infer the response of τ to climate and for constraining contemporaneous https://doi.org/10.5194/essd-2019-235 O pe n A cc es s Earth System Science Data D icu ssio n s Preprint. Discussion started: 22 January 2020 c © Author(s) 2020. CC BY 4.0 License.

initial photosynthetic fixation until respiratory or non-respiratory loss (Bolin and Rodhe, 1973;Barrett, 2002;. Ecosystem turnover time an emergent property that represent the macro-scale turnover rate of terrestrial carbon that emerges from different processes such as plant mortality and soil decomposition. Alongside photosynthetic fixation of carbon, τ is a critical ecosystem property that co-determines the terrestrial carbon storage and the carbon sink potential. As a 145 result of the balance between inputs and outputs of carbon, the terrestrial carbon pool can be approximated to reach the steady-state condition (inputs equal outputs) when long timescales are considered. This simplifies the calculation of τ to the ratio between the total terrestrial carbon storage and the influx or the outflux of carbon. The approach is advantageous to represent the highly heterogeneous intrinsic properties of the terrestrial carbon cycle as an averaged apparent ecosystem property which is more intuitive to infer large scale sensitivity of τ to climate change. Instead of focusing on the 150 heterogeneity of individual compartment turnover times we show the change of carbon cycle on the ecosystem level using τ as an emergent diagnostic property.
The magnitude of τ and its sensitivity to climate change is central to modelling carbon cycle dynamics. Therefore, τ has been used as an emergent ecosystem property to evaluate and constrain Earth system model (ESM) simulations of the carbon cycle. The current ensemble of ESMs shows a large spread in the simulation of soil and vegetation carbon stocks and its 155 spatial distribution, mostly attributed to the differences in τ among ESMs (Friend et al., 2014;Todd-Brown, 2013,2014Wenzel, et al., 2014Thurner et al., 2017). The large uncertainties in the simulated total carbon stock of soil and vegetation represent potential missing key processes that lead to diverse or even opposite response of τ to the rising temperature. Thus, it is instrumental to use observational-based estimations of carbon turnover times and their associated uncertainty in order to constrain the models and better predict the response of carbon cycle to climate change. 160 Current understanding of the factors that drive changes in τ are unclear due to the confounding effects of temperature and moisture even though it is well perceived that temperature and water availability are the main climate factors that affect root respiration and microbial decomposition (Raich, J. and W. H. Schlesinger,1992;Davidson and Janssens, 2006;. Therefore, it is difficult to implement local temperature sensitivity of τ into carbon cycle models due the large discrepancy between intrinsic and apparent sensitivity of τ to temperature. As the soil environment and climate are 165 highly heterogeneous in space, the temperature sensitivity of τ is substantially affected by other factors as spatial scale decreases (Jung et al., 2017).
Model simulations and observations do not agree in how the global distribution of τ is related to climate. Carvalhais et al. (2014) combined observational datasets that cover both low latitudes and circumpolar regions to estimate global τ and compared with CMIP5 simulations. They found a divergent result of global simulated total terrestrial carbon stocks that 170 range from 1101 Pg C to 3374 Pg C (mean difference of 36%) leading to a wide range of turnover times from 8.5 to 22.7 years (mean difference of 29%). The models also exhibit a large discrepancy in the τ-temperature and τ-precipitation relationships across different latitudes compared to observations. Koven et al. (2017) illustrated a higher sensitivity of τ to temperature in cold regions than in warm regions using an observational-based soil dataset. They found that most of the Deleted: It is ESMs fail to capture the global τ -temperature pattern. The difficulty of evaluating the response of soil carbon to climate 180 change is partly due to the fact that the dynamical observations at relevant timescales e.g. multi-decadal to centennial are lacking and the magnitude of projected change of τ to climate change is still poorly constrained (Koven et al., 2017).
There are not only large differences between simulations and observation-based estimates of τ, studies also show differences in the current observation-based estimates themselves. Specifically, estimates of global total carbon stock are characterized by large uncertainties as different in-situ measurements and methods were used to derive total carbon stocks (Batjes, 2016;185 Hengl et al., 2017;Sanderman et al., 2017). Alongside recent soil carbon datasets (Tifafi et al., 2018), there are also several different global vegetation biomass Avitabile et al., 2016;Saatchi et al., 2017;Santoro et al., 2018) and GPP (Gross Primary Production) (Jung et al., 2017) products which may lead to substantial differences in the global τ distribution and its relationship with climate. Thus, there is an urgent need to construct an ensemble of global τ estimates derived from different products and to quantify the uncertainty of the τ response to climate. 190 This study thus aims at developing an ensemble global estimation of τ in the spatial resolution of 10 kilometres (0.083º) which is derived from different observation-based products. Specifically, we will (1) update τ estimations with multiple state-of-the-art datasets; (2) quantify the contribution of the different components of τ to the global and local uncertainties; (3) identify the robust patterns across different ensemble members.

Datasets 195
The attributes of the τ dataset provided in this study, and the key external datasets that were used to estimate τ are summarized in Table 1. Details for each dataset are described in the below subsections.

Soil organic carbon datasets
Estimation of global soil carbon stock (Csoil) was based on five datasets that are derived from different approaches of estimating soil carbon: 200 a. SoilGrids is an automated soil mapping system that provides consistent spatial predictions of soil properties and types in the original spatial resolution of 250m . Global compilation of soil profiles is used to produce automated soil mapping based on machine learning algorithms. The data contains global soil organic carbon content at intervals of 0, 5, 15, 30, 60, 100 and 200 cm. In addition, chemical and physical properties such as bulk density and carbon concentration are provided. 158 remote-sensing based covariates including land cover classes, long-term 205 averaged surface temperature was used to fit the model. According to Hengl et al. (2017), the new version of the dataset can explain more of the variance (68.8%) in soil carbon stock than the previous version (22.9%) (Hengl et al., 2014).
However, it has also been recognized that the current version of SoilGrids (released on 2017.08.01, ftp://ftp.soilgrids.org/data/recent) may overestimate carbon stocks due to high values of bulk density (Tifafi et al., 2018).
In general, the estimation of Csoil is hampered by the availability of field data, especially in the circumpolar regions. 210 Deleted: available measurements Even though in-situ measurements had a large spatial extent and cover most of the continents, the regions that are characterized by severe climate or remoteness were much less sampled.
b. The dataset of soil carbon provided by Sanderman et al. (2017, hereafter Sanderman) used the same method as SoilGrids but different input covariates. The main difference between SoilGrids and Sanderman is that in addition to topographic, 215 lithological, climatic covariates, Sanderman also incorporated land use and forest fraction into the model fitting. The relative importance analysis based on Random Forest model shows that soil depth, temperature, elevation and topography are the most important predictors which is similar to the model fitting results of SoilGrids. Land use types such as grazing and cropping land area also contributes significantly to the variance. The Sanderman dataset provides soil carbon stocks at soil depths of 0-30 cm, 30-100 cm and 100-200 cm and at a spatial resolution of 10 km. 220 c. Harmonized World Soil Database (HWSD) was also used in this study which utilized over 16000 standardized soilmapping units worldwide which are harmonized into a global soil dataset (Batjes et al., 2016). HWSD is a 30 arc-second raster database that provides soil properties including organic carbon and water storage capacity at topsoil (0-30 cm) and subsoil (30-100 cm). HWSD combined regionally and nationally updated soil information worldwide to estimate soil properties in a harmonized way, and yet reliability of the data varies due to the different data sources. The database 225 which was derived from the Soil and Terrain (SOFTER) database had the highest reliability (Central and Eastern Europe, the Caribbean, Latin America, Southern and Eastern Africa) while the database that derived from the Soil Map of the World (North America, Australia, West Africa and Southern Asia) has a relatively lower reliability. d. We used the Northern Circumpolar Soil Carbon Database (NCSCD), which quantified soil organic carbon storage specifically in the northern circumpolar permafrost area (Hugelius et al., 2013). The dataset contained northern 230 circumpolar soil organic carbon content for depths of 0-30, 0-100, 100-200, 200-300 cm. The soil samplings included pedons from published literature, existing datasets and unpublished material. The 200 and 300 cm depth soil data was obtained by extrapolating the lowermost available values for bulk density and carbon content for a specific pedon to the full depth if the field data were only available in the first 50 cm of the full soil depth. However, the deep soil carbon (100-300 cm) showed the lowest level of confidence due to lack of in-situ measurements and much lower spatial 235 representativeness. The data was downloaded from https://bolin.su.se/data/ncscd/. e. The soil carbon stock and properties produced by the LandGIS maps development team (hereafter LandGIS) were also used in this study (Wheeler and Hengl, 2018). The soil profiles that were used in the training had a wide geographic coverage of America, Europe, Africa and Asia. One unique feature of LandGIS is that it included the soil profiles of Russia from the Dokuchaev Soil Science Institute/Ministry of Agriculture of Russia, which significantly improved the 240 predictions of Csoil in Russia. Different machine learning methods including random forest, gradient boosting and multinomial logistic regression were used to upscale the soil profiles to a global gridded dataset. Continuous 3D soil properties were predicted at 6 standard depths: 0, 10, 30, 60, 100 and 200cm. In comparison with the SoilGrids dataset, LandGIS added new remote sensing layers as covariates in the training and used 5 times more training points (360000 soil profiles) than SoilGrids (70000 soil profiles). The data was downloaded from 245 https://zenodo.org/record/2536040#.XYs1wpP7TUI.

Vegetation biomass datasets
Four different datasets of biomass at global scale were used to produce the total vegetation biomass (Cveg). 255 a. Thurner et al. (2014) estimated the above-ground biomass (AGB) and below-ground biomass (BGB) for northern hemisphere boreal and temperate forests (0.01° resolution, representative for the year 2010) based on satellite radar remote sensing retrievals of growing stock volume (GSV) and field measurements of wood density and biomass allometry. The carbon stocks of tree stems were estimated based on GSV retrieved with the BIOMASAR algorithm using remote sensing observations from the ASAR instrument on Envisat Satellite (Santoro et al., 2015), which was then 260 converted to biomass using wood density information. The other tree biomass compartments (BC) including roots, foliage and branches were estimated from stem biomass based on field measurements of biomass allometry. The total carbon content of the vegetation was derived as the sum of the biomass of the different compartments, which was then converted to carbon units using carbon fraction parameters. Comparison between the biomass map and inventory-based data shows good agreement at regional scales in Russia, the United States and Europe . Since data 265 from Thurner et al only covered northern boreal and temperate forests (30-80°N), we used data from Saatchi et al (2011) to cover the lower latitudes. b. We also incorporated a map of forest biomass carbon stocks for the tropical regions provided by Saatchi et al., (2011).
The map was derived using lidar, optical and microwave satellite imagery, trained using 4079 in-situ forest inventory plots (Saatchi et al., 2011). The method used GLAS Lidar observations to sample forest structure and used a power-law 270 functional relationship to estimate biomass from the Lidar-derived Lorey's height of the canopy. This extended sample of biomass density is then extrapolated over the landscape using MODIS and radar imagery, resulting in a pantropical AGB map. BGB was estimated as a function of AGB and the two were used together to derive total forest carbon stock at a 1 km spatial resolution.
c. The GlobBiomass map (Santoro et al., 2018) estimateed GSV and AGB density at global scale for the year 2010 at 100 275 m spatial resolution. The AGB was derived from GSV using spatially explicit Biomass Expansion and Conversion Factors (BCEF) obtained from an extensive dataset of wood density and compartment biomass measurements. GSV was estimated using space-borne SAR imagery (ALOS PALSAR and Envisat ASAR), Landsat-7, ICESAT LiDAR and auxiliary datasets, using the BIOMASAR algorithm to relate SAR backscattered intensity with GSV (Santoro et al., 2018b) 280 d. A pantropical AGB map (Avitabile et al., 2016) that combined two existing AGB datasets (Saatchi et al., 2011;Baccini et al., 2012) was also incorporated in the data ensemble. This map used a large independent reference biomass dataset to calibrate and optimally combine the two maps. The fusion approach was based on the bias removal and weightedaverage of the input maps, which incorporated the spatial patterns presented by the reference data in the fused map. The resulting map presents a total AGB stock for the tropics which was 9-18% lower than the two input maps and gave 285 different spatial patterns over large areas. The fused biomass map has a spatial resolution of 1 km.

Soil depth dataset
A full soil depth dataset was obtained from the Global Soil Texture and Derived Water-Holding Capacities database (Webb, 290 et al., 2000). Standardized values of soil depth and texture on a global scale, which were selected for the same soil types for each continent, were contained in the database. The full soil depth depends on soil texture and water availability which is usually higher than 100cm. In permafrost soil, full soil depth can extend beyond 400cm ( Figure S6).

The FLUXCOM global gross primary productivity dataset
FLUXCOM is an initiative to upscale biosphere-atmosphere fluxes measurements from eddy covariance flux towers 295 (FLUXNET) to global scale (Jung et al., 2017). In this study, we used the mean annual GPP datasets based on remotesensing forcing and nine machine learning methods with two flux partitioning methods trained on daily carbon fluxes, that is, 18 members of GPP (Tramontana et al., 2016). In order to produce high resolution (0.083º) spatial grids of carbon fluxes, only high-resolution satellite-based predictors were used in model training. In this study, we derived the long-term mean annual GPP by averaging annual GPP from 2001 to 2014. We note that all the 18 members is used independently to 300 estimated τ (not averaged).

Climate datasets
A high spatial resolution (~1km) climate dataset WorldClim 2 (Fick and Hijmans, 2017) was used to investigate the relationship between τ and climate. The data included monthly maximum, minimum and average temperature, precipitation, solar radiation, vapor pressure and wind speed. The data was produced by assimilating between 9000 to 60000 ground-305 station measurements and covariates such as topography, distance to the coast, and remote-sensing satellite products including maximum and minimum land surface temperature, and cloud cover in model fitting. For different regions and climate variables, different combinations of covariates were used. The two-fold cross-validation statistics showed a very high model accuracy for temperature-related variables (r > 0.99), and a moderately high accuracy for precipitation (r = 0.86).

Estimation of ecosystem turnover times 320
The total land carbon storage can be estimated by summing soil carbon stocks derived from extrapolation and vegetation biomass. Assuming steady state in which the total efflux (autotrophic and heterotrophic respiration, fire, etc.) equals to influx (GPP). Then τ can be calculated as the ratio between carbon stock and influx: Here Csoil and Cveg are the total soil and vegetation carbon stocks, respectively. We combined three soil carbon stock at three soil depth (1m, 2m, full soil depth), four vegetation, eighteen GPP that is, 648 members in total.

Estimation of global vegetation biomass stock 330
The aboveground biomass datasets only contain biomass of trees, meaning the herbaceous part is not considered. To account for herbaceous biomass, we used the same method as Carvarhais et al. (2014) which assumed the live vegetation fraction has a mean turnover time of one year, then using a uniformly distributed probability distribution of respiratory costs between 25 to 75 percent, we were able to relate GPP with the carbon stock of vegetation as: Where CH is the carbon stock for the herbaceous vegetation biomass; GPP is based on the FLUXCOM estimations; α is the percentage of respiration cost and fH is the fraction of herbaceous part for each grid cell based on the SYNMAP database (Jung et al., 2006).

Two vegetation biomass datasets (GlobBiomass and the Avitabile dataset) do not include BGB, in contrast to Saatchi's and
Thurner's products. In order to make all Cveg products comparable, we estimated the BGB from the empirical relationship 340 bewteen AGB and BGB derived previously by Saatchi et al. (2011):

Extrapolation of soil datasets
Extrapolation is necessary to obtain the accumulated carbon stock from surface to full soil depth because the soil datasets only extend to 2 meters below the surface. However, a large amount of Csoil is stored below this depth, especially in peatland 345 where soil carbon content is much higher in deeper soil (Hugelius et al., 2013). To estimate the total carbon storage in the land ecosystem, different empirical mathematical models were used (Table S1). The Covariance Matrix Adaptation Evolution Strategy (CMA-ES) method was used to optimize parameters of the models which is based on an evolutionary algorithm which used the pool of stochastically generated parameters of a model as the parents for the next generation (Hansen et al., 2001). 350 Extrapolation, which involves using empirical numerical models, may cause arbitrary bias and higher uncertainty if the models are not appropriately chosen. Here we used the in-situ observational data from the World Soil Information Service (WOSIS) (Batjes et al., 2019) and the International Soil Carbon Network (ISCN) (Nave et al., 2017) to select the ensemble of the models that could best simulate soil carbon stocks at full depth. The approach (i) fit of each empirical model against cumulative Csoil with all data points up to 2m; then (ii) predicted the cumulative Csoil at full soil depth (see section 2.3 for the 355 data) for each soil profile independently. The ability of a particular empirical model or combination of models (see Section 3.2) was then evaluated by comparing the predictions of Csoil at full depth against the observations. This procedure was applied on the two different in situ datasets, WOSIS which covers most of the biomes and ISCN which has more coverage in circumpolar regions. Finally, after comparing different model averaging methods (see supplement Table S2) we chose two model ensembles that could best represent circumpolar and non-circumpolar regions based on observational datasets, 360 respectively. The performance of the chosen ensembles is synthesized in Figure S3 and S4.

Uncertainty analysis
We performed a N-way ANOVA on different variables in order to calculate the uncertainties that stem from different data sources. The method can provide the sum of square variance and the total variance which derived the contribution of each data source to the total uncertainty. We defined the final contribution from each component involved in the calculation of τ 365 Where Cn is the contribution of uncertainty from a certain variable SSn, SStotal is the sum of square variance of all variables. 370 The contributions of uncertainties from soil, vegetation, GPP and soil depth of all ensemble members to the target variable τ were calculated. The uncertainty analysis reflects the relative spread of each group and the effect on the spread of τ.

The analysis of zonal correlations
The local correlation between τ and climate across latitudes was obtained by using a zonal moving window approach in which the Pearson partial correlations between τ and MAT/MAP were calculated using a 360º (longitudinal span) ×2.5º (latitudinal span) moving window. This approach allowed for the assessment of the relative importance for each climate parameter. The lowest and highest 1% of data points in each moving window was removed to avoid the effect of potential 380 outliers. In order to investigate the effect of latitudinal span, we chose different band size of 0.5º, 2.5º and 5º and performed the correlation analysis in the same manner for each selection. Table 2 shows the estimates of Csoil Cveg and GPP. Globally, estimates of soil carbon stocks within the top 2-meters of soil are 2749 PgC, 3628 PgC and 3546 PgC for the datasets of S2017, SoilGrids and LandGIS, respectively (bulk density corrected, see Supplement). The significant differences among different datasets indicate a high uncertainty in our current estimation of global soil carbon storage. The extrapolation of Csoil to full soil depth (FD) shows that approximately 18% of carbon stored is below 2 meters. Compared to the previous generation of soil data HWSD (available only within the top 1 390 meter), the state-of-the-art datasets of the current study have significantly higher carbon stocks within the top 1-meter of soil (Table 2). On the other hand, the current datasets of vegetation biomass show global Cveg ranges from 407 to 451 PgC and has less relative uncertainty than Csoil. The estimation of the uncertainty that derived from different GPP members shows a narrow range of 99 to 106 PgC from different products. Overall, the results show the difference in Csoil is much larger than Cveg and GPP. We next access the spatial distribution of soil, vegetation carbon and GPP in order to understand the 395 contribution of each component to the spatial uncertainties of τ.  in the circumpolar region. The results show that the relative difference of Csoil in circumpolar regions is two times larger than 410 that in non-circumpolar regions among all datasets.

The spatial distribution of soil carbon stock
The spatial distribution of Csoil is more consistent across datasets in the non-circumpolar region than in the circumpolar region ( Figure 1). The correlation coefficients (r) between each pair of datasets in the non-circumpolar region are generally higher than in the circumpolar region. Our results show a moderate agreement among the datasets in the spatial distribution 415 of Csoil globally (r>0.65). However, there are significant differences in the spatial patterns between the HWSD and each of the recent datasets ( Figure 1) as the correlation coefficients are all below 0.3. In addition, there is 2-fold lower carbon storage in the HWSD than the other datasets. Ratios between the total Csoil in the top 100 cm (Fig 1: upper off diagonal plots) show that LandGIS, SoilGrids and S2017 are consistent in temperate regions but show poor agreement in the tropical and the boreal regions. The comparison also shows that the gradient in carbon stocks between Europe and the lower latitudes 420 diminished in the HWSD soil map. In addition, the spatial distribution and the amount of carbon stocks in Indonesia is significantly different in the HWSD.
Higher dissimilarities of spatial patterns across the datasets in the circumpolar region is shown in Figure 2. We included the NCSCD dataset, which specifically focuses on the circumpolar region. The spatial correlations between each pair of the four datasets show low r values, which range from 0.2 to 0.5. In contrast with the non-circumpolar region, the high spatial 425 dissimilarity in circumpolar region indicates higher uncertainty regarding the estimation of total carbon storage. However, there is no clear evidence on which dataset is more credible in terms of total carbon storage and spatial pattern. The large differences are possibly due to fewer observational soil profiles in the northern high-latitude regions, which are crucial is the model training process.
The comparison of all datasets shows that there is a good agreement in the vertical structure of terrestrial carbon stocks. The 430 Csoil in the top 1-meter is about half of the total terrestrial carbon and 80% for the top 2-meter Csoil regardless of region or data source. For the non-circumpolar region, all the datasets show significantly higher carbon storage in the top 1m (451-635 PgC higher) than that in the HWSD, while showing less divergence of carbon storage among these three datasets (Table 2).
In general, the current datasets show similar vertical distribution of Csoil with consistent values and ratios between 1m and 2m soil. The extrapolation results indicate that about 15% of carbon is stored below 2m in the non-circumpolar region. For 435 the circumpolar region, the four datasets show a clear trend that the difference of Csoil increases with soil depth, as shown in Table 2. The difference between the top 1m Csoil among datasets has a higher difference than that of 2m. However, the ratio between storage in 1m and 2m is similar across all datasets.

The spatial distribution of vegetation and GPP 440
In comparison with soil carbon, the results show much more consistency and convergent global number of carbon stock among the four global vegetation datasets (Figure 3). Our results show that global vegetation carbon stock is 10% to 25% of the global soil carbon stock, depending on soil depth. The high spatial correlations (r>0.75) between each pair of data indicate the current estimations of vegetation is consistent. Different from the spatial distribution of soil carbon, most vegetation carbon is located in the tropics whereas much less carbon in higher latitudes. As a matter of fact, the Cveg in 445 circumpolar region is only 10% of that in non-circumpolar region (Table 2). Here we address that the Cveg is consist of three components including AGB, BGB and herbaceous biomass among which the herbaceous biomass is estimated from mean annual GPP (see Methods). We show the herbaceous biomass is only 5% of the total Cveg and less than 1% of the total Csoil indicating the minor role of herbaceous biomass in affecting the spatial distribution of total carbon stock and the uncertainty. Our results show that the spatial pattern of the global GPP is similar to the Cveg where there is higher primary productivity in the tropics and lower in the higher latitudes (Figure 4). The different members of the GPP estimation (see Methods) show very high consistency globally except for arid and polar region. The relative uncertainties in arid and polar region range from 455 50% to 100% whereas there is less than 50% of uncertainties in other regions. Although the differences among different vegetation and GPP estimations, in general, are not as high as soil, we show the uncertainties can be regionally high. We thus next investigate the contribution of each component to the spatial distribution 460 and uncertainty.

The ecosystem carbon turnover times and associated uncertainties
Using Csoil, Cveg and GPP, we estimated the carbon turnover times with different combinations of datasets in order to quantify the uncertainty. We calculated τ in the same manner as the previous study  in which they used only 1m of soil in the circumpolar region and full soil depth in the non-circumpolar region. We compared the spatial 465 distribution of the τ estimations of the previous and our study ( Figure S5). Our results show a large range of relative difference and low spatial correlation (r = 0.51). We found the main differences are in the northern circumpolar region, which is caused by the differences in the Csoil estimations, indicating a large uncertainty of τ estimations the in this region.
Our estimation of global mean τ is 35 years with an interquartile range of 29 to 43 years, which is much longer than the previous study of 23 years with interquartile range of 19 to 30 years. In addition, we derived a global τ of 40 years with an 470 interquartile range of 29 to 47 years by assuming the maximum active layer thickness to be the full soil depth in the circumpolar regions instead of using only 1-meter Csoil as was done in the previous study. The incorporation of deep soil in the circumpolar region increased the global τ by 7 years. The global spatial distribution of τ ( Figure 5) shows great heterogeneity, which ranges from 5 years in the tropics to over 1000 years in northern high latitudes. The results show a Ushaped distribution of τ along latitudes where τ increases nearly three orders of magnitude from low to high latitudes. Figure  475 5b shows the map of relative uncertainty that is derived from different datasets. The higher relative uncertainty indicates more spread among the datasets used to estimate τ. Our result shows that peatland and arid regions generally have higher uncertainties than the rest of the world. We found several regions with very different estimations of τ among the datasets including north-east Canada, central Russia and central Australia where the relative uncertainties are over 100%. Csoil, Cveg and GPP contribute differently to the overall uncertainty of τ as shown in Figure 6. The difference among soil datasets is the dominating factor of τ uncertainty, especially in the circumpolar regions and the Indonesian peatland where there is large amount of soil organic carbon in subsoil. On the other hand, the uncertainties of τ in arid and semi-arid regions are controlled by the difference in GPP products. The contribution of vegetation to the uncertainty in τ is most significant in 495 the tropics and warm temperate regions where there is large vegetation biomass. It is worth to note that contributions from each component also vary with depth of carbon stock that was used to calculate τ. For instance, the uncertainty contribution from Cveg becomes smaller when the Csoil up to 2 meters is used compared to only using 1-meter in calculating τ. However, the fact that the difference in the soil products was the major contributor to the τ uncertainty remains no matter what soil depth is used. Globally, the uncertainty of τ is mostly derived from soil and GPP, which dominate 82% and 17% of the 500 global land area, while vegetation plays a minor role globally (1%).

The zonal pattern of turnover times
The latitudinal distributions of τ can be best represented by a second-degree polynomial function (Figure 7b). After fitting the data of all ensemble members, the rate of τ change with latitude can be obtained by taking the first derivative of the fitted polynomial function. We found that the rate of τ change (Figure 7c) has very consistent zonal patterns for different τ 505 ensemble members from different data sources. The result shows a consensus on the change of τ with latitude of different datasets. We also found that the zonal τ gradients were not significantly (P > 0.05) different from each other for different selections of soil depth, indicating soil depth has no significant effect on the τ gradient along latitude. It is worth to note that there is a significant difference in the zonal τ gradient between the northern and southern hemisphere (P < 0.0001) and that τ increases faster from low to high latitude in northern latitudes than in the southern latitudes. The results show that we have 510 high confidence in the zonal distribution of τ and that the difference across datasets does not affect the robustness of the pattern.

The zonal correlation between turnover time and climate
The correlations between τ and mean annual temperature and mean annual precipitation are analysed for all the ensemble members on global scale (see Method section). The correlation (Figure 8a) is the strongest in northern mid-to-high latitudes 515 between 25º N and 60 º N, and it decreases rapidly from 20º N to the equator. In the southern hemisphere, it increases until 40º S, albeit having a weaker gradient than in the northern hemisphere. The uncertainties originating from different data sources are shown by the shaded area (Figure 8). The result shows that there are high uncertainties in the transitional regions between the temperate and Arctic regions (50 -70º N) as well as tropical regions (20º N to 20º S). Similar to the previous result of uncertainty contribution where soil is the dominating factor, the differences in Csoil also cause the spread in τ -T 520 correlation. However, the patterns of correlation along latitude do not change regardless of the data source and the soil depth.
All ensemble members agree that τ is negatively associated with temperature, with stronger associations in cold regions than in warm regions.
The correlation between τ and precipitation, in general, has larger variability across latitude and a higher uncertainty due to differences in data (Figure 8b). Contrary to the τ -T relationship, the uncertainty of the τ -P relationship derived from both different data sources and soil depths are smaller in the tropics than in high latitudes. Negative correlations dominate the 535 latitudes between 20 and 50º N as well as between 20 and 40º S, while there is a stronger positive correlation in the tropics.
There is a shift in the sign of the correlation coefficient from negative in temperate zone to positive in tropics, indicating the role of water changes from water-limited regions to water-excessive regions. We found the pattern of correlation between τ and precipitation is different from the previous study , especially in the tropics. We therefore investigated the possible cause of the difference by mixing all components (Csoil, Cveg and GPP) between the previous and 540 current study. By examining the correlation between each mixed τ estimation and climate factors, our results show that positive correlation in the tropics is caused by the Csoil ( Figure S7). This is consistent with the previous results ( Figure 1) which shows large difference in the spatial distribution of Csoil in the tropics between the three soil datasets we used in this study and HWSD soil dataset.

Discussion 545
In this section, we will discuss the robustness of the current state-of-the-art estimation on global terrestrial carbon turnover times and their response to climate change. We first show the variation of spatial and vertical distribution of carbon stock in different regions and the possible reason for the difference, and we then discuss the robustness of zonal distribution of turnover times and zonal changing rates across different datasets. Finally, we focus on the sensitivity of turnover times to climate and implications. 550

Estimation of global carbon stock
Accurate estimation of terrestrial carbon storage and turnover time are essential for understanding carbon cycle-climate feedback (Saatchi et al., 2011;Jobbágy et al., 2000). Our analysis benchmarks soil, vegetation carbon storage and GPP from multiple state-of-the-art observational based datasets at global scale and provide not only an estimate of the total carbon stock but also the vertical distribution and spatial variability of global carbon stock. We divide the global map into 555 circumpolar and non-circumpolar regions due to the different characteristics and uncertainty.
We found that there is a significant difference across the current soil carbon datasets in both circumpolar and noncircumpolar regions. The results show that the uncertainty of Csoil estimations in the circumpolar region is two times larger than that of the non-circumpolar region ( Figure 6). The spatial patterns of total ecosystem Csoil among the soil datasets are more consistent in the non-circumpolar region than in the circumpolar region. In contrast with the non-circumpolar region, 560 there is lower confidence in the circumpolar region in estimating Csoil due the fact that there is low spatial correlation across datasets. The difference can be caused by various reasons. As an important input to the machine learning method, in-situ soil profiles are very important factors that influence the final results of the upscaling. The sparse coverage of soil profiles in the circumpolar region may cause the large divergence in the northern circumpolar region. A major difference between S2017 and the other two soil datasets is that soil carbon stock was a direct target of upscaling in the former dataset, while in the 570 latter two datasets each component used to calculate Csoil (carbon density, bulk density and percentage of coarse fragments) was predicted individually. In addition, the climatic covariates that were used in the upscaling were different (see Method).
In contrast with the non-circumpolar region, the circumpolar Csoil does not have a decreasing trend up to 4 meters of soil depth ( Figure S1) which indicates that there is a significant amount of carbon stores in deep soil. The deep soil turnover is a key process to the global carbon cycle yet poorly understood (Todd-Brown et al., 2013). In this study, we extrapolated the 575 soil carbon stock to full soil depth. We chose the model ensembles from a framework to pick out the models that had a minimum distance between prediction and observations by using in-situ soil profiles (see Supplement). Two model ensembles were selected that can best represent the soil vertical distribution in circumpolar and non-circumpolar regions by comparing model simulations and in-situ observations. The final results depend on the information from the soil profiles and also the characteristics of the empirical models. The extrapolation gave us insights to the carbon storage and vertical 580 distribution in deep soil. The results of extrapolation show there is approximately 15% of carbon stored below 2 meters globally and over 20% of carbon stored below 2 meters in the northern circumpolar region. Although the total amount of carbon storage in the ecosystem shows a large divergence among different datasets, the ratio between different soil depths are quite consistent indicating a high confidence in the vertical structure of soil compare to the total amount.
The global soil carbon stocks across observational-based datasets are much less divergent than the current earth system 585 model (ESM) simulations. The CMIP5 results show the simulated carbon storage ranges from 500-3000 PgC making τ varies by a factor of 3.6, from 11 to 39 years (Todd-Brown et al., 2013). Our results show that the amount of carbon in the ecosystem is much higher than the estimation by ESMs. Even the lowest estimation (S2017 dataset) of total carbon storage is about 500 PgC higher than the highest ESM estimation (MPI-ESM-LR). The spatial distribution of carbon stocks among ESMs have a large variation across models while the observational-based datasets are more consistent in the non-590 circumpolar region. But we leave a question mark to the soil carbon in the circumpolar region, which is characterized by large uncertainty as shown by the current observational-based soil datasets.
Compared with soil, the higher level of consistency in the vegetation and GPP estimates indicate there is a consensus on the current estimations in the above-ground carbon stock. However, we note that the regional differences in the products can significantly affect the spatial distribution and uncertainty of τ. Nevertheless, vegetation and GPP contributes a little to the 595 global mean value of τ estimations.

The terrestrial carbon turnover time and uncertainty
The uncertainty analysis showed that our current estimation of τ has a considerable spread which derived from the state-ofthe-art observations of soil, vegetation and carbon fluxes. In this study, we showed the uncertainty is contributed mainly by the soil carbon stock and GPP, where the former dominates the vast areas in the circumpolar region and the tropical peatland 600 while the latter dominates the semi-arid and arid regions. We showed that GPP is the second largest contributor to the total uncertainty which potentially leads to significant differences in the estimation of τ considering different products of GPP (Zhang et al., 2017;Norton et al., 2019). This result is consistent with the previous study (Todd-Brown et al., 2013) that the bias in estimated primary productivity can affect the carbon turnover estimations to a large extent not only by using observational-based data but in the ESMs simulations. However, the uncertainty comes not only from the differences across 605 datasets but also from the soil depth we chose to estimate τ. The frozen permafrost soil in the circumpolar region, although containing a large amount of carbon is an important component in the process of turnover (Zimov et al., 2006). However, we do not know to what soil depth we should use in the τ estimation since currently our knowledge on the active layer thickness of frozen permafrost soil is still lacking. In addition, the active layer thickness of permafrost changes with climate, which adds more uncertainty to the estimation of τ. Thus, we argue that the current datasets cannot support robust estimation of 610 global τ. It is worth to note that our estimation of τ is based on the steady-state assumption, that is, the net exchange of carbon between the terrestrial ecosystem and the atmosphere equals to zero. In our study, the steady-state assumption is a proper assumption for that our analysis focused on the τ estimation at long-term temporal scale and large spatial scale.
However, this assumption is valid to a much less extent at site-level as the net exchange of carbon is, most of the time, not in balance (Ge et al., 2018). 615 Although the current estimation of τ has a large variation, we show that the zonal distribution of τ is a robust feature that changes little with different datasets, which indicates that the current state-of-the-art datasets all agree on the latitudinal gradient of the carbon turnover time. Another robust feature is that the zonal changing rate of τ does not change with the soil depth ( Figure 7). It has always been a problem of what soil depth should we use to represent the functional part of carbon in the ecosystem. The selection of soil depth is usually arbitrary and varies from study to study. For example, Koven et al. 620 (2017) and Wang et al. (2017) used the top one-meter of soil carbon to represent the total terrestrial carbon pool while Carvalhais et al. (2014) extrapolated soil to full depth and used it as the pool. Our results demonstrate that the selection of the soil depth does not affect the zonal pattern that we observed. This can be better seen in the next section with the response of τ to climate.

Robust associations of τ and climate 625
Despite the large uncertainty in the τ estimations, we identified robust response of τ to climate change . It is well recognized that the sensitivity of terrestrial carbon to climate is a major uncertainty, which is reflected by the spread of τ estimation by the different ESMs. However, we need reliable estimations of τ to quantify its climate sensitivity and provide robust constraints to improve the performance of the current ESMs. We showed the zonal correlation between τ and temperature varies with latitude where high correlations are found in the high latitude and low to moderate correlation in low latitude, 630 especially the tropics. The zonal pattern of τ-precipitation is more complicated in that water availability can cause local variability to a great extent. The correlation between turnover times and precipitation in the tropics is higher than that with temperature as shown in Figure 8d indicating a potentially more dominant role of precipitation in the tropics. The role shifted along latitude between temperature and precipitation in the pattern of τ due to the variation in the relative importance for each parameter. However, the temperature gradient shaped the zonal distribution of τ as it can be seen that τ increases with 635 latitude. All of these relationships are verified by each ensemble member of the data. We found the correlations, although they vary in strength, are very robust. The intimate interaction of energy and water along with other factors such as land use change all affect τ but on different spatial and temporal scales. It is worth mentioning that the τ -T relationship is similar when compared with previous results  whereas there are considerable differences in the τprecipitation relationship, specifically in tropical regions where the turnover times were always negatively correlated with 650 precipitation in previous study. The different τ -P zonal patterns of correlation between the previous and the current study, as shown before, is mainly caused by the difference in the soil carbon stock ( Figure S7). This finding indicates the response of τ to moisture is characterized by large uncertainty.

Conclusion 660
A full assessment of the global turnover times of carbon is provided using an observational-based ensemble of current stateof-the-art datasets of soil carbon stocks, vegetation biomass and GPP. At the global scale, the uncertainties in τ estimates are dominated by the large uncertainties in soil carbon stocks. The uncertainty of carbon stocks and τ estimation in the circumpolar region is significantly higher than that in the non-circumpolar region. Our results show that there is a consistent vertical distribution of soil carbon across datasets, and it is estimated that soils below 2 meters take up to 20% of total soil 665 carbon globally. A spatial analysis shows that both soil carbon and GPP are the major contributors of local uncertainties in τ estimation. The differences in soil stocks between datasets dominates the uncertainties of τ in the circumpolar region and in tropical peatlands, while the spread in GPP dominates the uncertainty in semi-arid and arid regions. The difference in vegetation data has a minor contribution to the uncertainty.
Despite the differences, we identified several robust patterns that change only marginally across different ensemble members 670 of τ that derived from different datasets or different soil depths. First, we found a consistent latitudinal pattern in τ that can be described by a second-degree polynomial function. The changing rate of τ with latitude can be described equally well for all ensemble members and the changing rate of τ with latitude is highly consistent across different datasets and does not change with soil depth. The same zonal correlations between τ and climate showed there is a robust association of τ with temperature and with precipitation. However, we note that association between temperature/precipitation and τ change with 675 latitude. Specifically, temperature mainly affects the τ variation in middle to high latitudes beyond 20ºN and 20ºS while precipitation affects τ not only in temperate zones but also in the tropical regions. Thus, the sensitivity of τ to a certain climate factor makes more sense to be calculated and interpreted in regions where the climate factor is the main driver of the τ variation. Overall, this study synthesizes the current state-of-the-art data on global carbon turnover estimation and argues Deleted: that the zonal distribution of τ and its response to climate is robust regardless of different datasets or assumptions on soil depth. This is a critical advancement since previous studies usually made arbitrary decisions on the soil depth that use to estimate τ, and supports exercises for benchmarking ESMs. Future studies should further investigate τ with regional spatial scale and its response to climate as well as other factors such as land use change, in order to have an in-depth understanding 685 of carbon cycle turnover.

Author contributions
NF is the main author who wrote the manuscript. NC and MR supervised the analysis of the data and development of the paper. SK intensively participate in the discussion and development of the paper. MT, VA, MS and BA provide feedbacks 690 on the scientific analysis and text of the paper. UW is technical support and responsible for data management.   In order to unify the spatial resolution and geographic coordinate system of dataset from different sources, we need to make sure that the total amount of stock for soil, vegetation, etc. doesn't change during aggregation and transformation, i.e. the variable need to be mass conservative. However, 'state' variable such as temperature, vegetation types do not need to fulfill the mass conservative requirement, nor they should. In our study, we developed a mass conservative method to maintain 905 mass for carbon stocks. We first multiply the variable that need to be aggregated (Xfine) by corresponding land area (Afine) at grid cell level represented by equation (1), then aggregate the product (XAfine) by summing the values in N×N grids cell depending on the target resolution (equation (2)). The land area is also sum to the target resolution (equation (3)). Finally, the area-weighted variable is derived by dividing aggregated product (XAcoarse) by corresponding land area (Acoarse) as illustrated by equation (4). 910 (1) We applied the method to all datasets that requires aggregation including soil, vegetation and GPP that were used in the 915 study.

Bulk density correction
The bulk density (BD) in SoilGrids and LandGIS are too high due to two reasons. First, the measurements of BD are less and missing in many horizons . And the measurements of BD in permafrost region, especially in Canada forest soil and Russian, are problematic (personal communication with Tomislav Hengl). In this study, we applied a 920 pedotransfer function from Köchy et al. (2015) to make correction based on organic carbon concentration (we only applied the function to the grid cells where carbon > 8%): In this section, we introduce the framework that we used to select the models for extrapolating soil from 0 -2m to full soil depth.

Different characteristics of permafrost and non-permafrost soil
The amount and vertical distribution of soil organic carbon are largely influenced by vegetation which fixes atmospheric 930 CO2 and transport carbon into the land ecosystem. However, the SOC stock have a much more complicated relationship with productivity of plants than a simple linearly one . The higher biomass, which implies more carbon sequestration by aboveground biomass, however, does not necessarily lead to increases in SOC storage. Although the processes of soil formation, accumulation, and stabilization have been intensively studied and debated, the mechanisms that determine the soil carbon stock, especially in deeper soil, are still unclear. Instead of using process modelling approach, we 935 chose statistical approach to extrapolate each soil profiles in the gridded dataset from 2m to full depth. The reason of performing soil carbon stock extrapolation is that we have little knowledge on how much the carbon stored in the soil that is deeper than 2m, although deeper soil is a crucial component in the climate-carbon cycle feedback. The other reason is the different dataset report SOC stock at different depths. The advantage of using statistical method is that we do not need to know the mechanisms that control the soil processes. Instead, we select simple empirical mathematical models that can 940 represent and predict the in-situ soil profiles.
We used 425 permafrost peatland profiles from ISCN soil database and 1000 profiles from WOSIS soil database to study the characteristics of vertical distribution of SOC. Figure S1 shows the accumulated SOC stock profiles with depths in permafrost and non-permafrost region. The vertical distribution of carbon with depth in permafrost soil has a distinguished feature that the SOC has a high linear relationship with depth. This fact implies the soil carbon keeps increasing even after 3 945 meters in permafrost soil ( Figure S1b). However, we have no idea to what depth can soil carbon keep increasing and the total amount of the storage in permafrost peatland due to the limited observational depth of SOC. In contrast, soil profiles in nonpermafrost region stop increasing mostly before 2 meters. The results demonstrate the necessity of extrapolating soil to full depth, especially for permafrost soil.

Selection of models 950
We included 12 models (Table S1) for predicting SOC stock to full soil depth. Figure S2 shows an example result in which the data points that is shallower than 1m were used to fit all the models and predict the point that is deeper than 2m for a typical soil profile. Due to the different mathematical characteristics of the models, the prediction has quite a spread.
Relatively 'conservative' models including model ensemble BHIJKL tend to underestimate the carbon stock while the more 'aggressive' ones ACDEFG tend to overestimate the stock. 955 In the sense that we do not know which one or group of models can best predict the accumulated carbon storage, we conducted a selection process (see Methods) by grouping all the models into all possible combination and rank the performance for all the model averaging results as shown in Table 2. All the models were used to fit the WOSIS data which covers most of the biomes and ISCN database which covers only permafrost soil. We conducted three batch of experiments in the same manner but used data points within different depths. The data points lower than 50cm, 100cm and 200cm were 960 used to predict the SOC that is higher than 200cm. Our goal is to find the ensemble of models that has the highest model performance, the best coverage, the minimum error and AIC.

Extrapolating soil with different method
The main goal of using several models is to search for the best group of models that can best predict the vertical distribution 965 of soil carbon stock and we compared the below approaches for that purpose: 1. The Bayesian Model Averaging (BMA) method is used in this study to find the best model ensemble for the prediction of soil carbon storage to full depth. The MODELAVG Matlab toolbox (Vrugt, 2016) which implemented many different model averaging techniques including BMA method. The advantage of BMA method is that it considers explicitly the uncertainty of prediction of a target variable which can provide a probabilistic distribution of weight for each model 970 instead of only a weighted-average, deterministic prediction. By maximize the likelihood function from the training dataset, the weights b = {b1, …,bk} and standard deviation s = {s1, …, sk} are estimated.
2. Equal weights averaging (EWA) which consider each the participating model have the same weight and the prediction is derived by equal-weighted averaging the model results.

975
The complete combination among different models are also compared and the best model ensembles are obtained by maximizes MEF, minimizes KL and minimizes AIC. The results show that EWA and BMA methods have similar performances (Table S2). We choose EWA method due to it have a slightly better coverage of observations. Two model ensembles were selected from the model selection framework that can best represent circumpolar and noncircumpolar region based on observational datasets in the two regions, respectively. The performance of the chosen 980 ensemble is synthesized in Figure S3. It shows the ensemble DIJKL overall can well predict the carbon stock in noncircumpolar region that is deeper than 200cm only using points lower than 50cm. Model efficiency is 0.83 and the residue between observation and prediction is little biasd in the prediction ( Figure S3c). The histogram ( Figure S3b) shows the prediction has the same distribution as the observations. The ensemble was also used to predict observations within different percentiles and over 70% of observations can be included in the uncertainty ([-σ, +σ]). The results show that the selected 985 models can well represent the vertical distribution of Csoil thus we used them to extrapolate the global gridded datasets in order to obtain the total soil carbon storage in the soil. The selected model ensemble ACDEF for the circumpolar soil have lower model efficiency and less well represent the soil in the region ( Figure S4). We then applied extrapolation on three global datasets which are Sanderman, SoilGrids and LandGIS. The averaged results of ensemble DIJKL is used to extrapolate non-circumpolar soil from 2m to full soil depth and ensemble ACDEF to extrapolate circumpolar soil. 990    Figure S15: The zonal pattern of the correlation between τ and climate factors. Each component (Csoil, Cveg and GPP) from the previous study  is mixed with each component of the current study. The prefix 'old' stands for the component from the previous study and the prefix 'new' stands for the component from the current study.