Articles | Volume 16, issue 4
Research article
27 Feb 2019
Research article |  | 27 Feb 2019

Evaluating the simulated mean soil carbon transit times by Earth system models using observations

Jing Wang, Jianyang Xia, Xuhui Zhou, Kun Huang, Jian Zhou, Yuanyuan Huang, Lifen Jiang, Xia Xu, Junyi Liang, Ying-Ping Wang, Xiaoli Cheng, and Yiqi Luo

One known bias in current Earth system models (ESMs) is the underestimation of global mean soil carbon (C) transit time (τsoil), which quantifies the age of the C atoms at the time they leave the soil. However, it remains unclear where such underestimations are located globally. Here, we constructed a global database of measured τsoil across 187 sites to evaluate results from 12 ESMs. The observations showed that the estimated τsoil was dramatically shorter from the soil incubation studies in the laboratory environment (median = 4 years; interquartile range = 1 to 25 years) than that derived from field in situ measurements (31; 5 to 84 years) with shifts in stable isotopic C (13C) or the stock-over-flux approach. In comparison with the field observations, the multi-model ensemble simulated a shorter median (19 years) and a smaller spatial variation (6 to 29 years) of τsoil across the same site locations. We then found a significant and negative linear correlation between the in situ measured τsoil and mean annual air temperature. The underestimations of modeled τsoil are mainly located in cold and dry biomes, especially tundra and desert. Furthermore, we showed that one ESM (i.e., CESM) has improved its τsoil estimate by incorporation of the soil vertical profile. These findings indicate that the spatial variation of τsoil is a useful benchmark for ESMs, and we recommend more observations and modeling efforts on soil C dynamics in regions limited by temperature and moisture.

1 Introduction

Carbon (C) cycle feedback to climate change is highly uncertain in current Earth system models (ESMs) (Friedlingstein et al., 2006; Bernstein et al., 2008; Ciais et al., 2013; Bradford et al., 2016), which largely stems from their diverse simulations of C exchanges among the atmosphere, vegetation, and soil (Luo et al., 2016; Smith et al., 2016; Mishra et al., 2017). Soil organic carbon (SOC) represents the largest terrestrial carbon pool, which stores at least 3 times as much as the atmospheric and vegetation C reservoirs (Parry et al., 2007; Bloom et al., 2016). However, a 5- to 6-fold difference in soil C stocks among ESMs or offline global land surface models has been found (Todd-Brown et al., 2013; Luo et al., 2016). It is difficult to reduce or even diagnose this uncertainty, as many processes collectively affect the time of C atoms' transit in the soil system (i.e., transit time, τsoil) (Sierra et al., 2017; Spohn and Sierra, 2018). Some recent attempts at evaluating and diagnosing the modeled SOC in ESMs have shown significant simulation uncertainties in the τsoil (Todd-Brown et al., 2013; Carvalhais et al., 2014; He et al., 2016; Koven et al., 2017). For example, there is a 4-fold difference in the simulated τsoil among the ESMs from the Coupled Model Intercomparison Project Phase 5 (CMIP5) (Todd-Brown et al., 2013). A recent data-driven analysis has suggested that the current ESMs have substantially underestimated the τsoil by 16–17 times at the global scale (He et al., 2016). Therefore, identifying the locations of such underestimations is critical to improve the predictive ability of ESMs for the terrestrial C cycle, and the construction of a benchmarking database of available observations is urgently needed (Koven et al., 2017).

The terms of transit time, turnover time, and age of soil C have been muddled in diagnosing the models (Sierra et al., 2017). The diagnostic times derived from observational data are based on the different assumptions and mainly derived from four approaches. The first approach is commonly defined as “turnover time”, calculated by the division of SOC stock by C fluxes such as net primary productivity (NPP) or heterotrophic respiration (Rh). It assumes the soil system is a time-invariant linear system in a steady state (Bolin et al., 1973; Sanderman et al., 2003; Six and Jastrow, 2012). The second approach is based on the shifts in stable isotopic C (13C) after successive changes in C3C4 vegetation, together with additional information from the disturbed and undisturbed soils (Balesdent et al., 1987; Zhang et al., 2015). The third approach is based on simulating soil C dynamics with linear models by assimilating the observational data from laboratory incubations of soil samples (Xu et al., 2016). The last approach derives the weighted inverse of the first-order cycling rate by fitting a one- or multiple-pool linear model to field observations of radiocarbon (14C) (Trumbore et al., 1993; Fröberg et al., 2011). The diagnostic times derived from the first three approaches indicate the transit times, which are the mean ages of C atoms leaving the carbon pools during a certain time (Rasmussen et al., 2016). Lu et al. (2018) has evaluated the deviation between C transit and turnover times with the CABLE model. Their results have shown that the global latitudinal pattern of C transit and turnover times is consistent under a steady-state assumption and autonomous conditions except for 8 % of divergence in the northern high latitudes (>60 N). However, the diagnostic time calculated by the radiocarbon signal indicates the average age of C atoms stored in the C pools. Although radiocarbon has been widely used to quantify the age or transit time of soil C, its validity has been challenged by some recent theoretical analyses (Sierra et al., 2017; Metzler et al., 2018). Rasmussen et al. (2016) has marked off the transit time and mean system age in a mathematic way and further applied it in the Carnegie–Ames–Stanford approach (CASA) model. Also, the methodological uncertainty is large, especially when these approaches are applied to estimate the τsoil of different soil fractions (Feng et al., 2016). Thus, this study mainly collects the τsoil values from the approaches of stock over flux, 13C changes, and lab incubations in further analyses.

In this study, we first construct a database from the literatures which reported τsoil (Fig. 1a, Supplement on Text S1). Then, the database is used to evaluate the simulated τsoil by the ESMs in CMIP5. The SOC τsoil values were calculated under a homogenous one-pool assumption at steady state for all studies. Data from observations and the CMIP5 ensemble were then used to calculate the τsoil based on both one-pool and three-pool models. Many ESMs, e.g., CESM, have released new versions in recent years, so we also evaluate whether the simulated τsoil has been improved. In the case of CESM, one of its major developments in soil C cycling is the vertically resolved soil biogeochemical scheme (Koven et al., 2013). Thus, we employ a matrix approach developed by Huang et al. (2017) to examine the impact of the vertically resolved soil biogeochemical scheme on the τsoil simulated by CESM.

2 Materials and methods

2.1 A global database of site-level τsoil

We collected the literatures that reported the τsoil based on measurements (Text S1 in the Supplement): (1) δ13C shifts after successive changes in C3-C4 vegetation, (2) measurements of CO2 production in laboratory SOC incubation over at least 7 months, and (3) simultaneous measurements of SOC stock and heterotrophic respiration (stock over flux). We constructed a database containing the measured τsoil from 187 sites across the globe (Fig. 1). Based on the homogenous assumption, the soil system is a time-invariant linear system at the steady state. The τsoil derived from this database is under a one-pool assumption. The information of climate (e.g., mean annual temperature and precipitation) was also collected from the literature or extracted from the WorldClim database version 1.4 (, last access: 1 June 2019, Hijmans et al., 2005) if literature was not available. The WorldClim dataset provided a set of free global climate data for ecological modeling and Geographic Information System analysis with a spatial resolution of 0.86 km2 (Hutchinson et al., 2004). We extracted the mean temperature and precipitation by averaging the monthly climate data over 1990–2000 for those observational sites with missing climate information. The classes of biomes were processed to match the seven biome classifications adopted by the MODIS land cover product MCD12C1 (NASA LP DAAC, 2008; Friedl et al., 2010) and Todd-Brown et al. (2013) (Fig. S1): (1) tropical forest including evergreen broadleaf forest between 25 N and 25 S; (2) temperate forest including deciduous broadleaf, evergreen broadleaf outside of 25 N and 25 S, and mixed forest south of 50 N; (3) boreal forest including evergreen needleleaf forest, deciduous needleleaf forest, and mixed forest north of 50 N; (4) grassland and shrubland including woody savanna south of 50 N and savanna and grasslands south of 55 N; (5) deserts and savanna including barren or sparsely vegetated open shrubland south of 55 N and closed shrubland south of 50 N; (6) tundra; and (7) croplands. Other land cover types like permanent wetland, urban, and bare land were not included in this study.

2.2 Outputs of Earth system models from CMIP5

The historical simulation outputs of 12 ESMs participating in CMIP5 from 1850 to 1860 (, last access: 11 January 2016) were analyzed in this study (Table S1). For each model, the SOC, litter C, NPP, and heterotrophic respiration (Rh) were extracted from the outputs in historical simulations (cSoil, cLitter, npp, and rh, respectively, from the CMIP5 variable list). The litter and soil carbon were summed as the bulk soil carbon stock. Among the 12 models, only the inmcm4 model did not output NPP, so we calculated it as gross primary production minus autotrophic respiration. Due to the diverse spatial resolutions among the models, we aggregated the results of different models to 1×1 with the nearest interpolation method (Fig. S2). The τsoil of SOC was calculated as the ratio of carbon stock over flux (NPP or Rh):

(1) τ soil = SOC flux .

Figure 1Spatial distributions of observational sites for estimates of SOC transit time (τsoil, year). (a) The site locations of measurements with different approaches. (b) Probability density functions of τsoil measured with different approaches. Note that the left axis is for the 13C and stock-over-flux approaches, and the right axis is for laboratory incubation studies.


2.3 Estimating the SOC τsoil with a three-pool model

To examine whether the major findings of this data–model comparison are affected by the one-pool homogenous assumption, we fitted a three-pool model with observational data and model ensemble outputs at the biome level. In this study, a three-pool C model consisted of fast, slow, and passive pools and carbon transfers among three pools (Fig. S3a). This model shares the same framework with the CENTURY and the terrestrial ecosystem models (Bolker et al., 1998; Liang et al., 2015). The dynamics of soil carbon pools follow first-order differential kinetics. The total C stocks and CO2 efflux from observations and the CMIP5 ensemble were separated into pool-specific decomposition rates by the deconvolution analysis (Fig. S3a, Liang et al., 2015). We assumed the total soil carbon input equals total soil respiration at the steady state.

Based on the theoretical analysis, the dynamics of the three-pool model can be mathematically described by the matrix equation (Luo et al., 2003; Xia et al., 2013) as

(2) d C t d t = I ( t ) - AKC t ,

where the matrix C(t)=(C1(t), C2(t), C3(t))T is used to describe soil carbon pool sizes. A is a matrix given by

(3) A = - 1 f 12 f 13 f 21 - 1 0 f 31 f 32 - 1 .

The elements fij are carbon transfer coefficients, indicating the fractions of the carbon entering the ith (row) pool from the jth (column) pool. K is a 3×3 diagonal matrix indicating the decomposition rates (the amounts of carbon per unit mass leaving each of the pools per year). The matrix of K is given by K=diag (k1, k2, k3).

The parameters in the three-pool model were estimated based on Bayesian probabilistic inversion (Eq. 4). The posterior probability density function P(θ|Z) of model parameters (θ) can be represented by the prior probability density function (P(θ)) and a likelihood function (P(Zlθ)) (Liang et al., 2015; Xu et al., 2016). The likelihood function was calculated by the minimum error between observed and modeled values with Eq. (5). In this study, we adopted the prior ranges of model parameters from Liang et al. (2015).

(4) P θ | Z P ( Z | θ ) P ( θ )

(5) P Z | θ exp - 1 2 σ i 2 ( t ) i = 1 n t obs ( Z ) [ Z i t - X i t ] 2

Here Zi(t) and Xi(t) are the observed and modeled transit times, and σi2(t) is the standard deviation of measurements. The posterior probability density function of the parameters was constructed with two steps: a proposing step and a moving step. In the first step, the dataset was generated based on the previously accepted data with a proposal distribution:

(6) θ new = θ new + d ( θ max - θ min ) D ,

where θmax and θmin are the maximum and minimum values of the given parameters, d is the random variable between −0.5 and 0.5 with uniform distribution, and D is used to control the proposing step size in this study. In the moving step, the new data θnew is tested against the Metropolis criteria to quantify whether it should be accepted or rejected. The parameters of posterior probability density function were constructed with the Metropolis–Hastings algorithm. The Metropolis–Hastings algorithm was run 50 000 times for observed data. Accepted parameter values were used in further analysis.

Based on the concepts of mean age and mean transit time published by Rasmussen et al. (2016) and Lu et al. (2018), the mean carbon age defined as the whole time carbon atoms are stored in the carbon pools and then the mean age of carbon ait in a certain carbon pool i could be calculated with Eq. (7):

(7) a i t = 1 + i = 1 3 ( a j t - a i t ) f i j t C i - a j t I i ( t ) C i ,

where the fij(t) values are the carbon fraction transfer coefficients from jth to ith pools, and Ii(t) is the external input into the ith carbon pool. The transit time τi(t) was defined as the mean age of carbon atoms leaving the carbon pool at a specific time:

(8) τ i t = i = 1 d f i ( t ) a i t ,

where the fi(t) is the fraction of carbon with mean age ai(t).

2.4 Matrix approach through CLM4.5 and CLM4.5_noV

The Community Land Model version 4.5 (CLM4.5) is the terrestrial component of the Community Earth System Model (CESM). This version mainly consists of exchanges among different carbon and nitrogen pools and other biogeochemical cycles and includes a vertical dimension of soil carbon and nitrogen transformations (Koven et al., 2013). The matrix approach was applied to extract the soil module from the original CLM4.5, which could evaluate which processes influence τsoil in the model (Huang et al., 2017). Once we obtain the total carbon pool and Rh in each pool, we can calculate the τsoil with Eq. (1). We represented the structure of SOC as seven carbon pools as (i) one coarse woody debris (CWD) pool, (ii) three litter pools (litter1, litter2, and litter3), and (iii) three soil carbon pools (soil1, soil2, and soil3). In this matrix, carbon is transferred from three litter pools and CWD to three soil pools with different transfer rates. In each layer, these transfer rates are regulated by the transfer coefficients and fractions. C inputs from litterfall were allocated into different C compartments by modifications by soil environmental factors (temperature, moisture, nitrogen, and soil oxygen) and vertical transfer process. To understand whether the incorporation of soil vertical profile affects the simulation of τsoil, we compared the results based on the matrix approach with (i.e., CLM4.5) or without (i.e., CLM4.5_noV) the soil vertical transfer process.

In CLM4.5, soil C dynamics was simulated with 10 soil layers, and the same organic matter pools among different vertical soil layers are allowed to mix mainly through diffusion and advection. The matrix approach determines the soil dynamic of each SOC pool by simulating the first-order kinetics as Eq. (9):

(9) d C t d t = B t I ( t ) - A ξ t KC t V ( t ) C ( t ) ,

where the C(t) is the organic carbon pool size at time t. I(t) is the total organic carbon inputs while B(t) is the vector of partitioning coefficients. K is a diagonal matrix that represents the intrinsic decomposition rate of each carbon pool. The decomposition rate in the matrix approach is modified by the transfer matrix A and environmental scalars ξ. The scalar matrix ξ shown in Eq. (10) is the environmental factor to modify the SOC intrinsic decomposition rate. Each scalar matrix combines temperature (ξT), water (ξW), oxygen (ξO), depth (ξD), and nitrogen (ξN) controlled scalars for SOC decay.

(10) ξ = ξ T ξ W ξ O ξ D ξ N

A is the horizontal carbon transfer matrix, which quantifies C movement among different carbon pools shown as matrix (10). The non-diagonal entries Aij shown in matrix (10) represent the fraction of carbon that moves from the jth to the ith pools. In CLM4.5 and CLM4.5_noV, transfer coefficients are the same in each soil layer.

(11) A = 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 f 44 0 0 0 0 f 52 f 53 0 f 55 f 56 f 57 0 0 0 f 64 f 65 f 66 0 0 0 0 0 f 75 f 76 f 77

V(t) is the vertical carbon transfer coefficient matrix among different soil layers, and each of the diagonal blocks is a tridiagonal matrix that describes transfer coefficients with Vij(t). In this section, CLM4.5_noV assumes no vertical transfers in all pools. Therefore, V(t) for CLM4.5_noV is a blank matrix in the simulation. In contrast, CLM4.5 was assigned by a matrix with vertical transfers in each C pool. The vertical transfer rates among different C pool categories in CLM4.5 are matrix (12).


2.5 Statistical analyses

The median and interquartile were used for the quantification of both observational and modeling results due to the fact that the probability distribution of τsoil is not normal. To test the difference in τsoil among three approaches, we first normalized the data with the log-transformation and then applied the one-way ANOVA with a multi-comparison technique (Fig. 1b insert). The linear regression and correlation analyses were performed in R (3.2.1; R development Core team, 2015).

The Gaussian kernel density estimation was used to obtain the distributions of observed transit times (Sheather and Marron, 1990; Saoudi et al., 1997). The Gaussian kernel density estimation is a nonparametric approach to estimate the probability density function of a random variable. Let (x1x2, xn) denote the observed SOC τsoil with density function f as below:

(13) f ^ h x = 1 n h i = 1 n K ( x - x i h ) ,

where K is the nonnegative function than integrates to one and has a mean of zero, and h>0 is a smoothing parameter called the bandwidth. The bandwidth for approaches of stable isotope 13C, stock over flux, and incubation are 48.61, 35.13, and 2.62, respectively.

3 Results and discussion

3.1 τsoil and its spatial variation using different approaches

The one-way ANOVA with multi-comparison analysis showed no significant difference in the log-transformed τsoil between the methods of 13C (median = 60 years; interquartile range = 8 to 29 years) and stock over flux (16; 3 to 156 years, Fig. 1b). The range of these field in situ measurements (31; 5 to 84 years) is comparable to a former estimate of mean SOC turnover time (48 with 24 to 107 years) across 20 long-term experiments in temperate ecosystems using the 13C labeling approach (Schmidt et al., 2011). However, the estimates of τsoil from laboratory studies (4; 1 to 15 year) were significantly shorter than those from the other two methods (Fig. 1b). It suggests that the τsoil could be underestimated by the measurements from the laboratory incubation studies. Thus, the τsoil values from the laboratory incubation studies were excluded in the following analyses.

Figure 2Global spatial variation of SOC transit time (τsoil) with climate and the difference of τsoil estimation between observations and models. (a) Spatial variation of τsoil with mean annual temperature (MAT) and mean annual precipitation (MAP). (b) Comparisons of modeled against observed τsoil. Details for the classification of biomes are provided in the method section.


We then integrated the estimates of τsoil based on the 13C and stock-over-flux approaches to examine the inter-biome difference. As shown by Fig. 2b, the longest τsoil was found in desert and shrubland (170; 58 to 508) and tundra (159; 39 to 649 years). Boreal forest (58; 25 to 170 years) has a longer τsoil than the temperate (44; 13 to 89 years) and tropical forests (15; 9 to 130 years). Grassland and savanna had short τsoil (35; 21 to 57 years) and croplands had moderate (62; 21 to 120 years) τsoil in comparison with other biomes (Fig. 2).

3.2 Modeled τsoil in the CMIP5 ensemble and its estimation biases

The longest ensemble mean τsoil values of multiple models were found in dry and cold regions (Fig. 2). In comparison with the integrated observations from 13C and stock over flux, the modeled τsoil values were significantly shorter across all biomes (Fig. 2b insert). The negative bias was larger in dry (desert, grassland, and savanna) and cold (tundra and boreal forest) regions than tropical and temperate forests. The longest modeled τsoil appeared in the tundra ecosystem with a median of 64 years. The modeled median τsoil values were also shorter than observations in tropical forest (9 years), temperate forests (13 years), boreal forest (24 years), grassland/savanna (25 years), desert and shrubland (58 years), and croplands (27 years) (Fig. 2). In comparison with the observations, the models obviously underestimated the τsoil in the cold and dry biomes (Fig. 2b). A recent global data–model comparison study at 0.5×0.5 resolution also detected a similar spatial pattern of underestimation bias in ecosystem C turnover time (Carvalhais et al., 2014), but its magnitudes of bias in the cold regions are much smaller than those found in this study.

By grouping the τsoil into different climatic categories, we found that the observed τsoil significantly covaried with MAT (y=-5.28x+156.04, r2=0.48, P<0.01) and MAP (y=-68.19x+1222.6, r2=0.60, P<0.01) (Fig. 3). These results support the previous findings of negative covariations between τsoil and temperature at both the site and global levels (Trumbore et al., 1996). Although there is no significant correlation between τsoil and MAP in the observations, the models produced negative correlations of τsoil with MAT (r2=0.24, P<0.05) and MAP (r2=0.44, P<0.05) (Fig. 3).

Figure 3Relationships between SOC transit time (τsoil) and climate factors in both observations and CIMP5 models. The black solid lines show the negative correlation between τsoil and (a) mean annual temperature and (b) mean annual precipitation. The black dots indicate the aggregated τsoil over each category of MAT (y=-5.47x+1971.5, r2=0.49, P<0.01) or MAP (y=-68.19x+1222.6, r2=0.60, P<0.01). The red and blue dots present the mean value of multiple models based on the ratios of carbon stock over NPP and Rh, respectively.


3.3 Estimation of the τsoil with a three-pool model

With the three-pool model, the total C stocks and CO2 efflux from observations and the CMIP5 ensemble were separated into pool-specific decomposition rates by the deconvolution analysis (Fig. S3a; Liang et al., 2015). Seven out of 11 parameters were constrained for tropical forest and cropland (Fig. S4, Fig. S9). Eight out of 11 parameters were constrained for temperate, boreal forest, and desert and shrubland (Figs. S5, S6, S8). Five out of 11 parameters were constrained for the tundra ecosystem (Fig. S7). For grassland and savanna, seven out of 11 parameters were constrained (Fig. S10).

Figure 4The SOC transit time (τsoil) calculated from the one- and three-pool models under the steady-state assumption.


The longest simulated τsoil appeared in tundra (167 years) and desert (135 years) (Fig. 4, Table S3). Temperate forest (79 years) has a longer τsoil than the boreal (66 years) and tropical forests (29 years). Grassland and savanna had short τsoil (53.8 years) and croplands had moderate (77 years) τsoil in comparison with other biomes. The τsoil calculated from the one- and three-pool models did not show a large difference across all biomes. Also, estimates based on these two model structures showed the largest underestimation of τsoil in the tundra and desert (Fig. 4).

Figure 5Simulated SOC transit time (τsoil) by CLM4 (a, median global τsoil=20.56 years), CLM4.5 (b, median global τsoil=127.50 years), and CLM4.5_noV (c, median global τsoil=22.24 years). Panel (d) shows the latitudinal spatial distribution of the mean τsoil of different models in desert and tundra. The inserts in (a)(c) compare the τsoil between models and observations. The bottom and top of the box represent the first and third quartiles.


3.4 Improved modeling of τsoil with vertically resolved SOC dynamics

Given that many ESMs have further developed their representations of the soil biogeochemistry in recent years, we also examined whether the τsoil estimates have been improved by one of the CMIP5 models (i.e., CESM). It is encouraging that the biases of τsoil in dry and cold regions have been substantially reduced in the new land version of CESM (i.e., version 4.5 of the Community Land Model, CLM4.5). One major improvement in CLM4.5 is the vertically resolved SOC dynamics (Koven et al., 2013). The soil organic carbon is allowed to transfer through diffusion and advection up to 3.8 m within 10 layers. In each layer, the transfer rates are regulated by the environmental scalars (i.e., temperature, soil moisture, and available oxygen). The τsoil values simulated by CLM4.5 are longer than those in CLM4 (with median value 137 years and 21 years), especially in northern high-latitudinal regions. By turning off the vertical C movements with a matrix approach (i.e., there is no vertical C transfer; thus, the vertical matrix is a zero matrix in Eq. 12), we showed a similar pattern of underestimation on τsoil by CLM4.5 (i.e., CLM4.5_noV in Fig. 5). Huang et al. (2017) also reported the longer τsoil and high carbon storage capacity in northern high latitudes. These results suggest that the vertically resolved soil biogeochemistry is promising in improving the τsoil estimates by ESMs. However, it should be noted that the spatial variation of τsoil is still largely underestimated by the CLM4.5 (Fig. 5b insert).

Higher NPP values simulated by ESMs in the cold and dry regions have been reported by previous studies (Shao et al., 2013; Smith et al., 2016; Xia et al., 2017). The models produce high NPP in cold regions largely because they overestimate the efficiency of plants transferring assimilated C to growth (Xia et al., 2017). The CMIP5 models overestimate the precipitation and underestimate the dryland expansion 4-fold during 1996–2005 (Ji et al., 2015), which could lead to high NPP and fast SOC turnover rates. These results suggest that once the NPP simulation is improved without the correction of the τsoil underestimation, the models will produce smaller SOC stock in the cold and dry ecosystems.

This study shows that adding the vertical resolved biogeochemistry is a promising approach to correct the bias of τsoil in current models. However, other processes such as microbial dynamics, SOC stabilization, and nutrient cycles could affect the estimation of τsoil but are so far fully considered by the CMIP5 models (Luo et al., 2016). For example, adding soil microbial dynamics could increase τsoil in cold regions by lowering the transfer proportion of decomposed SOC to the atmosphere (Wieder et al., 2013). By contrast, the incorporation of nitrogen cycles might shorten τsoil by increasing plant C transfers to short-lived litter pools (e.g., O-CN and CABLE models) (Gerber et al., 2010) or reducing litter C transfers to the slow soil C pools (e.g., LM3V model) (Xia et al., 2013).

Large challenges still exist in using observations derived from different methods to constrain the modeled τsoil. Laboratory incubation studies report much shorter τsoil than other methods, mainly due to the optimized soil moisture and/or temperature during the soil incubation (Stewart et al., 2008; Feng et al., 2016). This suggests that the ESM models will largely underestimate τsoil if its turnover parameters are derived from laboratory incubation studies. It should be noted that the observations from the 13C and the stock-over-flux approaches in this study are derived for the bulk soil. However, SOC is commonly represented as multiple pools with different cycling rates in most of the CMIP5 models (Luo et al., 2016; Sierra et al., 2017, 2018; Metzler and Sierra, 2018). As synthesized by Sierra et al. (2017), the observations of τsoil are useful for a specific model once its pool structure is identified. This study also detects a difference in the estimated τsoil between the one- and three-pool models (Fig. 4). Thus, model databases, such as bgc-md (, last access: 2 February 2019), are a useful tool to improve the integration of observations and soil C models. An enhanced transparency of C-cycle model structure in ESMs is highly recommended, especially when they participate in the future model intercomparison projects such as CMIP6 (Jones et al., 2016).

4 Conclusions

This study detected large underestimation biases of τsoil in ESMs in cold and dry biomes, especially the tundra and desert. Improving the modeling of SOC dynamics in these regions is important because the cold ecosystems (e.g., the permafrost regions) are critical for global C feedback to future climate change (Schuur et al., 2015) and the dry regions strongly regulate the interannual variability of land CO2 sink (Poulter et al., 2014; Ahlström et al., 2015). The current generation of ESMs represents the soil C processes with a similar model formulation as first-order C transfers among multiple pools (Sierra and Markus, 2015; Luo et al., 2016; Metzler and Sierra, 2018). Thus, tremendous research efforts are still required to attribute the underestimation biases of τsoil in current ESMs to their sources, such as model structure, parameterization, and climate forcing. Reducing these biases would largely improve the accuracy of ESMs in the projection of the future terrestrial C cycle and its feedback to climate change. Recent modeling activities aiming to increase the soil heterogeneity, for example, soil vertical profile (Koven et al., 2013, 2017) and microbial dynamics (Allison et al., 2010; Wieder et al., 2013), are promising. Overall, this study shows the great spatial variation of τsoil in natural ecosystems, and we recommend more research efforts to improve its representation by ESMs in the future.

Data availability

The data are available from the corresponding author by request.


The supplement related to this article is available online at:

Author contributions

JX designed the study. JW collected and organized the data. LJ provided the CMIP5 and HWSD data. XX provided the laboratory incubation data. YH provided the CLM4.5 matrix module. JW and JX wrote the first draft, and all other authors contributed to revision and discussion of the results.

Competing interests

The authors declare that they have no conflict of interest.


We appreciate the anonymous reviewers for their valuable suggestions. We also appreciate Katherine Todd-Brown for her support of the soil data in CMIP5 and Deli Zhai for valuable comments. The model simulations analyzed in this study were obtained from the Earth System Grid Federation CMIP5 online portal hosted by the Program for Climate Model Diagnosis and Intercomparison at Lawrence Livermore National Laboratory (, last access: 4 October 2018). This work was financially supported by the National Natural Science Foundation (31722009, 31800400, 41630528), the National Key R&D Program of China (2017YFA0604603), the Fok Ying Tong Education Foundation for Young Teachers in the Higher Education Institutions of China (grant no. 161016), and the National 1000 Young Talents Program of China.

Edited by: Jens-Arne Subke
Reviewed by: two anonymous referees


Ahlström, A., Raupach, M. R., Schurgers, G., Smith, B., Arneth, A., Jung, M., Reichstein, M., Canadell, J. G., Friedlingstein, P., Jain, A. K., and Kato, E.: The dominant role of semi-arid ecosystems in the trend and variability of the land CO2 sink, Science, 348, 895–899,, 2015. 

Allison, S. D., Matthew, D. W., and Mark, A. B.: Soil-carbon response to warming dependent on microbial physiology, Nat. Geosci., 3, 336–340,, 2010. 

Balesdent, J., Mariotti, A., and Guillet, B.: Natural 13C abundance as a tracer for studies of soil organic matter dynamics, Soil Biol. Biochem., 19, 25–30,, 1987. 

Bernstein, L., Bosch, P., Canziani, O., Chen, Z., Christ, R., and Riahi, K.: IPCC, 2007: Climate Change 2007: Synthesis Report, 2008. 

Bloom, A. A., Exbrayat, J. F., van der Velde, I. R., Feng, L., and Williams, M.: The decadal state of the terrestrial carbon cycle: Global retrievals of terrestrial carbon allocation, pools, and residence times, P. Natl. Acad. Sci. USA, 113, 1285–1290,, 2016. 

Bolin, B. and Henning, R.: A note on the concepts of age distribution and transit time in natural reservoirs. Tellus, 25, 58–62,, 1973. 

Bolker, B. M., Pacala, S. W., and Parton Jr., W. J.: Linear analysis of soil decomposition: insights from the century model, Ecol. Appl., 8, 425–439, 1998. 

Bradford, M. A., Wieder, W. R., Bonan, G. B., Fierer, N., Raymond, P. A., and Crowther, T. W.: Managing uncertainty in soil carbon feedbacks to climate change, Nat. Clim. Change, 6, 751–758,, 2016. 

Carvalhais, N., Forkel, M., Khomik, M., Bellarby, J., Jung, M., Migliavacca, M., Mu, M., Saatchi, S., Santoro, M., Thurner, M., and Weber, U.: Global covariation of carbon turnover times with climate in terrestrial ecosystems, Nature, 514, 213–217,, 2014. 

Ciais, P., Sabine, C., Bala, G., Bopp, L., Brovkin, J., and Thornton, P.: Climate Change 2013: the physical science basis. Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change, edited by: Stocker, T. F., Qin, D., Plattner, G. K., Tignor, M., Allen, S. K, Boschung, J., Nauels, A., Xia, Y., Bex, V., and Midgley, P. M., Cambridge Univ. Press, 465–570, 2013. 

FAO/IIASA/ISRIC/ISSCAS/JRC, Harmonized World Soil Database (version 1.10), FAO, Rome, Italy and IIASA, Laxenburg, Austria, 2012. 

Feng, W., Shi, Z., Jiang, J., Xia, J., Liang, J., Zhou, J., and Luo, Y.: Methodological uncertainty in estimating carbon turnover times of soil fractions, Soil Biol. Biochem., 100, 118–124,, 2016. 

Friedl, M. A., Sulla-Menashe, D., Tan, B., Schneider, A., Ramankutty, N., Sibley, A., and Huang, X.: MODIS Collection 5 global land cover: Algorithm refinements and characterization of new datasets, Remote Sens. Environ., 114, 168–182,, 2010. 

Friedlingstein, P., Cox, P., Betts, R., Bopp, L., von Bloh, W., Brovkin, V., Cadule, P., Doney, S., Eby, M., Fung, I., and Bala, G.: Climate-carbon cycle feedback analysis: Results from the C4MIP model intercomparison, J. Clim., 19, 3337–3353,, 2006. 

Fröberg, M., Tipping, E., Stendahl, J., Clarke, N., and Bryant, C.: Mean residence time of O horizon carbon along a climatic gradient in Scandinavia estimated by 14C measurements of archived soils, Biogeochemistry, 104, 227–236,, 2011. 

Gerber, S., Hedin, L. O., Oppenheimer, M., Pacala, S. W., and Shevliakova, E.: Nitrogen cycling and feedbacks in a global dynamic land model, Global Biogeochem. Cy., 24, GB1001,, 2010. 

He, Y., Trumbore, S. E., Torn, M. S., Harden, J. W., Vaughn, L. J., Allison, S. D., and Randerson, J. T.: Radiocarbon constraints imply reduced carbon uptake by soils during the 21st century, Science, 353, 1419–1424,, 2016. 

Huang, Y., Lu, X., Shi, Z., Lawrence, D., Koven, C.D., Xia, J., Du, Z., Kluzek, E., and Luo, Y.: Matrix approach to land carbon cycle modeling: A case study with Community Land Model, Glob. Change Biol., 24, 1394–1404,, 2017. 

Hutchinson, M. F. and Xu, T.: Anusplin version 4.2 user guide. Centre for Resource and Environmental Studies, The Australian National University, Canberra, 54, 2004. 

Ji, M., Huang, J., Xie, Y., and Liu, J.: Comparison of dryland climate change in observations and CMIP5 simulations, Adv. Atmos. Sci., 32, 1565–1574,, 2015. 

Jones, C. D., Arora, V., Friedlingstein, P., Bopp, L., Brovkin, V., Dunne, J., Graven, H., Hoffman, F., Ilyina, T., John, J. G., Jung, M., Kawamiya, M., Koven, C., Pongratz, J., Raddatz, T., Randerson, J. T., and Zaehle, S.: C4MIP – The Coupled Climate-Carbon Cycle Model Intercomparison Project: experimental protocol for CMIP6, Geosci. Model Dev., 9, 2853–2880,, 2016. 

Koven, C. D., Riley, W. J., Subin, Z. M., Tang, J. Y., Torn, M. S., Collins, W. D., Bonan, G. B., Lawrence, D. M., and Swenson, S. C.: The effect of vertically resolved soil biogeochemistry and alternate soil C and N models on C dynamics of CLM4, Biogeosciences, 10, 7109–7131,, 2013. 

Koven, C. D., Hugelius, G., Lawrence, D. M., and Wieder, W. R.: Higher climatological temperature sensitivity of soil carbon in cold than warm climates, Nat. Clim. Change, 7, 817–822,, 2017. 

Liang, J., Li, D., Shi, Z., Tiedje, J. M., Zhou, J., Schuur, E. A. G., Konstantinidis, K. T., and Luo, Y.: Methods for estimating temperature sensitivity of soil organic matter based on incubation data: A comparative evaluation, Soil Biol. Biochem., 80, 127–135,, 2015. 

Lu, X., Wang, Y.-P., Luo, Y., and Jiang, L.: Ecosystem carbon transit versus turnover times in response to climate warming and rising atmospheric CO2 concentration, Biogeosciences, 15, 6559–6572,, 2018. 

Luo, Y., White, L. W., Canadell, J. G., DeLucia, E. H., Ellsworth, D. S., Finzi, A., Lichter, J., and Schlesinger, W. H.: Sustainability of terrestrial carbon sequestration: A case study in Duke Forest with inversion approach, Global Biogeochem. Cy., 17, 12101–2113,, 2003. 

Luo, Y., Ahlström, A., Allison, S. D., Batjes, N. H., Brovkin, V., Carvalhais, N., Chappell, A., Ciais, P., Davidson, E. A., Finzi, A., and Georgiou, K.: Toward more realistic projections of soil carbon dynamics by Earth system models, Glob. Biogeochem. Cy., 30, 40–56,, 2016. 

Metzler, H. and Sierra, C. A.: Linear autonomous compartmental models as continuous-time Markov chains: transit-time and age distributions, Math. Geosci., 50, 1–34, 2018. 

Metzler, H., Müller, M., and Sierra, C. A.: Transit-time and age distributions for nonlinear time-dependent compartmental systems, P. Natl. Acad. Sci. USA, 22, 201705296,, 2018. 

Mishra, U., Drewniak, B., Jastrow, J. D., Matamala, R. M., and Vitharana, U. W. A.: Spatial representation of organic carbon and active-layer thickness of high latitude soils in CMIP5 earth system models, Geoderma, 300, 55–63,, 2017. 

NASA LP DAAC Land Cover Type Yearly L3 Global 0.05 Deg CMG (MCD12C1), USGS/Earth Resources Observation and Science (EROS) Center, Sioux Falls, South Dakota, available at: products table/land cover/yearly l3 global 0.05 deg cmg/mcd12c1 (last access: 14 April 2014), 2008. 

Parry, M., Parry, M. L., Canziani, O., Palutikof, J., Van der Linden, P., and Hanson, C.: Climate Change 2007: Impacts, Adaptation and Vulnerability, Assessment Report of the Intergovernmental Panel on Climate Change, Cambridge Univ. Press, Cambridge, UK, 211–272, 2007. 

Poulter, B., Frank, D., Ciais, P., Myneni, R. B., Andela, N., Bi, J., Broquet, G., Canadell, J. G., Chevallier, F., Liu, Y. Y., and Running, S. W.: Contribution of semi-arid ecosystems to interannual variability of the global carbon cycle, Nature, 509, 600–603,, 2014. 

Rasmussen, M., Hastings, A., Smith, M. J., Agusto, F. B., Chen-Charpentier, B. M., Hoffman, F. M., Jiang, J., Todd-Brown, K. E., Wang, Y., Wang, Y. P., and Luo, Y.: Transit times and mean ages for nonautonomous and autonomous compartmental systems, J. Math. Biol., 73, 1379–1398,, 2016. 

Sanderman, J. Ronald, G. A., and Dennis, D. B.: Application of eddy covariance measurements to the temperature dependence of soil organic matter mean residence time, Glob. Biogeochem. Cy., 17, 301–3015,, 2003. 

Saoudi, S., Ghorbel, F., and Hillion, A.: Some statistical properties of the kernel – diffeomorphism estimator, Appl. Stoch. Model Data Anal., 13, 39–58,<39::AID-ASM292>3.0.CO;2-J, 1997. 

Schmidt, M. W., Torn, M. S., Abiven, S., Dittmar, T., Guggenberger, G., Janssens, I. A., Kleber, M., Kögel-Knabner, I., Lehmann, J., Manning, D. A., and Nannipieri, P.: Persistence of soil organic matter as an ecosystem property, Nature, 478, 49–56,, 2011. 

Schuur, E. A. G., McGuire, A. D., Schädel, C., Grosse, G., Harden, J. W., Hayes, D. J., Hugelius, G., Koven, C. D., Kuhry, P., Lawrence, D. M., and Natali, S. M.: Climate change and the permafrost carbon feedback, Nature, 520, 171–179,, 2015. 

Shao, P., Zeng, X., Sakaguchi, K., Monson, R. K., and Zeng, X.: Terrestrial carbon cycle: climate relations in eight CMIP5 earth system models, J. Clim., 26, 8744–8764,, 2013. 

Sheather, S. J. and Marron, J. S.: Kernel quantile estimators, J. Am. Stat. Assoc., 85, 410–416,, 1990. 

Sierra, C. A. and Markus, M.: A general mathematical framework for representing soil organic matter dynamics, Ecol. Monogr., 85, 505–524,, 2015. 

Sierra, C. A., Müller, M., Metzler, H., Manzoni, S., and Trumbore, S. E.: The muddle of ages, turnover, transit, and residence times in the carbon cycle, Glob. Change Biol., 23, 1763–1773,, 2017. 

Sierra, C. A., Ceballos-Núñez, V., Metzler, H., and Müler, M.: Representing and understanding the carbon cycle using the theory of compartmental dynamical systems, J. Adv. Model. Earth Sy., 10, 1729–1734,, 2018. 

Six, J. and Jastrow, J. D.: Organic matter turnover, Encycl. of Soil Science, 2002, 936–942, 2002. 

Smith, W. K., Reed, S. C., Cleveland, C. C., Ballantyne, A. P., Anderegg, W. R., Wieder, W. R., Liu, Y. Y., and Running, S. W.: Large divergence of satellite and Earth system model estimates of global terrestrial CO2 fertilization, Nat. Clim. Change, 6, 306–310,, 2016. 

Spohn, M. and Sierra, C. A.: How long do elements cycle in terrestrial ecosystems?, Biogeochemistry, 139, 69–83,, 2018.  

Stewart, C. E., Paustian, K., Conant, R. T., Plante, A. F., and Six, J.: Soil carbon saturation: evaluation and corroboration by long-term incubations, Soil Biol. Biochem., 40, 1741–1750,, 2008. 

Todd-Brown, K. E. O., Randerson, J. T., Post, W. M., Hoffman, F. M., Tarnocai, C., Schuur, E. A. G., and Allison, S. D.: Causes of variation in soil carbon simulations from CMIP5 Earth system models and comparison with observations, Biogeosciences, 10, 1717–1736,, 2013. 

Trumbore, S. E.: Comparison of carbon dynamics in tropical and temperate soils using radiocarbon measurements, Glob. Biogeochem. Cy., 7, 275–290,, 1993. 

Trumbore, S. E., Chadwick, O. A., and Amundson, R.: Rapid exchange between soil carbon and atmospheric carbon dioxide driven by temperature change, Science, 272, 393–396,, 1996. 

Wieder, W. R., Bonan, G. B., and Allison, S. D.: Global soil carbon projections are improved by modelling microbial processes, Nat. Clim. Change, 3, 909–912,, 2013. 

Xia, J., Luo, Y., Wang, Y. P., and Hararuk, O.: Traceable components of terrestrial carbon storage capacity in biogeochemical models, Glob. Change Biol., 19, 2104–2116,, 2013. 

Xia, J., McGuire, A. D., Lawrence, D., Burke, E., Chen, G., Chen, X., Delire, C., Koven, C., MacDougall, A., Peng, S., and Rinke, A.: Terrestrial ecosystem model performance in simulating productivity and its vulnerability to climate change in the northern permafrost region, J. Geophys. Res., 122, 430–446,, 2017. 

Xu, X., Shi, Z., Li, D., Rey, A., Ruan, H., Craine, J. M., Liang, J., Zhou, J., and Luo, Y.: Soil properties control decomposition of soil organic carbon: results from dataassimilation analysis, Geoderma, 262, 235–242,, 2016. 

Zhang, K., Dang, H., Zhang, Q., and Cheng, X.: Soil carbon dynamics following landuse change varied with temperature and precipitation gradients: evidence from stable isotopes, Glob. Change Biol., 21, 2762–2772,, 2015. 

Short summary
Soil is critical in mitigating climate change mainly because soil carbon turns over much slower in soils than vegetation and the atmosphere. However, Earth system models (ESMs) have large uncertainty in simulating carbon dynamics due to their biased estimation of soil carbon transit time (τsoil). Here, the τsoil estimates from 12 ESMs that participated in CMIP5 were evaluated by a database of measured τsoil. We detected a large spatial variation in measured τsoil across the globe.
Final-revised paper