Quantifying uncertainties in soil carbon responses to changes in global mean temperature and precipitation

. Soil organic carbon (SOC)

the maintenance of SOC is important for global and social sustainability (e.g., Mol and Keesstra, 2012). In climate systems, because of the vast carbon pool of SOC, the behavior of SOC is key for understanding the feedback of terrestrial ecosystems to atmospheric CO 2 concentrations in a warmer world (Heimann and Reichstein, 2008;Thum et al., 2011). However, a large number of uncertainties exist in the observation and modeling of SOC dynamics (e.g., Post et al., 1982;Todd-Brown et al., 2013). For example, in the Coupled Model Intercomparison Project Phase 5 (CMIP5), Todd-Brown et al. (2013) reported that the (simulated) presentday global SOC stocks range from 514 to 3046 Pg C among 11 Earth system models (ESMs). Soil processes in terrestrial ecosystem models are significantly simpler than actual processes or above-ground processes, and thus exist structural uncertainties in SOC dynamics in ESMs.
Temperature and precipitation are critical factors for the feedback of terrestrial ecosystems to atmospheric CO 2 (Seneviratne et al., 2006). Similarly, SOC dynamics are strongly affected by temperature and precipitation, because SOC dynamics in biome models are parameterized as a function of soil temperature, moisture, and other factors (e.g., Davidson and Janssens, 2006;Ise and Moorcroft, 2006;Falloon et al., 2011). The differences in these functions and their parameters have important effects on the projection of global SOC stocks and their behavior (Davidson and Janssens, 2006;Ise and Moorcroft, 2006).
In this study, we examined the SOC dynamics simulated by seven biome models as part of the Inter-Sectoral Impact Model Intercomparison Project (ISI-MIP) (Warszawski et al., 2014), which were forced using the bias-corrected outputs of five global climate models (GCMs) in newly developed climate scenarios, i.e., representative concentration pathways (RCPs). We aimed to investigate the impact of climate change on the global SOC stock with respect to changes in global mean temperature and precipitation and explore the uncertainties in future global SOC stock projections.
In order to analyze the first-order behavior of the simulated global SOC-dynamics, we focused on the interannual responses of the biome models under the assumption that SOC is one-compartment of Earth's system. First, we considered global SOC dynamics as the following simple, differential equation: where Input is carbon derived primarily from photosynthesis products via chemical and microbial humification (Wershaw, 1993), and k is the global SOC turnover rate. In most conventional models (Li et al., 2014), SOC decomposition functions as a first-order decay process as in Eq. (1). SOC dynamics are regulated by the balance between the input from vegetation biomass carbon and SOC decomposition. In this study, we examined a simple hypothesis: can global mean temperature and precipitation anomalies ( T ( • C) and P (%), respectively) be used as explanatory variables of global SOC decomposition dynamics in future (projections over the 21st century). If true, this would mean that T and P can explain k during a projection period in biome models. This simplification enables us to review the global impact of climate change on SOC dynamics and identify the characteristics of biome models especially in global SOC behavior. Subsequently, we assessed whether the time evolution of the estimation uncertainties for SOC can be explained by T and P sensitivities during the 21st century for each biome model. Furthermore, we compared the spatial distributions of global SOC pools and their changes to evaluate regional differences, focusing on detailed processes in the interaction with vegetation dynamics.

Method and models
In this study, we examined SOC processes using seven biome models obtained from the ISI-MIP. The biome models are Hybrid4 (Friend and White, 2000), JeDi (Jena Diversity-Dynamic Global Vegetation model) (Pavlick et al., 2013), JULES (Joint UK Land Environment Simulator; Clark et al., 2011;, LPJmL (Lund-Potsdam-Jena managed land Sitch et al., 2003), SDGVM (Sheffield Dynamic Global Vegetation Model; Woodward et al., 1995), VISIT (Vegetation Integrative Simulator for Trace gases) (Ito and Oikawa, 2002;Ito and Inatomi, 2012), and ORCHIDEE (Organizing Carbon and Hydrology in Dynamic Ecosystems; Krinner et al., 2005). In this study, Hybrid4, JeDi, JULES, and LPJmL are dynamic global vegetation models, and the others are fixed vegetation models, in this study. General information about SOC processes is summarized in Table 1.
In the ISI-MIP framework, these models were run with 5 GCM × 4 RCP scenarios and a fixed CO 2 control was also run with RCP8.5 climate condition scenarios. In this study, for the biome model forcing, we used climate variables in HadGEM2-ES (HadGEM -Hadley Centre Global Environmental Model) with bias correction for temperature and precipitation from Hempel et al. (2013). For the spin-up of each model, we used de-trending forcing data for the years 1951-1980 repeatedly until reaching equilibrium of VegC (vegetation carbon) and SOC. For CO 2 , we used the CO 2 concentration for 1950 while running the 30 yr spin-ups. The global climate variables (atmospheric CO 2 concentration, global mean terrestrial temperature anomaly T ( • C), and global terrestrial precipitation anomaly P (%)) in each RCP scenario for HadGEM are summarized in Fig. 1. T and P were set to 0 as the averages of their values between 1980 and 2000. In addition, there was no anthropogenic land-use change for the entire simulation period in this study. More detail about the experimental setup is available in the literature (Warszawski et al., 2014).

Estimation of T and P sensitivity of global SOC
We used a state-space model (more properly vector autoregression) (Sims and Zha, 1998) to evaluate the sensitivity of global SOC decomposition to global temperature and precipitation anomalies in each biome model. This vector autoregression model considers only process uncertainty, not observation uncertainty in a state-space model. We applied this analysis to annual global SOC time-series data in each biome model simulated in the five scenarios (three scenarios for ORCHIDEE), i.e., the four RCPs and the fixed CO 2 experiment with RCP8.5 climate conditions in HadGEM (Figs. 1, 2).
We first modeled the likelihood function using the following equation. The model outputs were archived for each year; therefore, we discretized the equation as the annual time step t.
where SOC [n,t] is the global SOC stock at time t (year) in scenario n, and σ ps is the process error. µ [n,t−1] is defined as follows: where VegC [n,t] indicates the global vegetation biomass C stock at time t in scenario n, and α is the fraction of VegC transformed into SOC per year, which is assumed to represent the annual input of SOC. k is the turnover rate for global SOC (yr −1 ) under standardized global mean temperature and precipitation conditions (averages between 1980 and 2000). β 1 and β 2 are the global SOC sensitivities to T and P , respectively (units: yr −1 T −1 and yr −1 P −1 ). The priors of these parameters are defined as follows: k ∼ uniform (0, 1), β 2 ∼ normal (0, 100).  We used vague priors for β 1 and β 2 to estimate the T and P effect on k. For α and k, we used uniform priors, which are sufficiently broad theoretically.
Then, the joint posterior is given by following equation.

Evaluation of stimulated global SOC decomposition in T from posteriors
Using posteriors in the steady-state model, we simulated the global SOC decomposition stimulated by increased global mean temperature at 2, 3, and 4 • C.

Global SOC and VegC projection in HadGEM
The increase of T depends on the RCP scenario, with the maximum increase in RCP8.5 being 7.5 • C in 2099 in HadGEM2. In RCP2.6, the maximum T was 1.9 • C during the entire simulation period and showed signs of leveling off in 2050. In all RCP scenarios, P increased to 11 (RCP4.5) and 16 % (RCP8.5). However, there were high amplitudes of P within each RCP scenario; thus, there were no obvious differences between RCPs.
For 2000, in HadGEM, the global SOC stocks varied from 1090 (Hybrid4) to 2646 Pg C (JULES) between the biome models ( Fig. 2). The mean global SOC stock in the six models was 1772 Pg C (standard deviation; 568 Pg C). An estimated empirical global SOC stock was 1255 Pg C (Todd-Brown et al., 2013). However, global VegC stocks in 2000 ranged from 510 (VISIT) to 1023 Pg C (JULES). The mean global VegC among the seven biome models was 809 Pg C (SD (standard deviation); 223 Pg C) (Fig. 2). The global VegC stocks in most models were comparable with the VegC (493 Pg C) estimated by the IPCC Tier-1 method (Ruesch and Gibbs, 2008).
In the projection period , the SOC stock in the six models (except for Hybrid4) increased in all RCPs compared to that in 2000. The global SOC stock in Hybrid4 continuously decreased in all RCPs during the projection period (Fig. 2). Under the RCPs, the maximum SOC stock increase for the projection period was observed in JeDi with RCP8.5, with a value of 347 Pg C. In the fixed CO 2 scenarios, the global SOC stocks continuously decreased in most biome models, showing global SOC changes from −299 to 65 Pg C at the end of the simulation period.
The global VegC stocks increased in nearly all RCPs and biome models compared to the global VegC in 2000. However, the global VegC stocks in Hybrid4 and LPJmL with RCP8.5 did not continuously increase in the projection period and were not the largest stock at the end of the simulation period during the projection period. In the fixed CO 2 scenarios, the global VegC stocks also continuously decreased, and global VegC changes ranged from −517 to −40 Pg C at the end of the simulation period (Fig. 2).
The rank order of the SOC stock over each RCP at the end of the simulation (2099) is in good agreement with the rank order of each corresponding VegC stock in the same period in JeDi, JULES, LPJmL, and SDGVM. However, the orders of the SOC stock in the other biome models are different than those of the global VegC stocks. These stock changes are attributed to the different SOC decomposition processes.

Posteriors of the state-space model; global SOC sensitivity to T and P
The Gelman and Rubin convergence statistics (R) of all parameters were lower than 1.01 in all models; therefore, the parameters represented successful convergences (data not shown). The posterior distributions of the parameters for each biome model are summarized in Table 2. α, which is the fraction of annual translation of VegC to SOC, among the biome models varied from 0.721 % in Hy-brid4 to 3.860 % in VISIT. The SOC turnover rate k (yr −1 ) ranged from 2.51 × 10 3 in LPJmL to 16.10 × 10 3 yr −1 in VISIT.
The 95 % credible intervals (CI) in sensitivity of global SOC to T (β 1 ) in each biome model did not cover 0 in all models (Table 2). And the 95 % CI of β 1 in each model was not partially duplicated, which means that the sensitivity to T could be statistically distinguished between the biome models. The highest β 1 was observed in VISIT, with a median value of 1.225 × 10 −3 yr −1 T −1 (or • C −1 ). The lowest β 1 was observed in JeDi and was approximately 0 yr −1 T −1 .
The sensitivity of global SOC to P (β 2 ) in the biome models was lower compared to the SOC turnover rate k and β 1 . Their values (yr −1 P −1 ) were nearly one order of magnitude less than β 1 . Considering the range of the values of P in the projection period, the impact on global SOC stock dynamics is small in all biome models. Furthermore, the 95 % CIs of β 2 in each model were partially duplicated.  On the basis of the posterior parameters, we estimated the stimulated global SOC decomposition for 2, 3, and 4 • C, assuming that each global SOC stock is at the 2000 level (Fig. 3). A statistical difference was observed among the 2, 3, and 4 • C in five biome models (i.e., Hy-brid4, JULES, LPJmL, VISIT, and ORCHIDEE). However, the magnitudes of the stimulated global SOC decomposition varied. At 4 • C, it ranged from 1.9 (in LPJmL) to 8.1 Pg C yr −1 (in JULES). In SDGVM, there were no statistical differences in the stimulated global SOC decomposition between 3 and 4 • C. There were also no differences in this term among 2, 3, and 4 • C in JeDi.

Latitudinal δSOC (2099-2000 and CO 2 -fixed CO 2 )
in HadGEM RCP8.5 Latitudinal SOC stock in the HWSD (Figs. 4, 5a) displays a double peak in both the northern high latitudes and low latitudes. The most SOC stock is found around 60 • N. In all biome models, large SOC stocks were also observed in high-latitude zones (50-75 • N; Figs. 5a and S1 in the Supplement). However, the range of simulated SOC change during this century (kg C m −2 ) in each biome model was different. The upper 99 percentile of SOC accumulation in each biome model varied from 23.8 in SDGVM to 97.6 kg C m −2 in LPJmL ( Fig. S1 in the Supplement). For differences between 2099 and 2000 in HadGEM RCP8.5, a large variance among biome models was observed between 30 • S and 10 • N (tropic region) and between 40 and 75 • N (boreal to Arctic region) (Figs. 4,5a) in the biome models. There were four types of latitudinal changes: (i) SOC increase in both regions (JeDi, SDGVM, ORCHIDEE), (ii) SOC increase in boreal to arctic regions and decrease in the tropics (JULES, VISIT), (iii) SOC increase in the tropics and decrease in boreal to arctic regions (LPJmL), and (iv) SOC decrease in both regions (Hybrid4). The maximum difference was observed in the boreal regions, where it reached more than 20 Pg 2.5 • −1 .
There were also differences between the increasing CO 2 scenario (RCP8.5) and the fixed CO 2 scenario with the RCP8.5 climate condition in SOC ( SOC CO 2 −fixed CO 2 ) (Fig. 5c). This suggests that the increases of plant production and biomass due to CO 2 fertilizer effects in the increasing CO 2 scenario (RCP8.5) contributed to the SOC stock increases because of the increase of C input to soil (indirect CO 2 effect). We observed bimodal increases in six biome models, and the peaks were between 30 and 70 • N and between 30 • S and 10 • N. In Hybrid4, the large SOC increase due to CO 2 was unimodal around the boreal regions. The maximum difference between the increasing CO 2 scenario and the fixed CO 2 scenario was observed around 60 • N, which was approximately 10 Pg 2.5 • −1 The different values of SOC CO 2 −fixed CO 2 / VegC CO 2 −fixed CO 2 (Fig. 5d) indicate a different turnover rate of vegetation carbon to SOC (via litter) among the biome models and regions. This is because of the assumption of almost the same states except in VegC dynamics between RCP8.5 and fixed CO 2 scenarios. SOC CO 2 −fixed CO 2 / VegC CO 2 −fixed CO 2 varied with latitude and among the biome models. In almost all the models, SOC CO 2 −fixed CO 2 / VegC CO 2 −fixed CO 2 was the highest in the higher latitude regions. In the Hybrid4 model, the SOC CO 2 −fixed CO 2 / VegC CO 2 −fixed CO 2 was relatively low in all regions, compared with other model results.

Global mean temperature and precipitation impact(s) on global SOC decomposition and projection uncertainties
During the projection period , the global SOC stock changes in all RCPs (without the fixed CO 2 scenario) ranged from −6 to 280 Pg C under RCP2.6 (mean ± SD: 89 ± 104 Pg C). Under RCP8.5, the SOC changes varied from −124 to 392 Pg C (113 ± 176 Pg C) (Fig. 2) at the end of the projection period. These global SOC stock changes are equivalent to −185 to +58 ppmv in atmospheric CO 2 concentration. Thus, in higher radiative forcing scenarios, uncertainties associated with future global SOC projection increase. These ranges of the global SOC stock changes by 2099 were comparable with the VegC changes (Fig. 2). However, in the projection period, the global VegC stocks primarily act as sinks for atmospheric CO 2 , while the global SOC   stocks act as either sinks or sources depending on the biome model. There were similar SOC projections in the same period (2000-2100) from multiple model simulations in previous studies. In the C4MIP study, for example, the global SOC stock changes ranged from approximately −50 to 300 Pg by the end of the simulation period among the 11 coupled climate-carbon models (Friedlingstein et al., 2006;Eglin et al., 2010). It has also been predicted that SOC stocks in 2100 differ by approximately 200 Pg among five DGVMs under forced A1FI and B1 scenarios (Sitch et al., 2008), which is the highest forcing scenario in the AR4 assessment. Compared with these studies, the SOC changes simulated in this study varied comparably or showed slightly higher uncertainty than those of previous projections. The magnitude of global SOC decomposition and the response to T primarily depend on the amount of the  global SOC stock and a turnover rate of SOC decomposition process. As has been reported in a CMIP5 experiment (Todd-Brown et al., 2013), our study has also shown that simulated global present-day SOC stocks in seven ecosystem models show high variation (1090-2646 Pg C) compared to the variation of global present-day VegC stocks (Fig. 2). There were some estimations available for global SOC stock, ranging from 700 (Bolin, 1970) to 3000 Pg C (Bohn, 1976). The most widely cited studies (Post et al., 1982;Batjes, 1996)      that the sensitivity of global SOC to T varied among the biome models and that the present-day global SOC stock can be used to make more reliable SOC projections. Although actual global SOC stock estimation still has significant uncertainty, global SOC stock constraints are essential for reducing uncertainty in global SOC projections in ecosystem models.
Our simplified global dynamic model for the global SOC stock revealed that the balance of the global SOC stock turnover and input from VegC is quite different among the biome models, which further implies the different sensitivities to T of the global SOC stocks among the biome models (Table 1). Hybrid4-simulated global SOC stocks decrease by 2099 in all RCPs because of the relatively high T sensitivity in addition to the low turnover rate (high residence time) in VegC to SOC (Table 2, Friend et al., 2014). Although temperature is the most significant regulation factor of SOC dynamics (Raich and Schlesinger, 1992), discussion of the effect of increasing global mean temperature on SOC stocks is still lacking. According to our statistical analysis (Table 2), most biome models had adequate resolution to describe the global SOC stock change among the 1 • C (or 2 • C for SDGVM) difference in the projection period. In these models, the global mean temperature T could be a measure of the robustness of global SOC stock projection. However, the global SOC in JeDi was not sensitive to T in this projection period. According to our estimation, the highest global SOC sensitivity was observed in VISIT, in which the rate of global SOC stock change was enhanced by −6.95 Pg C yr −1 in 4 • C (Fig. 3). However, the highest magnitude of SOC decomposition stimulated by increasing T was observed in JULES (−8.13 Pg C yr −1 in 4 • C) due to high global SOC stock in JULES. The Carnegie-Ames-Stanford approach model showed global SOC decomposition sensitivity of 2.26 Pg C yr −1 • C −1 , which is nearly equivalent to results obtained from JULES when the 4 • C value was derived from simple extrapolation (Zhou et al., 2009). There is still a lack of observation-based estimation of global SOC response intensity to T . Both global SOC stocks and data-oriented parameters such in Raich et al. (2002) could represent important information for the constraint and validation of global SOC dynamics.
However, β 2 was not effective for global SOC dynamics in all ecosystem models in our analysis, which does not mean Earth Syst. Dynam., 5, 197-209 that precipitation is not important in SOC dynamics. Precipitation trends are globally heterogeneous; therefore, the representative P might not be a useful index of SOC stock dynamics at a global scale in this projection period. However, precipitation is quite important in both soil decomposition (Falloon et al., 2011) and vegetation processes (Seneviratne et al., 2006), which considerably contribute to regional SOC dynamics.

SOC stock changes from vegetation dynamics and regional aspect
There were consistent latitudinal (geographic) patterns among the biome models (Figs. 5a, 1), and the highest SOC stock was observed between 40 and 75 • N. However, we found that the amount of SOC stocks among the biome models significantly vary in this region. The models' SOC densities are different, possibly because of the balance of input and decomposition and the consideration of depth in the biome models (1 to 3 m or not explicit, Table 1). Tarnocai et al. (2009) estimated SOC stock depth up to 3 m, with a value of 1672 Pg C in permafrost-affected regions only. Thus, the SOC stock of this region and the global SOC stock in the biome models may be significantly underestimated. From a regional perspective, the biome models showed quite different spatial patterns of SOC changes under HadGEM RCP8.5 (Figs. 4,5), while the spatial patterns of VegC changes were generally more consistent among the biome models (Friend et al., 2014). We found that this spatial heterogeneity among the biome models was also present in the SOC stock changes in different scenarios (data not shown). In particular, in boreal to arctic regions, SOC acts as a sink and source of C depending on the biome model (Fig. 5). This result indicates that there is an underlying mechanistic difference among the biome models in these regions. Two models show decreased SOC stocks by 2099 in this region in HadGEM RCP8.5. LPJmL shows unique features in SOC stocks and changes in this region. This implies that high SOC accumulations (over 80 kg-C m −2 ) (Figs. 5 and S1 in the Supplement) will be reduced with decreasing VegC by 2099 (Fig. S2 in the Supplement) in this region. This trend would result in low water availability in the permafrost regions, because the prediction is based on a mechanistic permafrost scheme (Beer et al., 2007;Schaphoff et al., 2013). Because LPJmL incorporated a freeze-and-thaw thermodynamics explicitly in discrete layers, it can simulate vertical water and carbon distributions in the model. This scheme enables LPJmL to describe the surface soil water deficit due to permafrost melting. Whereas, in Hybrid4, SOC decomposition is the main factor contributing to reduced SOC in this region. Dynamic vegetation and freeze-thaw schemes are important for SOC dynamics in permafrost zones, because they provide more accurate prediction of the balance of C input from successive vegetation and old soil carbon decomposition (Schuur et al., 2008Schaphoff et al., 2013). However, in this study, dynamic vegetation and freeze-thaw schemes are only implemented in LPJmL. The potential release from SOC in permafrost regions could have a large impact on the global C cycle (Koven et al., 2011;Burke et al., 2012;MacDougall et al., 2012), and further model development is essential for the modification of projections for this region.
Previous extensive field research has shown that the CO 2 fertilizer effect on plant growth in higher CO 2 concentrations could also result in the accumulation of SOC (De Graaff et al., 2006). For the RCP8.5 climate forcing, the fixed CO 2 experiment suggested that the CO 2 fertilizer effect on plant production contributed considerably to the global SOC stock increase in all biome models. The indirect CO 2 fertilizer effect on the global SOC stock varied from 93 (Hybrid4) to 264 Pg C (VISIT) (mean ± SD; 196 ± 60 Pg C) at the end of the simulation period, while VegC stock increased from 295 to 645 Pg C (275 ± 150 Pg C) by 2099 because of increasing CO 2 (Figs. 2 and S2 in the Supplement). Thus, the CO 2 fertilizer effect on global SOC accumulation strongly affects the biome models, and further quantitative assessment might be needed. For example, Friend et al. (2014) focused their attention on the effects of CO 2 fertilizers on biomass production and turnover rate of biomass. In addition to the indirect CO 2 effects, other nutrient limitations (e.g., nitrogen and phosphorus) and their sensitivities could be large sources of uncertainty in SOC projection via vegetation production (Goll et al., 2012;Exbrayat et al., 2013). In our study framework, we cannot adequately validate these issues since only a few models consider them in their current versions (e.g., Hy-brid4). Therefore, further interactions must be validated to more comprehensively understand the uncertainty sources in SOC projection.
A large variance in SOC CO 2 −fixed CO 2 / VegC CO 2 −fixed CO 2 was observed among the biome models (Fig. 5d), suggesting that the vegetation-soil interactions including the vegetation turnover rate (Friend et al., 2014) and litter decomposition rate also had large uncertainties. This variance might cause an SOC projection difference among the biome models. To reduce these uncertainties, a more observation-based validation is desirable. For the litter decomposition process, for example, global database of the long-term intersite decomposition experiment team (LIDET) is one useful validation case study (Bonan et al., 2013). In addition, the process in SOC formation from alteration of litter via decomposition process (i.e., humification) and in their stabilization have not yet been implemented robustly in biome models when compared with actual SOC formation processes (Sollins et al., 1996;Six et al., 2002). This is another major process missing from the vegetation-soil interaction in biome models. To comprehensively address the biome model uncertainties in each successive process, the traceability framework developed by Xia et al. (2013) could be helpful.

SOC modeling issues
The accurate estimation of the present-day global SOC stock remains difficult because of a lack of appropriate broad and non-destructive investigation techniques to measure SOC stock, such as satellite-based remote sensing. In fact, current SOC was formed in slow turnover fractions over thousands of years (Trumbore, 2000). Therefore, when getting an initial SOC by the spin-up phase in biome models, there may not be enough information on the historical climate conditions and vegetation dynamics to duplicate in the entire SOC formation history. This is potentially one of the biggest issues for accurate estimation of SOC stock in biome models. In addition, observations of global long-term SOC stock dynamics for model validation are limited. Thus, it is very difficult to assess projected global SOC trends in each biome model. Therefore, in addition to quantitatively understanding the SOC stock, deductive inferences based on the extensive understanding of the processes are essential for minimizing uncertainties in SOC stock prediction. For example, the apparent variability in global SOC sensitivity to T may result from differences in model structures and parameters. Regarding temperature sensitivity and the magnitude of response to rising temperatures, the following topics require improvement: (i) SOC compartments and their turnover rates (Jones et al., 2005;Conant et al., 2011), (ii) the temperature sensitivity parameter (e.g., Q 10 ) (Davidson and Janssens, 2006;Allison et al., 2010), and (iii) soil temperature prediction (radiation, heat production by microbes) (Luke and Cox, 2011;Khvorostyanov et al., 2008). In addition, microbial dynamics are a key component for the temperature acclimation of SOC decomposition (Todd-Brown et al., 2012;Wang et al., 2013). The acclimation response of SOC decomposition by microbial physiology is not included in the biome models used in this study. For SOC accumulation, soil mineralogical properties control soil C turnover (Torn et al., 1997). However, the biome models do not exploit global soil classification information (i.e., volcanic or non-volcanic soils), which still has significant uncertainties (Guillod et al., 2012;Hiederer and Köchy, 2011). In this study, peat and wetland soils are not explicitly simulated because of the large simulation grid size. Because of large carbon stock and water regime changes in future climates in such ecosystems, the SOC and soilwater-holding capacity feedback should also be considered in the SOC process in biome models (Ise et al., 2008). The interactions between SOC decomposition and nutrients (nitrogen) are also influential factors for global SOC projection (Manzoni and Porporato, 2007).
However, the details of these processes are beyond the scope of this study; therefore, we did not explore these issues in depth. A more specific model intercomparison, such as an environmental-response-function-based assessment (e.g., Falloon et al., 2011;Sierra et al., 2012;Exbrayat et al., 2013) is recommended. Furthermore, land-use change is not included in our projection; however, the effect of land-use changes on SOC dynamics is critical (Eglin et al., 2010). Estimating land-use change with high confidence is essential for accurate global SOC stock projections and could be used as a basis for policies that moderate the impacts of climate change.

Conclusions
The uncertainties associated with SOC projections are significantly high. The projected global SOC stocks by 2099 act as CO 2 sources or sinks depending on the biome model, even though models have similarly simulated historical SOC trends. The uncertainties of the SOC changes increase with higher forcing scenarios, and the global SOC stock change varies from −157 to 225 Pg C in HadGEM under RCP8.5 across biome models.
By adopting the simplified approach of global SOC as one compartment in the Earth system we can understand the comprehensive characteristics of each biome model on a global scale. The magnitude of SOC responses to global mean temperature increase considerably differed depending on the biome model. Our results confirmed that the SOC process implementations are dissimilar among the biome models at the global scale. In addition, global precipitation anomalies could not explain the simulated future global SOC stock changes. Moreover, the indirect CO 2 fertilizer effect contributed strongly to global SOC stock changes and projection uncertainties. For more reliable projections, both SOC dynamics and vegetation processes require reliable global SOC stock estimation and region-based improvements.