Interactive comment on “ Convergent modeling of past soil organic carbon stocks but divergent projections ” by Z .

We appreciate the positive comments and the suggestions on microbial biomass from H. Ibrahim. The response to each of the three points described by H. Ibrahim follows: 1.In this study, we did not directly compare the ASPIM model to other models. As mentioned in Page4250 Line5, we stated that the APSIM is “. . .similar to other SOM models like RothC and Century. . .”. The following sentence described that each pool is “. . .treated as first-order decay process. . .”. As the first-order decay process is commonly used in soil carbon models, and is generally well-known and accepted by carbon modelling community, we did not go to the details about the equation. 2.Thanks for this question. Actually, incorporating microbial processes into soil carbon


Introduction
Soil is the largest carbon (C) reservoir in the terrestrial biosphere, and CO 2 emission from soil organic matter (SOM) decomposition accounts for ∼ 35 % of the global CO 2 emissions (Schlesinger and Andrews, 2000).Due to the large amount of soil organic carbon (SOC), carbon sequestration in soils represents a great potential for mitigating greenhouse gas emissions and climate change as well as maintaining soil fertility (Lal, 2004).Accurate projections of future change in SOC are therefore needed for C and greenhouse gas (GHG) inventories to guide the development of future policies and land management practices (Janssens et al., 2003).Due to the complex and dynamic interactions between SOC, climate, soil and land management practices, processbased SOM models have become an important tool to investigate SOC change and project SOC trends under different land uses (Jenkinson et al., 1991;Friedlingstein et al., 2006;Smith et al., 2007;Piao et al., 2009).Some studies have suggested that the uncertainties in such projections should be systematically addressed in order to judge the credibility of the underlying projections and develop appropriate polices for carbon sequestration and climate change mitigation (Friedlingstein et al., 2006;Tang and Zhuang, 2008;Todd-Brown et al., 2013;Nishina et al., 2014).A better understanding of these Z. Luo et al.: Uncertainty in soil C projections uncertainties and their drivers will help identify knowledge gaps and improve process-based models (Luo et al., 2015).
Uncertainty in simulation results derived from dynamic models can arise from inaccuracies in input data, deficiencies in model structure and inappropriate optimization of model parameters.For SOM models, initialization of the SOM pools can also be a major cause of divergent model projections.Most SOM models divide SOM into several conceptual pools (e.g.fast, slow and recalcitrant pools) and simulate the decomposition of each pool as a first-order decay process (Smith et al., 1997;Davidson and Janssens, 2006;Schmidt et al., 2011).In many cases, measurements are only available for total SOC, and there is no agreed-on procedure for initialization of these model pools using total SOC (Basso et al., 2011).As a result, model optimization was often conducted based on limited SOC measurements (usually at temporal scales less than decades) together with empirical initialization.The optimized model was then used to project SOC change at wider spatiotemporal scales (Friedlingstein et al., 2006;Thornton et al., 2007).Such projection is subject to unknown uncertainty (Friedlingstein et al., 2006;Tang et al., 2008;Luo et al., 2013), because it does not properly address the inaccuracies in both model initialization and model parameters, with the latter potentially caused by imperfect knowledge and model structure (Schmidt et al., 2011).
To illustrate the uncertainty propagation in SOC projections caused by initialization and parameterization and to understand what correlates with the change in the patterns of projection uncertainty, we used the Agriculture Production System sIMulator (APSIM) (Keating et al., 2003;Wang et al., 2002;Holzworth et al., 2014) together with data from 90 agricultural experiments at 26 sites across the Australian cereal-growing regions.The data include measurements of total SOC stock (0-30 cm), C input (i.e.amount of residue retention), crop yield, and records of management practices.The APSIM model uses a very similar SOM pool structure and first-order decay approach to simulate SOM dynamics to other common Earth system models (Smith et al., 1997;Friedlingstein et al., 2006;Thornton et al., 2007).We firstly conducted sensitivity analysis to identify the model parameters whose change impacted most on simulated SOC dynamics.We then used Bayesian optimization approach to derive the posterior joint distribution of the identified parameters that enabled best match between measured and observed SOC.These ensembles of parameters were used to run AP-SIM for each of the 90 experiments, and simulations were continued for a further 100 years after the end of the experiment to produce SOC projections for uncertainty analysis.We quantified the uncertainty in SOC projections induced by both initialization of SOC pools and parameterization of algorithms for simulation of process dynamics.While the uncertainty obviously increases with years of projections, we further hypothesized that it is also influenced by site-specific climate, soil and management conditions, in addition to the impact of model initialization and parameterization.We fur-ther investigated how the projection uncertainty can be quantified by using these drivers, so that future SOC projections can become more useful with attached and well-quantified uncertainties.

Study sites and data sets
Data from a total of 90 experimental plots located within 26 different sites (Fig. S1 in the Supplement) and compiled and described by Skjemstad and Spouncer (2003) were used in this study.The experimental duration of these trials ranged from 5 to 82 years; the experiments cover diverse climate, soil and agricultural management conditions and are representative of Australian cereal-growing regions (Table S1 in the Supplement).The data set included detailed records on crop sequence, crop yield, crop residue production (estimated according to harvest index) and agricultural management practices such as residue management (removal or retention) and fertilizer application over each year.SOC stock was determined for representative 0-30 cm soil samples at least at the beginning and end of the each experiment, with some experiments having as many as six temporal measurements.Other soil properties at the start of the experiment were also measured, including total nitrogen content, bulk density, clay content and pH, and were used to initialize the APSIM model.

The APSIM model
APSIM was developed to simulate biophysical process in agricultural systems, and has been comprehensively verified and used to study productivity, nutrient cycling and environmental impacts of farming systems as influenced by climate variability and management practice (Keating et al., 2003;Wang et al., 2002;Holzworth et al., 2014).APSIM simulates crop growth and soil processes on a daily time step in response to climate (i.e.temperature, rainfall, and radiation) and soil conditions (water availability, nutrient status, etc.).The model allows flexible specification of management options like crop and rotation type, tillage, residue management, fertilization and irrigation.The ability of APSIM to simulate SOC dynamics under different cropping and management practices has been verified (Probert et al., 1998;Luo Z. et al., 2011).
APSIM simulates the dynamics of both soil C and N stocks in each soil layer.Similar to other SOM models like RothC and Century, SOM in APSIM is divided into six conceptual pools (i.e.microbial biomass, humic organic matter and inert organic matter, together with three fresh organic matter pools; Fig. S2).Inert organic matter is considered to be nonsusceptible to decomposition, i.e. indecomposable, due to physicochemical and/or biological protections.The amount of inert organic C is initialized at the start of the simulation and dos not change during the simulation.The decomposition of other pools is treated as a first-order decay process modified by soil temperature, moisture and nitrogen availability (for fresh organic matter pool only), leading to the release of CO 2 to the atmosphere and transfer of the remaining decomposed C to other pools.Microbial carbon use efficiency (CUE), i.e. the efficiency of microbial community to assimilate the decomposed SOC, determines the fraction of decomposed C transferred to other pools.The model assumes a constant CUE for all C pools.The flow of N depends on the C : N ratio of the receiving pool.The C : N ratio of each pool is assumed to be constant through time.The decomposition of surface residues is modified by the degree of contact of the residue with soil (Thorburn et al., 2001).
The model requires values for initial SOC content, total soil N content, bulk density, and soil hydraulic parameters for each soil layer simulated.In the Skjemstad and Spouncer (2003) data set, measured values for SOC content, bulk density and total soil nitrogen content were provided for the 0-30 cm layer.For the deeper soil layers and hydraulic parameters in the whole soil profile, values from a measured soil profile nearest to the site were selected from the Agricultural Production Systems Research Unit (AP-SRU) reference sites soil database (http://www.asris.csiro.au/mapping/hyperdocs/APSRU/).Daily weather data (from 1889 to present) for each site including radiation, maximum and minimum temperatures, and rainfall were obtained from the SILO Patched Point Dataset (https://www.longpaddock.qld.gov.au/silo/).
The APSIM model was first set up for each experiment.Agricultural management including crops, residue management and fertilizer application was set according to available historical records.Crops were sown depending on rainfall (> 20 mm in five successive days) and soil water content (90 % of saturation water content in the top 20 cm soil).Crop cultivars were assigned according to sowing date, i.e. the earlier the sowing date, the later the maturity type of the crop cultivar.For simplification, three cultivars for each crop representing early, middle and later maturity cultivars were selected from the default cultivars in the files released with the APSIM model.For pasture, however, there was no record on the species and cultivar.Here, perennial lucerne (Medicago sativa, a commonly used species in Australian pasture) was used to represent pasture and only one cultivar -Trifectawas used in the simulation.Lucerne was sown and removed after harvesting and before sowing of annual crops in the corresponding rotations, respectively.Harvest to the height of 10 cm was assumed each time lucerne reached the flowering stage to mimic possible grazing and/or haying.
In the experiments included in this study, C from assimilation of crop growth was the only source of C input to the soil.In the APSIM model, crop growth is simulated using light interception and radiation use efficiency, modified by water and nitrogen supply.In order to achieve credible simulation of crop growth, plant available water capacity (PAWC) of the soil was adjusted.This adjusted PAWC at each site was used throughout the simulations.Despite the reliability of the APSIM model to simulate crop growth (both belowground and aboveground), we did not use the simulated aboveground C input during the simulation.Alternatively, the recorded aboveground C input (as crop residue) was manually incorporated into the model at the time of crop harvesting, whilst the simulated crop residue was removed.This manipulation eliminated the effect of imperfect match of modelled with observed crop residue on SOC dynamics.

Sensitivity analysis of SOC dynamics
A total of eight parameters (Table S2) that directly link to the SOC dynamics in the model were selected for sensitivity analysis in order to identify the most important ones regulating SOC dynamics.One model input for model initialization, i.e. the fraction of inert organic carbon in the total SOC at the start of the simulation (finert), was also included in the sensitivity analysis, due to a lack of observed data of finert and its great effect on simulated soil C changes.To inspect the response of simulated SOC to variations of those parameters (finert was also called as a parameter for convenience hereafter), a univariate local sensitivity analysis was conducted by looking at the impact of one parameter at a time and fixing all other parameters.As the purpose was to identify the most influential parameter(s), a continuous wheat system with 100 % residue retention (the dominant crop in the studied experiments; see Table S1) and a nitrogen application of 200 kg N ha −1 yr −1 were used and simulated for 100 years.The default model parameters were first used (Table S2), and then each parameter was sequentially increased by 10 % of its default value.For each parameter, the sensitivity function (S i ) was calculated to represent the sensitivity of model output, y, (i.e. total 0-30 cm SOC stock) to changes in a single parameter, θ i (Soetaert and Herman, 2008): where θ i is the default parameter value, y| θ i is the model output using θ i , θ * i is the altered parameter value (increased by 10 %) and y| θ * i is the model output using θ * i .Finally, the importance index of the ith parameter (I i ), i.e. the overall sensitivity of the output with respect to this parameter, was calculated by summarizing the sensitivities for the 100-year outputs (n = 100): where S ij is the sensitivity function for parameter i at the j th year of the n (n = 100) years of each simulation.The greater the magnitude of I is, the more sensitive the model output was to the parameter (Soetaert and Herman, 2008).The importance indices were compared among the nine parameters, and the most important parameters were identified and optimized to obtain the best agreement between simulated and observed SOC dynamics for each of the 90 experiments.As the relative importance of those parameters was independent of soil and climate conditions, the typical soil and climate at Wagga Wagga (a major cropping area in Australia, and one of the 26 sites used in the main text), New South Wales, Australia, were selected to conduct above analyses.

Model optimization
The differential evolution (DE) algorithm (which belongs to the class of genetic algorithms) was used to optimize the most influential parameters identified.The optimization was performed in R 3.0.3using the DEoptim function in the "DEoptim" package (Mullen et al., 2011).DE is a global optimization algorithm for continuous numerical minimization problems, which use biology-inspired operations of crossover, mutation, and selection on population in order to minimize an objective function over the course of successive generations.
To use DE, each parameter was first assumed to exhibit a uniform distribution bounded within a range (i.e. the prior distribution; see Table S2) that was biologically and physically possible based on previous knowledge about the process, thereby eliminating solutions in conflict with prior knowledge.The optimization performed a quasi-random walk through the multi-dimensional parameter space to find the parameter set that caused the model to generate the best match between predicted and observed SOC.The "best match" was defined as the model output that minimized the criteria selected for model evaluation (Table S3).Seven criteria that are commonly used in the literature were selected to assess the possible effects of criterion selection on modelling results.Using each criterion, the optimization was conducted 100 times (i.e. 100 ensembles of initial parameter values through quasi-random walk), which generated 100 ensembles of parameters (i.e. the joint posterior distribution of the most influential parameters), giving simulation results with approximately equally good matches to the observed data.Consequently, 700 ensembles of parameters (from using seven criteria) for each experiment were produced.The optimizing procedure and related simulations were operated on Bragg and Dell CPUs of CSIRO clusters.
However, the required computing time (∼ 2 days for one experiment and one selection criterion using 100 computer cores) has posed a significant challenge even using the highperformance computing clusters (Bragg and Dell CPUs) for this multi-parameter optimization of the process-oriented APSIM model.To complete all optimizations using seven criteria for the 90 experiments, a run time of 4 months was expected assuming that 1000 cores could be continuously available on the clusters.For this reason, the global optimization DE was only applied for two sites, i.e.Brigalow and Tarlee, providing two cases of DE optimization as compared to an alternative and faster Bayesian sampling approach as described below.
For all the experiments, a Bayesian sampling approach was substituted for the DE optimization in order to complete the work within a reasonable time but without much sacrificing of model performance.The APSIM model was run for each experiment for 100 000 times using 100 000 ensembles of parameters that were randomly sampled from their prior distributions.The best 100 ensembles of parameters were selected as their posterior distributions through using each criterion listed in Table S3.At Brigalow and Tarlee, the distributions of parameters "optimized" through this Bayesian sampling approach were compared with those optimized through DE optimization.The identified parameter ensembles by Bayesian sampling approach were referred to as "optimized parameters" in the following text and used to assess the uncertainty in projected SOC.

Uncertainty in projected SOC
After obtaining the 700 ensembles of optimized parameters (i.e. after "optimization period"), the APSIM model was run continuously from the start to the end of each experiment and then for an additional 100 years after the end of each experiment using each parameter set (i.e.700 simulations for each experiment).For the last 100-year simulations (i.e.projection period), a continuous wheat system was assumed together with 100 % residue retention, which is the same as that used in sensitivity analysis.Carbon input through crop residue retention was expected to be an important factor regulating SOC dynamics in the projection period.As residue (or biomass) production is dominantly controlled by fertilizer application rates under natural rainfall condition at each site, scenarios with nitrogen application rates ranging from 0 to 300 kg N ha −1 yr −1 with increments of 20 kg N ha −1 yr −1 were modelled.These scenarios made it possible to mimic different management practices that influence C input to the soil and to assess its impact on the uncertainty of simulated SOC induced by model initialization and parameterization.
Climate data from the start year of each experiment through to 2013 were used for the corresponding simulation period.For all years from 2014 onwards, the corresponding years of the latest historical climate data were used.For example, for the possible simulations from 2014 to 2104 (91 years), the historic climate data of 91 years from 1923 to 2013 were used.As we focused on the potential uncertainty induced by model parameterization and initialization, we did not consider the uncertainty related to climate change.
SOC content in the 0-30 cm soil layer was output at the start of projection (excluding the optimization period) and at the end of each year of projection (C i ).For the ith year of projection, the mean (M SOCi ) of C i of the 700 estimates was calculated, and the range (R SOCi ) of the 95 % confidence interval was calculated as the difference between 97.5th and 2.5th percentile of the 700 estimates.Then, the percentage uncertainty (U P i ) for that year of projection was estimated based on half of the R SOCi divided by the M SOCi : × 100 %, i = 1, 2, 3, . .., 100. (3)

Attributes controlling the variability in the uncertainty
After estimating U P , we further addressed the following question: how and why does the uncertainty (i.e.U P ) in projected SOC change across space and time?We hypothesized that U P is associated with the management in terms of residue C inputs.At the same time, we assumed that the detailed relationship between U P and C inputs is different not only across experiments but also across time periods of the projection.As the hierarchy of the relationships (i.e.individual-level C inputs group in experiments and time periods of the projection), a hierarchical regression model, also called a multilevel model (Gelman and Hill, 2006), was fitted to estimate U P i (y i ) on C input (x i ), applied to the J = 90 experiments and K = 100 time periods of projection.The multilevel model was written as a data-level model (the predicted U P i belonging to experiment j with k years of projection), allowing the model coefficients (α and β) to vary by experiment (j = 1, . .., J ) and time period of projection (k = 1, . .., K) (Gelman and Hill, 2006): and a decomposition of its intercepts and slopes into terms for experiment, the time period of projection and their interaction, and models for variation, where was the 2 × 2 covariance matrix representing the variation in the intercepts and slopes in the population of groups (experiments and time periods of projection).In essence, there is a separate regression model for each experiment and time period combination with the coefficients estimated by the weighted average of pooled (which do not consider groups) and unpooled (which consider each group separately) estimates, i.e. partial pooling.This hierarchical structure of the model allows the assessment of the variation in individual-level coefficients across groups and accounting for group-level variation in the uncertainty for individuallevel coefficients.
To assess the variation in individual-level coefficients (α expt j and β expt j ) across different experiments, a classic linear regression was conducted to identify the effects of different sources of variation.At the experiment level, we assumed that two groups of attributes influence α expt j and β expt j : (1) uncertainty in model parameters, i.e. the three optimized parameters based on experiment-specific data set, and (2) climate including mean annual rainfall and temperature, which are predominant factors controlling SOC dynamics during model optimization as well as during projection.The generalized variance (GV) was calculated as an indicator of the overall variation in model parameters, which is defined as the determinant of the variance-covariance matrix of the three parameters and is a scalar measure of overall multidimensional scatter.The two groups of attributes including all interactions were selected through a stepwise regression model selection by the Akaike information criterion.Before fitting the model, GV was logarithmically transformed to satisfy additivity and linearity assumptions and then centred by subtracting the mean of the data, and rainfall and temperature were also centred.For coefficients over the time spans of projection (α year k and β year k ), their relationship with the time span of projection were presented.All the statistical analyses including the multilevel model fitting were conducted using R software version 3.0.3(R Core Team, 2013).

Sensitivity analysis and model performance
Three parameters were identified as most influential on simulated SOC (Fig. S3).Microbial carbon use efficiency (CUE) had the biggest impact.This highlights the key role of microbial process to control SOM decomposition, and the need for better capturing the dynamics and impact of microbial process in SOM models (Allison et al., 2010;Singh et al., 2010;Sinsabaugh et al., 2013;Xu et al., 2014).As CUE was treated as a constant in most SOM models, a framework is needed to incorporate microbial data (e.g.community, activity, and their responses and feedbacks to biotic and abiotic factors) into SOM models to provide robust estimations and predictions.Potential decomposition rate constant of humic organic matter (k hum , day −1 ) ranked the second, followed by the fraction of the humic carbon that is recalcitrant to decomposition (finert).This result further indicates the importance of better www.biogeosciences.net/12/4373/2015/Biogeosciences, 12, 4373-4383, 2015 quantifying the decomposability of the heterogeneous SOM (Schmidt et al., 2011;Sierra et al., 2011).It should be noted that the actual decomposition rate is simulated through modifying k hum by a series of biotic and abiotic variables at different spatiotemporal scales, and different models simulate the responses differently (Todd-Brown et al, 2013;Exbrayat et al., 2014).Although we did not quantify the relative importance of these modifiers (e.g.soil moisture ad temperature), the results indicated that k hum has to be constrained, implying the importance of determining how it responses to environmental factors.The wide distributions of CUE, k hum and finert parameters (derived by constraining the model against the measurement data, Fig. 1b) imply deficiencies in our understanding of the microbial community and its activity and how they change with environmental conditions to modulate the SOM decomposition processes.
Our optimization procedure enabled accurate simulation of SOC change during the optimization period (Fig. 1a) using distinct ensembles of model parameters for each experiment (Fig. 1b).Pooling together all data of the 90 experiments, the modelled average SOC of the 700 simulations could explain 96 % (P < 0.001) of the variance in observed SOC (Fig. 1a).
For each experiment, model performance was nearly identical (Fig. 1a) when the simulations (using different parameter sets) were intercompared.At the Tarlee site (Fig. 2a), for example, the RMSE between modelled and observed SOC ranged from 0.44 to 0.52 t ha −1 , compared with the range of 3.11 to 3.12 t ha −1 at the Brigalow site (Fig. 2b).This high level of consistency highlights significant equifinality, i.e. different parameter ensembles leading to similar simulation results (Figs. 1b, 2c and d), in process-based SOM models, which must be addressed in modelling studies aimed at enhanced process understanding and hypothesis testing (Tang and Zhuang, 2008;Luo Y. et al., 2011).

Uncertainty in SOC projections
The accurate simulations of past SOC, however, do not guarantee convergent projections beyond the model optimization period.In contrast, running the model with the same parameter ensembles generated very divergent future projections (Fig. 2a and b), indicating significant uncertainty propagation with time of projection (Luo Y. et al., 2011;Tang and Zhuang, 2008).Furthermore, the uncertainty is also related to management in terms of C input level and site conditions.At Brigalow (Fig. 2b), for example, the 95 % confidence interval of projected SOC under optimal N input (i.e.no N stress for crops) ranged from 37 to 56 t ha −1 10 years after the model optimization period, which increased to 26-68 t ha −1 for the projected SOC after 50 years.Under the low N input scenario (0 kg N ha −1 ), the uncertainty was smaller.At Tarlee (Fig. 2a), the uncertainty propagation followed a similar pattern to that at Brigalow, but the uncertainty under the low N input scenario was much smaller.At Brigalow, in addition, we found that the choice of criterion (objective functions) in- fluenced the distributions of the derived parameters (Fig. 2d) because a specific criterion only focuses on a specific aspect (e.g.mean or variance) of the data and the model results, of which the consequence for SOC simulations (e.g. the bifurcation pattern of projected SOC shown in Fig. 2b) ought to be carefully considered in future studies.
It is important to note that the posterior distributions of model parameters were apparently different across experiments (Fig. 1b, c and d, and S4), confirming that model parameters are sensitive to the data constraining the model (Keenan et al., 2012;Hararuk et al., 2014;Luo et al., 2015) Our results indicate that CUE was likely higher for the site with a longer cultivation history (the Tarlee site) than for the newly cleared site (the Brigalow site, Fig. 2c vs. 2d), implying the potential importance of land use history for constraining model parameters such as microbial carbon use efficiency because land use history has a direct effect on the quantity and quality of carbon input as well as on soil properties.However, such impact needs further confirmation with more data.The distributions of the optimized model parameters were also influenced by the choice of criteria to evaluate model performance (Figs.2d, S5).The differences in parameter distributions subsequently impact on the SOC projections as shown in Fig. 2b, albeit with near-identical model performance in simulating historical SOC.In addition, finert and k hum were positively related (Fig. 2c and d), implying the importance of the interactions and/or feedback between different C pools and their impacts on soil C projection.These highlight the need for (1) improving the science for capturing process interactions in the model such as the role of microbial processes and conceptualization of hetero-Biogeosciences, 12, 4373-4383, 2015 geneous C pools and their transformation (Manzoni et al., 2012;Luo et al., 2014), (2) conducting model optimization conditioned on all observed data from experiments together with Bayesian inference technique, and (3) quantifying uncertainty in SOC projections with ensemble model simulations (Post et al., 2008;Weng and Luo, 2011;Xia et al., 2013;Hararuk et al., 2014;Luo et al., 2015).
If a continuous wheat system was practiced for 100 years after the end of each experiment at the 26 sites, optimal N management was predicted to result in an average increase in SOC (Fig. 3a), whereas a SOC decline would occur under zero N input (Fig. 3b).The amount of potential SOC change depends on not only the management level (N input) and the climate and soil conditions that determine the potential productivity of crops but also the initial SOC level at the start of the projections.Across the 90 experiments, the percentage uncertainty in the SOC projections ranged from 2 to 140 % with an average of 53 % under optimal N management (Fig. 3c), and from 0.8 to 108 % with an average of 40 % under zero N input (Fig. 3d).From applying this result to Australia's cereal-growing regions, the simulated potential SOC stock of ∼ 7.5 Pg (Luo et al., 2013) could be subject to 53 % uncertainty under no N deficiency and 100 % residue retention.

Attributes controlling the variability in the uncertainty
The uncertainty propagation with time of prediction and across experiments could be explained using a linear model by linking the percentage uncertainty (U P ) to the C input from crop residue (C R ), i.e.U P = α + β C R .However, both α and β changed significantly across experiments (Fig. 4a) and years of projections (Fig. 4b), and were also impacted by their interactions.Across the time periods of projection, the uncertainty increased with the number of years for projection, reflected by the linear increase in α (model intercepts) and asymptotic increase in β (model slope, Fig. 4b).The asymptotic increase in β (model slope) also implies that the relative contribution of C input to prediction uncertainty reduces with time.Across experiments, there was a marked variation in the effect of C input on U P , indicating impact of site-specific conditions (e.g.climate and soil as described later).Across sites and years of projections, the majority of positive β implies increased uncertainty in SOC projections with increasing C input, which has not been properly addressed in previous modelling studies (Joos et al., 2001;Jones et al., 2005;Smith et al., 2005;Ogle et al., 2010).The fate of C input has a direct effect on the amount of soil C.   The variance in model parameters (GV) across experiments had a major effect on the intercepts (positive at P < 0.001) and slopes (positive at P < 0.001) of the regression model linking U P to C input (Table 1).As GV was logarithmically transformed when fitting the model, the increase in U P with GV was exponential across experiments.This result highlights the crucial role to improve the representation of the sensitive microbial processes (Zhou et al., 2012;Xu et al., 2014) and the heterogeneous SOM (Sierra et al., 2011) in biogeochemical SOM models, and to constrain the space of relevant model parameters.For example, we assumed a relatively wide range of CUE (0.2-0.8) as the prior information for the Bayesian optimization.Sinsabaugh et al. (2013) suggested that CUE prediction should consider resource composition, stoichiometry constraints and biomass composition, as well as environmental drivers.A more informative prior of CUE could help reduce the uncertainty in soil C projections.
Rainfall and temperature, together with their interaction, had a significant impact on SOC projection uncertainty through their impact on the fitted model intercepts across experiments (Table 1).α expt j increased with temperature, but tended to decrease with rainfall, implying increased uncertainty in SOC projection under future warming and drying conditions.Based on the results, the uncertainty in projected SOC will be increased by 4.95 % if average temperature is increased by 1 • C under global warming.For the slopes β expt j , rainfall and its interaction with GV had a significant negative effect.These effects may reflect the impact of rainfall on both primary productivity (and thus C input) and soil moisture conditions (and thus microbial activity and decomposition rate of SOC), emphasizing the importance of understanding the interactions between soil processes and their re-sponses to external drivers and management such as temperature and rainfall (Davidson and Janssens, 2006;Carvalhais et al., 2014).

Conclusions
Our results demonstrate that great uncertainty exists in soil C projections from process-based SOM models, due to deficiency in model initialization and parameterization in capturing the process interactions, such as microbial C use efficiency and its drivers, as well as a lack of detailed information to initialize the model, e.g. the heterogeneous SOM with different decomposability.The prediction uncertainty propagates with extended years of projections and C input into soil.It is also influenced by site-specific climate (temperature and rainfall) and soil conditions together with management inputs, which determine both the C input (through primary productivity) and the SOM decomposition processes.The results also suggest that C projection into warming and drying future climate will be subject to even greater uncertainty.For agricultural land uses, uncertainty caused by management practices has to be carefully considered due to its impact on microbial activity and subsequent projected SOC.For any future predictions of SOC change, ensemble simulations conditioned on total observed data sets together with a Bayesian inference technique should be used in order to quantify the uncertainties in modelling results.Based on our results, future improvement in SOM modelling should focus on how the microbial community and its carbon use efficiency change in response to environmental changes, better quantification of heterogeneous SOM and the effects of its change on total soil C turnover.
The Supplement related to this article is available online at doi:10.5194/bg-12-4373-2015-supplement.
Author contributions.Z. Luo collected data, ran simulations, and performed data analysis; Z. Luo, E. Wang, J. A. Baldock designed the study; H. Zheng and Q. Shao were involved in statistical analysis; Z. Luo, E. Wang and O. J. Sun wrote the paper.All authors discussed the results and commented on the manuscript.

Figure 1 .
Figure 1.Model performance in simulating soil organic carbon (SOC) dynamics (a) and the corresponding optimized model parameters (b) across the 90 studied experiments.Circles and bars (a) indicate the average and 95 % confidence interval of the simulations for each experiment using different parameter ensembles.Red and blue symbols in (a) highlight the data at Tarlee and Brigalow, respectively, corresponding to the data in Fig. 2. Dashed line is the 1 : 1 line in (a).The parameter ensembles at Tarlee and Brigalow are highlighted in (b).See Fig. 2 for the means of the coloured symbols in (b), showing the different ranges of optimized fraction of inert organic carbon (finert).

Figure 2 .
Figure 2. Projected soil organic carbon dynamics at two case sites Tarlee (a) and Brigalow (b) and the correspondingly used parameter ensembles (c and d).Black symbols show the observations.Seven criteria (RMSE, MAE, pMAE, IoA, rIoA, NSE and rNSE; see Table 3 in the Supplement for details) are used to derive the posterior joint distribution of model parameters (CUE, k hum and finert).CUE, microbial carbon use efficiency; k hum , the potential decomposition rate of humic organic carbon; finert, the fraction of inert organic carbon.

Figure 3 .
Figure 3. Projected SOC (a and b) and its percentage uncertainty (c and d) under high (a and c) and low (b and d) carbon input scenarios after 100-year simulations in 90 experiments across 26 sites.Concentric circles show the different experiments at the same site.The sizes of the circles correspond to the projected average of SOC content (a and b) and the corresponding percentage uncertainty (c and d).Blue and red circles indicate that the average of the 700 simulations is increased and decreased, respectively, compared with the SOC content at the start of the projection.Blue and red sectors of the circles in (c) and (d) indicate the fraction of 700 bootstrapping simulations that shows an increase and a decrease in the projected SOC, respectively, compared with the SOC content at the start of the projection period.
0.21 * * * * * * P < 0.001; * * P < 0.01; * P < 0.05; * * * * P < 0.1.a GV, generalized variance of the identified three model parameters including microbial carbon use efficiency, decomposition rate of humic organic carbon and the fraction of inert organic carbon; R, the annual average rainfall; T , the annual average temperature.GV was logarithmically transformed and centred, and R and T were also centred when fitting the model.

Figure 4 .
Figure 4. Coefficients (estimate ± SD) for the regression model: U P = α+β C R .The model is fitted to estimate the effects of carbon input (C R ) on the percentage uncertainty (U P ) in soil organic carbon projections, applied to 90 experiments (a) and 100 time spans of projection (b).α and β show the data-level coefficients (i.e.averaging over experiments (a) and time spans of projection (b)) and σ represents model error.In (a), experiments are sorted according to α expt j .The coefficients at the experiment × time span level are not shown.See more details in the "Materials and methods" for the regression model.

Table 1 .
The effects of experiment-specific variance of model parameters and climate on individual-level coefficients (i.e.α