Mapping global distributions, environmental controls, and uncertainties of apparent topsoil and subsoil organic carbon turnover times

The turnover time (τ) of global soil organic carbon is central to the functioning of terrestrial ecosystems. Yet our spatially explicit understanding of the depth-dependent variations and environmental controls of τ at a global scale remains incomplete. In this study, we combine multiple state-of-the-art observation-based datasets, including over 90 000 geo-referenced soil profiles, the latest root observations distributed globally, and large numbers of satellite-derived environmental variables, to generate global maps of apparent τ in topsoil (0–0.3 m) and subsoil (0.3–1 m) layers, with a spatial resolution of 30 arcsec (∼1 km at the Equator). We show that subsoil τ (${\mathrm{385}}_{\mathrm{20}}^{\mathrm{3485}}$ years (mean, with a variation range from the 2.5th to 97.5th percentile)) is over 8 times longer than topsoil τ (${\mathrm{15}}_{\mathrm{11}}^{\mathrm{137}}$ years). The cross-validation shows that the fitted machine learning models effectively captured the variabilities in τ, with R² values of 0.87 and 0.70 for topsoil and subsoil τ mapping, respectively. The prediction uncertainties of the τ maps were quantified for better user applications. The environmental controls on topsoil and subsoil τ were investigated at global, biome, and local scales. Our analyses illustrate the ways in which temperature, water availability, physio-chemical properties, and depth jointly exert impacts on τ. The data-driven approaches allow us to identify their interactions, thereby enriching our comprehension of mechanisms driving nonlinear τ–environment relationships at global to local scales. The distributions of dominating factors of τ at local scales were mapped for purposes of identifying context-dependent controls on τ across different regions. We further reveal that the current Earth system models may underestimate τ by comparing model-derived maps with our observation-derived τ maps. The resulting maps, with new insights, as demonstrated in this study, will facilitate future modelling efforts relating to carbon cycle–climate feedbacks and support effective carbon management. The dataset is archived and freely available at https://doi.org/10.5281/zenodo.14560239 (Zhang, 2025a).

1 Introduction

As the largest active reservoir of organic carbon in terrestrial ecosystems, soils are integral to the global carbon cycle (Schimel, 1995; Batjes, 1996; Smith et al., 2024). Plants capture CO₂ from the atmosphere through photosynthesis and transfer carbon into soils through litterfall and root exudates. This carbon is then cycled back to the atmosphere through heterotrophic respiration by organic matter decomposers, a process governed by decomposition rates (Balesdent et al., 2018). The turnover time of soil organic carbon (SOC), denoted as τ (in years), is the average time that organic carbon molecules remain in the soil (Six and Jastrow, 2002; Sierra et al., 2017). This turnover time is a critical factor in determining the size of soil carbon pools (Sierra et al., 2017; Crowther et al., 2019b). Understanding the spatial variation in τ and the underlying environmental drivers is therefore crucial for comprehending the scale and dynamics of terrestrial carbon storage under current and future climate change scenarios (Torn et al., 1997, 2009; Field et al., 2014).

Several studies have estimated the global apparent carbon turnover times in whole terrestrial ecosystems (Carvalhais et al., 2014; Fan et al., 2020). Nevertheless, the detailed patterns of τ in soil systems at the global scale remain to be elucidated. Much work has focused on shallow soil horizons due to the higher carbon content and the greater availability of data in topsoils (e.g. Crowther et al., 2016; Luo et al., 2017; Viscarra Rossel et al., 2019; Wu et al., 2021). However, capturing the biogeographic variability of carbon dynamics in deeper soil layers is emerging as a critical area of research (Hicks Pries et al., 2023) as, there, the environmental sensitivities of SOC can differ substantially from those in surface layers (Rumpel and Kögel-Knabner, 2011; Hicks Pries et al., 2017; Luo et al., 2019; Soong et al., 2021; Zosso et al., 2023). To enhance our understanding of soil carbon turnover and to address these issues, there is a clear need for a global analysis of τ that considers both topsoil and subsoil layers and for the provision of spatially explicit τ maps. Such an analysis should leverage the latest and most comprehensive global soil profile datasets available, enabling more reliable assessments of and insights into terrestrial carbon sink potential.

The ensemble of the latest datasets also supports a more nuanced understanding of the environmental controls on τ, particularly with regard to their potentially nonlinear or distinct effects across different spatial scales and soil depths. Previous studies have shown evidence that the turnover time of SOC is negatively correlated with temperature and precipitation (Davidson and Janssens, 2006; Chen et al., 2013; Wang et al., 2018) and that subsoil carbon may be particularly sensitive to temperature fluctuations (Jia et al., 2019; Soong et al., 2021; Chen et al., 2023). Yet, other studies could not confirm this strong climatic dependency (Giardina and Ryan, 2000; Doetterl et al., 2015), suggesting a predominant influence of soil properties on subsoil carbon turnover times in certain regions or over decadal timescales (Luo et al., 2019). This inconsistency underscores the need for a multifaceted approach to quantify the effects of multiple factors on τ, particularly the interactions between climate and edaphic factors across different spatial scales and soil depths (Schmidt et al., 2011). Comprehensive assessments at the global, biome, and local levels will be crucial to identify the primary controls of τ in both topsoil and subsoil layers. Moreover, observation-based global estimates of τ are essential for simulating the global carbon cycle (Todd-Brown et al., 2013; Friend et al., 2014; Varney et al., 2022). An accurate representation and deeper understanding of the environmental controls of τ – spanning diverse spatial scales and soil depths – will be integral to benchmarking current Earth system models (ESMs) and reducing bias in future carbon cycle projections.

This study aims to develop a global estimation of τ by integrating the state-of-the-art soil and root profile databases with satellite-derived environmental observations. The collected datasets allowed us, firstly, to estimate τ at over 90 000 global sampling sites, and then we used machine learning methods to generate a spatially explicit understanding of global SOC turnover times in the topsoil (0–0.3 m) and subsoil (0.3–1 m). To comprehend the interactive mechanisms among multiple environmental drivers that have shaped variations in topsoil and subsoil τ at the global, biome, and local scales, we used data-driven approaches to characterize the directional contributions of climate, topography, and the physical and chemical properties of soil to explain τ patterns. We further quantified the uncertainty maps for better user applications and compared our observation-derived τ with ESM-derived τ in terms of both spatial variabilities and climatic dependencies.

2 Materials and methods

2.1 Estimation of τ in topsoil and subsoil layers

The SOC turnover time (τ) we estimated in this study is the mean transit time of carbon (i.e. the time that the newly entered carbon spends in soils until it leaves) (Six and Jastrow, 2002; Sierra et al., 2017). When assuming a steady state and homogeneity of the system, in which all particles have the same probability of leaving at any time, τ can be defined as the ratio between SOC stock (SOCS) to the input or output flux of carbon ($\mathit{\tau }=\text{SOCS}/\text{flux}$), which can be called the apparent turnover time (Carvalhais et al., 2014; Fan et al., 2020). We estimated the SOC stock (see Sect. 2.1.1) and vertical allocation of carbon input in two soil layers (see Sect. 2.1.2) and then combined them to calculate topsoil and subsoil τ at all soil profile locations. We then adopted the point-level estimates to map τ variation across the globe.

2.1.1 Estimation of SOC stocks based on soil sample databases

Soil sample data were mainly collected from the standardized soil profile database provided by the World Soil Information Service (WoSIS) (Batjes et al., 2020). The Northern Circumpolar Soil Carbon Database (NCSCD) (Hugelius et al., 2013) was incorporated to supplement samples across the high latitudes of the Northern Hemisphere. In addition, the soil sample data from the project of the National Soil Survey of China were incorporated (Zhang et al., 2013; Liu et al., 2022). To minimize bias due to erroneous and/or uncertain measurements, we removed samples (i) with low accuracy of the geographical coordinates (i.e. the information of degree, minute, and second is not fully provided in the source), (ii) located in areas with exceptionally low NPP (net primary productivity) estimates (below 10 $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$), (iii) with low-quality NPP estimates (if the quality control value is larger than 50 %), and (iv) with layer observations flagged as surficial litter. Finally, a total number of 95 200 and 66 807 geo-referenced sampling locations providing topsoil and subsoil information were collected for this study (Fig. 1). The estimation of SOC stock (SOCS, kg C m⁻²) for each soil sample site at a certain depth interval between the upper depth (D_u, m) and lower depth (D_l, m) can be computed as follows:

$$\begin{array}{}\text{(1)}& {\text{SOCS}}_{D_{\mathrm{u}}-D_{\mathrm{l}}}={\text{SOC}}_{D_{\mathrm{u}}-D_{\mathrm{l}}}\cdot \text{BD}\cdot \left(\mathrm{1}-{\displaystyle \frac{\text{CF}}{\mathrm{100}}}\right)\cdot (D_{\mathrm{u}}-D_{\mathrm{l}}),\end{array}$$

where ${\text{SOC}}_{D_{\mathrm{u}}-D_{\mathrm{l}}}$ is the SOC content (g kg⁻¹) in a certain range of soil depths, BD is the soil bulk density (g cm⁻³), and CF is the percentage of coarse fragments in the entirety of the soil (%). We adopted an approach described in Sect. S1 in the Supplement to fill the missing values by using a specific pedotransfer function depending on the SOC content and the considered layer. When CF was missing, we extracted the corresponding data from their respective locations on the latest version of the SoilGrids maps (Poggio et al., 2021). The equal-area spline algorithm was employed to fit observations along depths (Bishop et al., 1999; Malone et al., 2009). The average of fitted values was adopted to estimate the SOCS in the topsoil and subsoil layers for each profile (Fig. S1 in the Supplement). The probability distributions of SOCS, BD, and CF in different biomes are shown in Figs. S2–S4.

Figure 1 Spatial distributions of geo-referenced soil samples used in the study across global biomes. A total of 95 200 (a) (for topsoil) and 66 807 (b) (for subsoil) sampling sites were collected from multiple soil profile databases. The bar plots show the sample size within each biome.

2.1.2 Estimation of carbon allocation belowground

Carbon influx at each soil layer comes from vertical allocation of the amount of net primary productivity (NPP). This involves the estimation of aboveground and belowground NPP and NPP allocation into different soil layers. The annual NPP ($\mathrm{kg}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$) produced by the Moderate Resolution Imaging Spectroradiometer (MODIS) was collected at each sample location (Running and Zhao, 2019). The MOD17 version 6 product (https://ladsweb.modaps.eosdis.nasa.gov, last access: 1 January 2025) was used in this study. The mean of annual NPP from 2001 to 2019 was computed to represent the general status of each sample location. Then, we used the biomass ratio between roots and shoots (root-to-shoot ratio, RSR) as a proxy to estimate the amount of NPP allocated to soils as this is generally a realistic estimation for the mean long-term NPP partitioning (Gower et al., 1999). The fraction of belowground NPP can be estimated as the root mass fraction (RMF):

$$\begin{array}{}\text{(2)}& \text{RMF}={\displaystyle \frac{\text{RSR}}{\text{RSR}+\mathrm{1}}}.\end{array}$$

The RSR value at each sample location was collected from the harmonized global maps of aboveground and belowground biomass carbon density (Spawn et al., 2020). The biome-specific variabilities of RMF are shown in Fig. S5.

We used the vertical root biomass distribution to represent the belowground NPP partitioned into different soil depths (Luo et al., 2019; Xiao et al., 2022, 2023). The root distribution information was obtained from the 564 global root profiles compiled by Schenk and Jackson (2002). The logistic dose–response curve function was used to estimate the cumulative amount of root mass r(D) above a certain soil depth D (m):

$$\begin{array}{}\text{(3)}& r\left(D\right)={\displaystyle \frac{R_{\text{max}}}{\mathrm{1}+{\left(\frac{D}{D_{\mathrm{50}}}\right)}^c}},\end{array}$$

where R_max represents the total root mass, D₅₀ is the depth at which r(D) equals half of R_max, and c is a dimensionless shape parameter; regarding the latter, the reader can refer to Schenk and Jackson (2002) for details. The fraction of roots in a certain soil layer between D_u and D_l (${\text{fr}}_{D_{\mathrm{u}}-D_{\mathrm{l}}}$) can be estimated as follows:

$$\begin{array}{}\text{(4)}& {\displaystyle }\begin{array}{rl}{\text{fr}}_{D_{\mathrm{u}}-D_{\mathrm{l}}}& ={\displaystyle \frac{r_{D_{\mathrm{l}}}}{R_{\text{max}}}}-{\displaystyle \frac{r_{D_{\mathrm{u}}}}{R_{\text{max}}}}\\ & ={\displaystyle \frac{\mathrm{1}}{\mathrm{1}+{\left(\frac{D_{\mathrm{l}}}{D_{\mathrm{50}}}\right)}^c}}-{\displaystyle \frac{\mathrm{1}}{\mathrm{1}+{\left(\frac{D_{\mathrm{u}}}{D_{\mathrm{50}}}\right)}^c}}.\end{array}\end{array}$$

This root profile dataset has been used to analyse the belowground NPP allocation in several previous studies (Luo et al., 2019; Shi et al., 2021). In our study, we extended these data by using the latest Root Systems of Individual Plants (RSIP) database and the Global Root Traits (GRooT) database, which include information regarding the maximum belowground extents of terrestrial plants (Guerrero-Ramírez et al., 2021; Tumber-Dávila et al., 2022). Assuming that the root distribution is proportional to its morphological distribution (Bardgett et al., 2014; Tumber-Dávila et al., 2022), we used the form trait of roots, the maximum rooting depth (D_max), to generate an alternative representation of the root distribution. A total of 1732 geo-referenced measurements were collected, with descriptions of rooting depth from the above two databases. The vertical root distribution can be estimated according to a commonly used asymptotic equation (Gale and Grigal, 1987; Jackson et al., 1996; Zeng, 2001):

$$\begin{array}{}\text{(5)}& \text{CRF}\left(D\right)=\mathrm{1}-{\mathit{\beta }}^D,\end{array}$$

where CRF is the cumulative root fraction from the surface to soil depth D in centimetres, and β is the estimated parameter which controls the decreasing rate of root mass with increasing soil depth (Gale and Grigal, 1987; Zeng, 2001). Finally, the fraction of roots in a certain soil layer between the upper and lower depths (${\text{fr}}_{D_{\mathrm{u}}-D_{\mathrm{l}}}$) can be quantified as follows: CRF(D_l)−CRF(D_u). An illustration of the root distribution estimation using the above approach is shown in Fig. S6.

Hence, the global distributions of all root profiles that characterize the fraction of roots in the two layers were calculated and are shown in Fig. S7. Given that root distribution is generally related to the biome, vegetation type, and soil conditions (Jackson et al., 1996; Schenk and Jackson, 2005), we apply the arithmetic mean of ${\text{fr}}_{D_{\mathrm{u}}-D_{\mathrm{l}}}$ from root profile observations within the same terrestrial ecoregions (Dinerstein et al., 2017) and soil types (FAO–UNESCO, 1990) as the soil sample (Figs. S8 and S9).

It is also necessary to consider the vertical physical transportation of organic carbon, such as through leaching and/or bioturbation. We followed the function designed in a previous model that includes vertical transport of SOC (Braakhekke et al., 2011; Koven et al., 2013; Sierra et al., 2024):

$$\begin{array}{}\text{(6)}& V\left(D\right)=\text{Diff}{\displaystyle \frac{{\partial }^{\mathrm{2}}C}{\partial D^{\mathrm{2}}}},\end{array}$$

where V is the transported SOC stock (V) in a certain layer; C is the organic carbon content, which is defined volumetrically (kg m⁻³), at the depth D; and Diff is the diffusivity, which is constant and set to be $\mathrm{1}\times {\mathrm{10}}^{-\mathrm{4}}$ m² yr⁻¹ in accordance with previous studies (Koven et al., 2013; Sierra et al., 2024). Thus, the belowground NPP (BNPP) in a certain soil layer can be estimated as follows:

$$\begin{array}{}\text{(7)}& {\text{BNPP}}_{D_{\mathrm{u}}-D_{\mathrm{l}}}=\text{NPP}\cdot \text{RMF}\cdot {\text{fr}}_{D_{\mathrm{u}}-D_{\mathrm{l}}}+V_{D_{\mathrm{u}}-D_{\mathrm{l}}}.\end{array}$$

For topsoil, carbon inputs also contain a portion of carbon from surface litterfall, which should also be considered. First, the RSR dataset was adopted to support us in obtaining the aboveground NPP (ANPP). As a decent fraction of NPP was removed as the harvest product in croplands, we performed a specific calculation procedure to estimate the aboveground carbon input in this region, as described in Sect. S2. Second, we used two databases including measurements of aboveground litterfall across the globe (Holland et al., 2015; Jia et al., 2016) to estimate the fraction of aboveground NPP (ANPP) converted into litterfall. Then, we determined the fraction of litterfall allocated as the carbon input into topsoil according to the decomposition processes described in the Community Land Model (Oleson et al., 2013). The details of the calculations related to the fraction of ANPP that was transferred as carbon input into the topsoil (denoted as fr_a) are described in Sect. S3.

Therefore, the carbon input flux between upper (D_u) and lower depths (D_l) can be estimated as follows (the biome-specific variabilities of this variable are shown in Fig. S10):

$$\begin{array}{}\text{(8)}& {\text{flux}}_{D_{\mathrm{u}}-D_{\mathrm{l}}}={\text{BNPP}}_{D_{\mathrm{u}}-D_{\mathrm{l}}}+{\text{fr}}_a\cdot \text{ANPP}.\end{array}$$

Combining the estimated SOC stock (Eq. 1) and carbon allocation belowground (Eq. 8) with the equation of apparent τ, the values of τ in the topsoil and subsoil layers were calculated at all sample locations across the globe and are shown in Fig. S11. Table S1 in the Supplement shows the details of the datasets used to calculate τ for all samples.

2.2 Geospatial mapping of τ

The machine-learning-based method was adopted to generate the spatially explicit maps of τ in the two soil layers. For the procedure of geospatial predictive mapping illustrated in Fig. S12, we used the random forest (RF) model to establish the relationship between τ and potentially related geographical variables (i.e. environmental covariates) for mapping the global distribution of τ with a quantification of uncertainty. The advantage of RF is the incorporation of randomized feature selection and training sample selection (Breiman, 2001), such that it can reduce the overfitting risk and lead to a good ability in terms of generalization. This model has been effectively applied to large datasets for geospatial mapping tasks at large scales (van den Hoogen et al., 2019; Ma et al., 2021; Poggio et al., 2021; Guo et al., 2025). There are two important user-defined parameters when using the RF model. The first is the number of covariates that are randomly selected for each tree-building process. We used the rounded-down square root of the total number of covariates as this parameter value by default. The second parameter is the number of trees to be learned in RF. We set this value to be 200 considering the fact that previous soil mapping studies have shown that a number of trees larger than 150 is sufficient to obtain stable results (Wadoux, 2019; Zhang et al., 2021). The covariates, including climate, soil physical properties, soil chemical properties, and topography, were collected (Table S2 and Figs. S13–S18). These covariate maps were overlaid and thus allowed us to predict τ at the global scale with a spatial resolution of 30 arcsec (∼1 km² at the Equator). Considering the fact that some regions covered by organic soils cannot be well represented when using our method to generate estimations, we masked these regions according to the definition of organic soil in Brady and Weil (1999).

2.3 Analysis of environmental controls

To investigate the environmental controls of the spatial distribution of τ in the two layers, the RF model was adopted based on all soil sample data to quantify the relative importance (RI) of each covariate at the global scale. We used permutation-based feature importance as the metric to assess RI (Altmann et al., 2010) (see Sect. S4).

The directional effects of important environmental factors on τ were investigated by using the partial-regression model to calculate the partial correlations of topsoil and subsoil τ with the six top-ranked environmental variables from RF modelling, as mentioned above. The analysed variables are the mean annual temperature (MAT) and precipitation (MAP), the fine-particle size fraction (CLAY+SILT), the ratio of organic carbon to total nitrogen (C:N), the cation exchange capacity of the soil (CEC), and the soil pH. The partial correlation of each influencing factor was calculated at the mean level while controlling for other factors.

The potential for non-significant effects of a variable (e.g. the global impact of precipitation) on τ in a global-level linear model does not necessarily imply that the respective variable has no influence on τ. Instead, this lack of significance may arise from its divergent effects when interacting with other factors. This has led to an interest in exploring the reasons behind such nonlinear and incoherent driving mechanisms. Considering the fact that the extrinsic climate effects on soil carbon turnover can interact with regionally intrinsic soil characteristics (Doetterl et al., 2015), we used the approach outlined in Sect. S5 to exemplify the interactive effects of MAT and MAP on topsoil and subsoil τ in response to the changes in other factors (namely interactive factors). In addition to the overall importance of covariates analysed at the large scale, the local-level importance for four categories of variables was also assessed (see Sect. S6).

2.4 Evaluation of mapping results

2.4.1 Assessment of accuracy by cross-validation

The mapping accuracy was assessed by using a 10-fold cross-validation. The total sample data were divided into 10 equally sized subsets, and samples of each biome in each subset had the same proportion as in the whole dataset. Nine folds were used as the training data to fit the model, and the prediction was validated based on the one remaining fold. This procedure was carried out 10 times, each time using a different fold for validation. The mean values of the coefficient of determination (R²) and the root mean squared error (RMSE) for all folds were computed as the final accuracy metrics for the assessment of mapping results.

2.4.2 Quantification of uncertainty

To evaluate the uncertainty of predictive maps of τ, we considered two sources of error. One is the model error, which refers to the covariates that do not fully explain the variations in τ and the error in the estimation of the model parameters. Another is the measurement error, which represents the difference between the actual and recorded values of input variables related to τ calculation. To account for the modelling uncertainty, we adopted quantile regression forests (QRFs) (Meinshausen, 2006; Hengl et al., 2018) to derive prediction intervals. QRFs first construct an RF model in the usual way, namely by developing multiple decision trees that use subsets of the training data, whereby the prediction of each tree equals the average of the observations in the end node of the tree in which the prediction point sits. The RF prediction is the average of all tree predictions. Since averaging is a linear process, the RF prediction boils down to a weighted sum over the n observations of the response variable:

$$\begin{array}{}\text{(9)}& \widehat{y}\left(x\right)=\sum _{i=\mathrm{1}}^n{w}_i\left(x\right)\cdot y_i.\end{array}$$

In QRFs, the weights w_i(x) are used to estimate the cumulative distribution $F\left(y\mathrm{|}x\right)=P(Y\le y\mathrm{|}x)$ of the response variable Y, given the covariate data x, as follows:

$$\begin{array}{}\text{(10)}& \widehat{F}\left(y\mathrm{|}x\right)=\sum _{i=\mathrm{1}}^n{w}_i\left(x\right)\cdot {\mathrm{1}}_{\mathit{\{}{y}_i\le y\mathit{\}}},\end{array}$$

where ${\mathrm{1}}_{\mathit{\{}{y}_i\le y\mathit{\}}}$ is the indicator function (i.e. it is 1 if the condition is true and 0 otherwise). Any quantile q of the distribution can then be derived by iterating towards the value of y for which $\widehat{F}\left(y\mathrm{|}x\right)=q$.

We used the R package ranger (Wright and Ziegler, 2017) with the function quantreg to build the QRF models. Using this function, the prediction values at each location can be obtained, and the 0.05th quantile (q_0.05), 0.50th quantile (q_0.50), and 0.95th quantile (q_0.95) can also be computed to derive the lower limit, median, and upper limit of a symmetric 90 % prediction interval. This interval has also been adopted for uncertainty assessment in GlobalSoilMap specifications (Arrouays et al., 2014) and the SoilGrids product (Poggio et al., 2021).

The uncertainty from the input data in the calculation of τ was also considered. When calculating τ for each sample, four input variables (SOCS, NPP, RMF, and ${\text{fr}}_{D_{\mathrm{u}}-D_{\mathrm{l}}}$) introduced uncertainty due to errors in the soil measurements and in the estimation of carbon allocation belowground. Errors in these input variables will propagate into the output of τ estimation. To incorporate the error propagation from these inputs into the uncertainty evaluation, the standard deviation (SD) of each input at each sample location needs to be quantified. Here, as the calculation of τ is a simple arithmetic function, the SD of the estimated τ for each sample can be calculated as follows, ignoring the cross-correlation among those inputs (Heuvelink, 1998):

$$\begin{array}{}\text{(11)}& {\text{SD}}_f=\sqrt{\begin{array}{c}{\left(\frac{\partial f}{\partial I_{\mathrm{1}}}\right)}^{\mathrm{2}}{\text{SD}}_{I_{\mathrm{1}}}^{\mathrm{2}}+{\left(\frac{\partial f}{\partial I_{\mathrm{2}}}\right)}^{\mathrm{2}}{\text{SD}}_{I_{\mathrm{2}}}^{\mathrm{2}}\\ +\mathrm{\dots }+{\left(\frac{\partial f}{\partial I_m}\right)}^{\mathrm{2}}{\text{SD}}_{I_m}^{\mathrm{2}}\end{array}},\end{array}$$

where SD_f represents the SD of the output value calculated by the function f, and I₁, I₂, and I_m represent the first, second, and mth input variables, respectively. For RMF, we directly used the uncertainty maps provided by Spawn et al. (2020). The SD of SOCS was obtained from Poggio et al. (2021) when SOCS was supplemented by extracting values from the SoilGrids maps. For quantifying the SD of ${\text{fr}}_{D_{\mathrm{u}}-D_{\mathrm{l}}}$, the SD of all selected root profile observation values for each soil sample site was adopted to represent its uncertainty. For NPP, most MODIS NPP values are within the mean ±1 SD of the respective values observed from the flux tower (Reichstein et al., 2002), based on which we define the uncertainty of NPP values. The relationships between the estimated values of τ and their corresponding SD for all samples are shown in Fig. S19. After quantifying the SD of the estimated τ at each sample location, a Monte Carlo approach was adopted to incorporate the SD into the estimated τ for each sample. That is, the value of τ at each sample location was randomly drawn 100 times from a normal distribution given the known mean and SD. Then, all generated samples were used to fit a QRF model and to produce the uncertainty maps.

To visualize the spatial distribution of the prediction uncertainty, we calculated the prediction interval ratio (PIR), defined as the ratio of the range between lower and upper limits over the median:

$$\begin{array}{}\text{(12)}& \text{PIR}={\displaystyle \frac{q_{\mathrm{0.95}}-q_{\mathrm{0.05}}}{q_{\mathrm{0.50}}}}.\end{array}$$

A lower PIR indicates higher confidence in the predictions, while a higher PIR suggests greater uncertainty.

We further assessed the quality of the estimated prediction uncertainty through an accuracy plot approach by calculating the prediction interval coverage probability (PICP) (Goovaerts, 2001). PICP is a metric used for evaluating the reliability of uncertainty quantifications by measuring the proportion of observed values that fall within a given prediction interval. It can effectively evaluate whether the probability assigned to the prediction intervals is equal to the frequency of empirical test data within the prediction intervals (Malone et al., 2011; Schmidinger and Heuvelink, 2023).

3 Results and discussion

3.1 Accuracy assessment

The cross-validation demonstrated that the machine learning models can effectively capture a substantial proportion of τ variations, achieving R² values of 0.87 and 0.70 and RMSE values of 14.60 and 175.05 for the topsoil and subsoil τ predictions, respectively (Fig. 2). More accuracy metrics are shown in Table S3. Relatively higher model performance was found in boreal and tundra areas for topsoil and in grasslands and shrublands in warm regions for subsoil, while low model performance was found in temperate forests, croplands and wetland areas for the subsoil layer (Figs. S20 and S21). This is probably because soils in these regions constitute some of the most diverse soil landscapes compared to the pedogenetically similar and more climate-driven soils in high-latitude regions and the more weathered homogenous soils dominating many tropical lowland areas.

Figure 2 Validation plots for predictions of soil organic carbon turnover time (τ, years) in topsoil (a) and subsoil (b) layers. Predictions were generated using random forest regression models. Grey dots are the estimated τ values based on observations for soil profiles and the predicted values for all validation samples using 10-fold cross-validation. Coloured points represent the biome-level validation plot. Error bars show 95 % confidence intervals of τ in each biome. Black lines indicate the regression lines between predicted and measured values. Axes are log₁₀-transformed to account for high skewness of τ.

Table 1 Global- and biome-level statistics of the soil organic carbon turnover time (τ, years) in topsoil (0–0.3 m) and subsoil (0.3–1 m) layers.

Note that P 2.5 and P 97.5 represent the range of τ between the 2.5th and 97.5th percentiles in each biome based on the aggregation of all estimated τ at the pixel level.

3.2 Global distributions of τ in topsoil and subsoil

The spatially continuous maps of global topsoil and subsoil τ are shown in Fig. 3. On average, the global τ in the topsoil was 45 years, ranging from 11 years (2.5th percentile) to 137 years (97.5th percentile). The average subsoil τ was 385 years, ranging from 20 to 3485 years (Table 1). On a global scale, subsoil SOC turnover time was over 8 times longer than that in topsoil. The values of τ generally increased from low to high latitudes, and this latitudinal pattern was more pronounced in the Northern Hemisphere than in the Southern Hemisphere, and it was more evident in the subsoil than in the topsoil layer (Fig. 3b and d). In tropical forests, the average turnover times were shortest, with approximately 25 and 74 years in the topsoil and subsoil, respectively. The turnover times in temperate, desert, and cropland areas for topsoil and subsoil layers were around 29–43 years and 78–170 years, respectively. The longest turnover times were found for tundra regions, with averages of longer than 90 and 1900 years for the topsoil and subsoil layers, respectively. Boreal forests and wetlands also show long carbon turnover times, with average values of over 50 and 500 years, respectively. As such, the differences in average turnover times between the warmest and coldest biomes were more than 60 and 1800 years for topsoil and subsoil layers, respectively.

Figure 3 Global patterns of topsoil and subsoil organic carbon turnover times (τ, years). Global distributions of τ in topsoil (0–0.3 m) (a) and subsoil (0.3–1 m) (c) layers. The predicted maps, with a spatial resolution of 30 arcsec (∼1 km at the Equator), were generated with τ–environment relationships using a machine learning model trained on global soil profile observations and their environmental covariates. (b, d) Latitudinal patterns of topsoil and subsoil τ. Orange and blue lines represent the average τ in topsoil and subsoil, respectively, with latitude. The shaded grey areas represent the variations in τ between the 2.5th and 97.5th percentiles with latitude. (e, f) Average τ of the two layers in different main biomes. Error bars show the 95 % confidence intervals of the spatial predictions within each biome.

The uncertainty maps quantified for the topsoil and subsoil are shown in Fig. 4. Uncertainties in subsoil τ predictions were generally higher than those in topsoil τ predictions, likely due to the fewer horizontal observations in the deep layers. Wide prediction intervals were mostly found in areas with low sampling density, such as deserts and permafrost regions. According to the PICP calculation results, the accuracy plots (Fig. S22) show that the proportion of observed τ values in the validation set as covered in a certain prediction interval approximately equal the size of that probability interval, which validates the unbiased quantification of uncertainty.

Figure 4 Uncertainty maps of predicted soil organic carbon turnover time (τ) in topsoil (a) and subsoil (b) layers. The uncertainty was quantified by using quantile regression forests. The values shown in the maps represent the prediction interval ratio (PIR), which is the ratio of the range between the lower and upper limits (90 % prediction interval bounded by the lower (0.05th quantile) and upper (0.95th quantile) limits) over the median (0.50th quantile) of the predictions.

3.3 Environmental controls of topsoil and subsoil τ

3.3.1 Global- and biome-level analyses

The relative importance and directional effects of climate, soil (divided into physical and chemical properties), and topographic factors in relation to τ variation at global and biome scales were analysed (Figs. 5, S23, and S24). The results suggest scale and depth dependency in the drivers of soil turnover times. While climate was found to be the primary driver (explaining nearly half of the total variation) of topsoil τ at the global scale, the intrinsic soil characteristics (61 %) were more important than climate (32 %) for the subsoil τ (Figs. S23a and S24a). This supports recent studies suggesting that soil turnover in deeper layers is less controlled by climatic factors (Luo et al., 2019; Chen et al., 2021; Han et al., 2022) and more controlled by soil properties (Mathieu et al., 2015; Luo et al., 2019). Importantly, by separating soil factors into soil physical (represented by the fine-particle size fraction (CLAY+SILT)) and chemical categories (including C:N, CEC, and soil pH), our study revealed that the topsoil τ is more influenced by soil physical properties than by chemical properties (33 % versus 15 %). However, for subsoil τ, soil physical properties are less important than soil chemical properties (27 % versus 39 %). These results support the finding that, in addition to the direct control of climate, the pedologic traits and geochemistry of soils exert an equivalent or even stronger control on τ at the global scale (Davidson and Janssens, 2006; Doetterl et al., 2015, 2016).

Figure 5 Relationships between environmental variables and soil organic carbon turnover time (τ) in topsoil and subsoil layers. Effects of six important variables on τ are shown for the global and biome levels. Values in the bar plots represent coefficients between τ and each variable from the partial-regression model, where positive or negative values indicate a positive or negative effect on τ, respectively. MAT denotes mean annual temperature, MAP denotes mean annual precipitation, CLAY+SILT denotes fine-particle size fraction (the sum of clay and silt content), C:N denotes ratio of organic carbon to total nitrogen, CEC denotes cation exchange capacity of the soil, and pH denotes soil pH. * P<0.001.

In different biomes and at different soil depths, the effects of MAT and MAP show different magnitudes or even opposite directions (Fig. 5a and b). Although a general negative effect of MAT on τ has been found in a previous meta-analysis (Chen et al., 2013), our results show that such an effect was reduced and even changed to be positive in tropical regions (Fig. 5a). This can be attributed to the contrasting relationships of NPP with temperature between tropical and extratropical regions, stemming from their distinct limiting factors affecting plant growth (Fig. S25) (Chapin et al., 2011; Slot and Winter, 2016). In addition, the effect of MAP on τ was not significantly linear at the global scale. This is inconsistent with a previously detected negative relationship (Schimel et al., 1994; Chen et al., 2013) but supports a recent global study that detected a nonlinear relationship between τ and hydrometeorological conditions (Fan et al., 2022). Such previous discrepancies can be explained by contrasting the effects of MAP across different biomes (Fig. 5b). While warm forests and croplands exhibit significant positive effects, notable negative effects are shown in tropical savannahs and grasslands and in arid regions. Additionally, it is shown that MAP negatively impacts subsoil τ in both boreal forests and temperate grasslands and shrublands, but this effect was not significant for topsoil. These τ–climate patterns highlight the nonlinear effects of temperature and water availability on τ globally and the different driving mechanisms across biomes.

Our analyses further revealed that soil physical and chemical properties have comparable or greater effects on τ than climatic effects at global and biome scales (Fig. 5c–f). For soil physical properties, as represented by CLAY+SILT, the fine-particle size fraction in soils was positively related to τ in general (Fig. 5c). This is mainly because soils dominated by finer particles tend to stabilize SOC by means of physical and organo-mineral stabilization mechanisms, while the well-drained sandy soils generally provide limited protection of SOC against microbial decomposition (Krull et al., 2003; Lützow et al., 2006; Cotrufo et al., 2019). The relatively higher effects of CLAY+SILT on subsoil τ than on topsoil τ (except in croplands) further indicate that mineral-associated and physically protected SOC in deep layers may play a crucial role in stabilizing SOC stocks and, consequently, may decrease climate sensitivity (Gillabel et al., 2010; Qin et al., 2019).

Among the effects of soil chemical properties, the C:N ratio plays a key role in regulating subsoil τ, ranking as the most influential soil factor (Fig. S24a). The positive impact of C:N was observed in all biomes for subsoil and specifically in boreal forests, temperate grasslands, deserts, and tundra for topsoil (Fig. 5d). However, the effects of C:N were not found to be significant for topsoil τ, and it is even shown that there is a negative impact in tropical and cropland regions. This finding implies that, in general, nutrient limitations decrease soil respiration and extend carbon residence times by favouring the establishment of slow-decomposing organisms or reducing organic matter quality (Crowther et al., 2019b; Li et al., 2012; Zhou et al., 2007). Whilst this effect is pronounced in the subsoil layer, it cannot be straightforwardly generalized to the topsoil in the tropics and in human-impacted regions. The factor of CEC showed a generally positive impact on τ in forests and croplands, but its effect was weaker in other biomes (Fig. 5e). Soil pH shows a non-significant linear correlation with τ globally, yet it exhibits divergent influences across different biomes and soil depths (Fig. 5f). It has a positive impact on τ in warm forests and on subsoil τ in temperate grasslands and croplands. Conversely, in boreal forests, tropical grasslands, and arid regions, it more negatively impacts τ in topsoil.

In addition to the climate and soil factors, topography can influence τ by affecting soil moisture distribution, drainage conditions, and permafrost dynamics. Our analysis (Figs. S23h and S24h) indicates that topographic variables contribute significantly to τ variability in tundra regions. This might be related to previous studies that found that terrain-driven hydrological differences strongly influence permafrost carbon stability (Schuur et al., 2008; Mishra et al., 2021). This also aligns with research showing that microtopographic variations, such as polygonal tundra structures, create distinct carbon accumulation and decomposition regimes, further influencing SOC turnover times (Jorgenson et al., 2010; Liljedahl et al., 2016).

Figure 6 The climatic effects on topsoil and subsoil organic carbon turnover time (τ) interacting with other factors. The lines show the mean annual temperature (MAT) effect and the mean annual precipitation (MAP) effect on τ in response to another climatic variable and four important soil physio-chemical (including the fine-particle size fraction (CLAY+SILT), carbon-to-nitrogen ratio (C:N), cation exchange capacity (CEC), and soil pH) variables. The variables on the x axis can be considered to be interactive factors, which influence the MAT and MAP effects on τ. The values on the lines above or below the dashed horizontal lines indicate the positive or negative climatic effects on topsoil or subsoil τ (the MAT and MAP effects were calculated by means of the partial-regression coefficient between MAT and/or MAP and τ) in response to the corresponding values of interactive factors. The shaded areas on the lines represent confidence intervals for the results of partial-correlation analysis.

3.3.2 Interactive effects of climatic and edaphic factors on τ

The detected relationships between environmental factors and τ prompted us to delve deeper into the nonlinear driving mechanisms. We aimed to uncover potential interactions among the primary climatic and edaphic drivers that might be masked by analysing the directional effects of only a single variable. We illustrate how fluctuations in one climatic variable (MAT or MAP) related to changes in another climatic variable (Fig. 6a and f) and four key soil properties (Fig. 6b–e and g–j). We found that the negative effect of MAT on τ was magnified with decreasing MAP (when MAP <1500 mm). Meanwhile, it was also shown that the MAT effect on τ changed to be smaller and could even be positive for subsoil with increasing MAP (when MAP>2000 mm) (Fig. 6a). Notably, the MAP effect transitioned from negative to positive with increasing MAT in colder environments. This effect then diminished after a critical point around 15 °C (Fig. 6f).

Beyond the interplay of climatic factors, soil property values significantly modulated the effects of MAT and MAP on τ. For instance, as CLAY+SILT content decreased, the negative impact of MAT was pronounced (Fig. 6b). The MAP effect was also accentuated in sandy soils but was restricted to the topsoil layer (Fig. 6g). This result suggests that the finer soil texture may reduce carbon turnover sensitivity to temperature or water fluctuations (Krull et al., 2003). Among the interactions of soil chemical properties with climatic impacts on τ, it is evident that extremely low or high C:N ratios in soils lead to smaller climatic effects on τ (Fig. 6c and h), suggesting that the sensitivity of τ to climate is strongly influenced by whether the C:N ratio falls within an optimal range. The negative impacts of temperature were intensified with rising CEC (when CEC<40 cmol_c kg⁻¹) (Fig. 6d), mainly because higher CEC typically enhances the nutrient-supplying capacity of soils, thus quickening carbon turnover in warm conditions (Crowther et al., 2019a). The different responses of topsoil and subsoil temperature sensitivity to CEC (Fig. 6d) may result from interactive and combined influencing factors. Higher CEC enhances nutrient availability for microbes and promotes mineral–organic compound formation, potentially increasing SOC temperature sensitivity. Concurrently, elevated CEC can promote plant productivity, leading to greater carbon inputs and influencing NPP sensitivity to temperature. The observed patterns are likely to reflect the combined effects of these processes. The influence of soil pH on MAT effects on τ demonstrated a strong nonlinear trend. A more pronounced negative MAT effect was found when soil pH was neutral for topsoil and mildly alkaline for subsoil (Fig. 6e). This is related to microorganisms having a physiologic optimum pH (Rousk et al., 2010; Don et al., 2017) and results in the warmer temperature significantly accelerating SOC turnover rates when pH is within an optimum range (Frostegård et al., 2022; Xiang et al., 2023). Nevertheless, there is currently a need to further explore the reasons for why these pH ranges are different with soil depth, as well as for their species specificity (Bahram et al., 2018). The impact of MAP appears to be more positive in acidic soils but becomes negative when soil pH exceeds 8 and 9 for topsoil and subsoil, respectively (Fig. 6j). The positive MAP effect on turnover time under acidic conditions could be partly explained by the enhanced effect of increased moisture on soil weathering processes (Porras et al., 2017). Under alkaline conditions, which are mostly seen in dry climates (Slessarev et al., 2016), the negative MAP effect may arise from the intensified microbial transformation of plant-derived organic matter as these soils provide pH and moisture that are favourable to microorganisms (Yang et al., 2022).

Table 2 Statistical summary of the local-scale analytical results of variable (classified into four categories) importance in terms of influence on topsoil and subsoil organic carbon turnover time (τ) at the local scale.

Figure 7 Local-level analysis of dominating factors of topsoil and subsoil organic carbon turnover time (τ). Global maps of dominant factors controlling τ in topsoil (a) and subsoil (b) at the local scale. Maps were created by interpolating all results of factor importance from modelling using the sample data within each local area stratified by the ecoregion map. The lighter colours show that the reliability of the local result is low, with a smaller number of available samples for analysis. Histograms of the local-level percentage contribution of each variable category in controlling topsoil (c) and subsoil (d) τ across all local areas. Dashed lines represent the mean values.

3.3.3 Distributions of dominating factors of τ at local scales across the world

We further analysed the effects of these covariates by generating τ–environment models using data from sampling locations within different local areas to assess the deviation of regional patterns from those at the large scale. Overall, we find that the dominant factors controlling τ vary across different local areas (Fig. 7a and b). Across approximately 41 % of the land surface area, climatic factors are the predominant controller of topsoil τ (Table 2). However, for subsoil τ, climatic factors predominantly mediate the variation in τ in only 29 % of the global area. It is noteworthy that, in 57 % of the global area, soil characteristics were the dominant driver of subsoil τ. Overall, soil chemical properties had a larger effect than soil physical properties. This is particularly notable in parts of the tropical forests of the Amazon; the broadleaf forests spanning from eastern Indian to Southeast Asia; most parts of North America; central Europe; and the grasslands, deserts, and xeric shrublands scattered across Africa and Australia (Fig. 7b). This phenomenon is related to previous studies suggesting that nutrient availability for soils and vegetation greatly affects τ (Cleveland and Townsend, 2006; Carvalhais et al., 2014). The contributions of four categories of environmental factors in the local-scale analysis differed less in subsoil compared to in topsoil (Fig. 7c and d and Table 2).

Figure 8 Comparisons of soil organic carbon turnover times (τ) from our data-derived estimates and from Earth system models (ESMs). (a) Biases in τ represented by ratios of data-derived to model-derived values across the global grid cells. Black points indicate locations where τ values from ESMs are outside of the 0.05th and 0.95th quantiles of the prediction uncertainty in our data-derived estimations. (b, c) Comparisons of associations of data- and model-derived (ESMs) τ to mean annual temperature (MAT) (b) and mean annual precipitation (MAP) (c) along latitudinal gradients. The y axis shows partial correlations (r) of τ with two climate factors, controlling for precipitation when calculating the correlation between τ and temperature (and vice versa).

3.4 Comparison with τ from Earth system models

Our observation-based τ estimates can help constrain ESM simulations and improve predictions of current and future carbon cycle dynamics. We analysed historical simulation outputs of selected ESMs from CMIP6 (Table S4). While the comparison of our empirical estimates with τ estimates from ESM outputs based on CMIP6 showed broad agreement in terms of the spatial patterns of τ (Pearson's r=0.53, P<0.001), it highlighted that ESMs are likely to underestimate τ across the majority of the globe (Fig. 8a). The soil carbon turnover times estimated from ESMs were, on average, more than 2 times shorter than our data-derived (using ground-sourced samples integrated with remote sensing observations) τ estimates (Table S4). This discrepancy was particularly pronounced in tropical forests, grasslands, and tundra areas (Table S5). In ∼30 % of the global land grids, the τ estimates from ESMs did not fall into the 90 % estimation intervals of our data-driven estimates, and ∼92 % of global land area showed a mean underestimation bias (ESM-derived τ < data-derived τ) (Fig. 8a, Table S4).

The discrepancies between τ from observations and from ESMs are also reflected in their associations with climate variables. The ESMs predicted stronger correlations with temperature and precipitation compared to our τ estimates (Fig. 8b and c). Additionally, there were notable differences in terms of the climate correlations between topsoil and subsoil, which ESMs do not yet capture. These findings emphasize that other factors, such as soil physico-chemical properties and soil depth, need to be accounted for in models to accurately project global soil carbon dynamics.

3.5 Implications and future perspective

Our quantitative maps of SOC turnover times in topsoil and subsoil layers based on an extensive soil profile dataset represent a key step towards a better understanding of global soil carbon stocks and dynamics. Our research illustrates how global variations in τ in both topsoil and subsoil layers are influenced by the interplay of climatic and edaphic factors and demonstrates their pronounced heterogeneity within and across biogeographic zones. These analyses revealed complex interactions between temperature, water availability, and soil physio-chemical properties, thereby enriching our comprehension of the complex mechanisms driving spatial variability and nonlinearity in τ–environment relationships from global to local scales. The distinct factors driving τ in topsoil versus subsoil underscore the importance of incorporating soil depth when assessing large-scale τ patterns.

The findings in this study also have the potential to improve the parameterization and future projections of ESMs by integrating more accurate global τ data. Our results support several previous studies that also identified an underestimation of SOC turnover times in ESMs through ¹⁴C observations (He et al., 2016; Shi et al., 2020). This discrepancy between data- and model-derived τ emphasizes the need to incorporate more detailed edaphic- and climate-dependent and depth-specific τ estimates into biogeochemical models to enhance predictive accuracy. Furthermore, there is evidence that the greatest exchange of soil carbon with the atmosphere occurs through relatively small and fast soil carbon pools on a short timescale, which can result in a “leaky-sink” response when carbon input is elevated (Bradford, 2017; van Groenigen et al., 2017). This response may conceal the longer turnover times in inert pools that are not precisely captured by models. Adjusting turnover rates and carbon transfer parameters in ESMs to align with the longer τ values we report here may extend the turnover time of “slow” or “passive” pools in models (especially for deep layers) that constitute the majority of soil carbon (Torn et al., 2009; He et al., 2016).

To improve the accuracy of carbon cycle models and to enhance the projections of future soil carbon sequestration rates and magnitudes, it is essential to incorporate edaphic- and climate-dependent and depth-resolved estimates of turnover times into these models. These insights across various climate zones, biomes, terrains, and soil properties contribute to reducing uncertainties related to context-dependent effects governing soil carbon stocks and dynamics, thus helping to inform strategies that enhance sustainable soil management and mitigate the impacts of climate change.

4 Data availability

The global maps of topsoil and subsoil organic carbon turnover times are available online at https://doi.org/10.5281/zenodo.14560239 (Zhang, 2025a).

5 Code availability

The code used for this study is publicly available at https://github.com/leizhang-geo/global_soil_carbon_turnover_time.git and is permanently deposited in a Zenodo repository at https://doi.org/10.5281/zenodo.15636310 (Zhang, 2025b).

6 Conclusions

This study provides a comprehensive assessment of global apparent SOC turnover times in both topsoil and subsoil layers. By integrating state-of-the-art datasets of soil profiles, plant roots, and satellite observations, we employed machine learning models to produce the spatially explicit maps of τ with quantified uncertainties. The results reveal pronounced spatial heterogeneity and nonlinearity in τ–environment relationships within and among different biomes and climatic zones. Our findings demonstrate the context-dependent effects of temperature and water availability on τ, which vary depending on soil attributes, mainly including soil texture, carbon-to-nitrogen ratio, cation exchange capacity, and soil pH. Considering the potentially large differences in the driving factors of τ from large to small scales, we further mapped the dominating factors of τ in two soil layers at local scales. Overall, this study synthesizes the multiple observation-based datasets that are currently available and provides new global maps of topsoil and subsoil organic carbon turnover times. The dataset and new insights of this study are expected to serve as a foundation for benchmarking biogeochemical models and supporting effective carbon management.