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

Tiantong National Forest Ecosystem Observation and Research Station & Research Center for Global Change and Ecological Forecasting, School of Ecological and Environmental Sciences, East China Normal University, Shanghai 200062, China Laboratoire des Sciences du Climat et de l’Environnement, 91191 Gif-sur-Yvette, France Center for ecosystem science and society, Northern Arizona University, Arizona, Flagstaff, AZ 86011, 10 USA College of Biology and the Environment, Nanjing Forestry University, Nanjing 210037, China Environmental Sciences Division & Climate Change Science Institute, Oak Ridge National Laboratory, Oak Ridge, Tennessee 37830, USA CSIRO Ocean and Atmosphere, PMB #1, Aspendale, Victoria 3195, Australia 15 Wuhan Botanical Garden, Chinese Academy of Sciences, Wuhan 430074, Hubei Province, China 8 Department of Earth System Science, Tsinghua University, Beijing 100084, China


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 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 differ-Published by Copernicus Publications on behalf of the European Geosciences Union.
ence 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 .
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 (R h ). 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 ( 13 C) after successive changes in C 3 -C 4 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 ( 14 C) (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;. 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 . Thus, this study mainly collects the τ soil values from the approaches of stock over flux, 13 C 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 onepool 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.

A global database of site-level τ soil
We collected the literatures that reported the τ soil based on measurements (Text S1 in the Supplement): (1) δ 13 C shifts after successive changes in C 3 -C 4 vegetation, (2) measurements of CO 2 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 World-Clim database version 1.4 (http://worldclim.org/, 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 km 2 (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  (7) croplands. Other land cover types like permanent wetland, urban, and bare land were not included in this study.

Outputs of Earth system models from CMIP5
The historical simulation outputs of 12 ESMs participating in CMIP5 from 1850 to 1860 (https://esgf-data.dkrz.de/ search/cmip5-dkrz/, last access: 11 January 2016) were analyzed in this study (Table S1). For each model, the SOC, litter C, NPP, and heterotrophic respiration (R h ) 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 R h ): (1)

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 as-sumption, 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 CO 2 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 where the matrix C(t) = (C 1 (t), C 2 (t), C 3 (t)) T is used to describe soil carbon pool sizes. A is a matrix given by The elements f ij are carbon transfer coefficients, indicating the fractions of the carbon entering the ith (row) pool from the j th (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 (k 1 , k 2 , k 3 ). 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).
Here Z i (t) and X i (t) are the observed and modeled transit times, and σ 2 i (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: 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 a i (t) in a certain carbon pool i could be calculated with Eq. (7): where the f ij (t) values are the carbon fraction transfer coefficients from j th to ith pools, and I i (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: where the f i (t) is the fraction of carbon with mean age a i (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 R h 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): 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.
A is the horizontal carbon transfer matrix, which quantifies C movement among different carbon pools shown as matrix (10). The non-diagonal entries A ij shown in matrix (10) represent the fraction of carbon that moves from the j th to the ith pools. In CLM4.5 and CLM4.5_noV, transfer coefficients are the same in each soil layer.
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 V ij (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).

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 (x 1 x 2 · · ·, x n ) denote the observed SOC τ soil with density function f as below: 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 13 C, stock over flux, and incubation are 48.61, 35.13, and 2.62, respectively.

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 13 C (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 13 C 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. We then integrated the estimates of τ soil based on the 13 C 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).

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 13 C 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.

Estimation of the τ soil with a three-pool model
With the three-pool model, the total C stocks and CO 2 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 crop-  = 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 R h , respectively. land (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).
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).

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 move- ments 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 . The CMIP5 models overestimate the precipitation and underestimate the dryland expansion 4fold 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 . 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) .
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 13 C 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 Sierra et al., 2017Sierra et al., , 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 (https: //github.com/MPIBGC-TEE/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 Ccycle model structure in ESMs is highly recommended, especially when they participate in the future model intercomparison projects such as CMIP6 (Jones et al., 2016).

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  and the dry regions strongly regulate the interannual variability of land CO 2 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;. 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 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.
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.