Sensitivity analysis and calibration of a soil carbon model ( SoilGen 2 ) in two contrasting loess forest soils

To accurately estimate past terrestrial carbon pools is the key to understanding the global carbon cycle and its relationship with the climate system. SoilGen2 is a useful tool to obtain aspects of soil properties (including carbon content) by simulating soil formation processes; thus it offers an opportunity for both past soil carbon pool reconstruction and future carbon pool prediction. In order to apply it to various environmental conditions, parameters related to carbon cycle process in SoilGen2 are calibrated based on six soil pedons from two typical loess deposition regions (Belgium and China). Sensitivity analysis using the Morris method shows that decomposition rate of humus ( kHUM), fraction of incoming plant material as leaf litter (fr ecto) and decomposition rate of resistant plant material ( kRPM) are the three most sensitive parameters that would cause the greatest uncertainty in simulated change of soil organic carbon in both regions. According to the principle of minimizing the difference between simulated and measured organic carbon by comparing quality indices, the suited values of kHUM , frecto andkRPM in the model are deduced step by step and validated for independent soil pedons. The difference of calibrated parameters between Belgium and China may be attributed to their different vegetation types and climate conditions. This calibrated model allows more accurate simulation of carbon change in the whole pedon and has potential for future modeling of carbon cycle over long timescales.


Introduction
The terrestrial ecosystem is one of the essential parts of the global carbon cycle.Significant variations of terrestrial carbon pool at geological timescales have played an important role in past atmospheric CO 2 concentration change (Falkowski et al., 2000;Post et al., 1990).The soil carbon pool is much larger than the biotic pool (Lal, 2004) and accounts for about two-thirds of the terrestrial carbon pool.Thus quantitative estimation of the soil carbon pool is the key to revealing the mechanism of past terrestrial carbon cycle and narrows the uncertainties in the global carbon cycle inventory.However, because only parts of the carbon pool are preserved in sediments, past soil carbon pool reconstruction is difficult by direct measurement.Modeling approaches then become the potential option for accurate estimation.
Currently, with the development of soil carbon models, quantitative simulation of soil carbon storage has been widely done, but mostly focuses on modern processes and aims to predict future atmospheric CO 2 concentration change (Coleman et al., 1997;Jensen et al., 1997;Kelly et al., 1997;Li et al., 1997).Changes of past soil carbon pools over long timescales have to account for changes in soil properties (e.g., particle size, pH) due to soil formation processes, which are seldom included in existing soil carbon models (Finke, 2012;Finke and Hutson, 2008;Mermut et al., 2000).Information on past soil formation factors for different regions is unavailable (Finke, 2012;Sauer et al., 2012), and few models consider the effect of all soil formation factors (Jenny, 1961;e.g., climate, organisms, relief, parent material and time) in simulation of soil formation (Kirkby, 1977; Y. Y. Yu et al.: Sensitivity analysis and calibration of a soil carbon model (SoilGen2) Minasny andMcBratney, 1999, 2001;Minasny et al., 2008;Parton et al., 1987;Salvador-Blanes et al., 2007).
SoilGen2, developed by Finke (Finke, 2012;Finke and Hutson, 2008), is a first attempt to reconstruct most aspects of soil evolution by taking all soil formation factors into account.The advantage of the model is that it simulates the organic and inorganic carbon cycle simultaneously and reveals the influences on carbon pool by other soil processes at long time scales.The model has been validated and applied in European soils developed over 15 000 yr from loess parent materials (Finke, 2012;Finke and Hutson, 2008), and the results show that clear sensitivity and plausible response of this model to the "climate", "organisms" and "relief" factors of soil formation exist.It also has been confirmed that reconstructions of realistic initial status of soil profiles (including carbon and other element contents) can be evaluated through simulating soil formation by SoilGen2 (Sauer et al., 2012).Therefore, the model offers an opportunity to reconstruct the past soil carbon cycle.
Because the verification and application of SoilGen2 is still at its preliminary stage, only parts of the soil processes included in the model have been calibrated (e.g., calcite leaching and clay migration) (Finke, 2012;Finke and Hutson, 2008).No calibration on parameters related to organic carbon (OC) cycle has been done yet.This work is necessary, and this activity should be preceded by an analysis of such a model to determine its most sensitive parameters (Skjemstad et al., 2004).
In this study, soil pedons from two typical loess deposition regions (Belgium and China) with distinct climate conditions are selected to calibrate OC cycle process in SoilGen2.Loess deposits have been continuously and widely deposited in Eurasia since 22 Myr ago (Guo et al., 2002;Kukla, 1987;Liu, 1985).More than 400 paleosols were developed in the loess-soil sequences in China (Guo et al., 2002), and these provide the best record for reconstruction of past carbon cycle through modeling soil formation processes in future studies.
In summary, the objectives of this study are as follows: (1) to use sensitivity analysis to assess which parameters in SoilGen2 potentially cause the greatest uncertainty in calculated change in soil OC in Belgium and China; and (2) to do calibration and validation of the parameters related to OC cycle in Belgian and Chinese soil pedons.We focus on forest vegetation on loess soils in this study.

Modeling soil carbon change with SoilGen2
In essence, SoilGen2 is an extended solute transport model solving the Richards equation for unsaturated water flow and the convection-dispersion equation for solute transport.Additionally, heat flow is calculated to estimate soil  temperature, which allows evaluation of the effect of soil temperature change on values of chemical constants, mineralization of OM (organic matter) and simulation of the effect of frozen soils on water flow.It simulates various aspects of pedogenesis including, e.g., OM accumulation, clay migration and CaCO 3 leaching.For detailed model description, refer to Finke and Hutson (2008).This article focuses on the description of the OC cycle in SoilGen2, which interacts with other soil formation processes (e.g., clay migration, (de-)calcification and bioturbation) through the change of soil physical properties (porosity and texture), hydraulic and thermic conductivity, and associated water and heat flow in the soil profile.One clear feedback mechanism between soil formation processes and the OC-cycle is that changes in water content, clay content and temperature due to pedogenetic changes in soil properties affect the degradation rates of organic matter.Basically, the soil profile is divided into a number of compartments with equal thickness, and the routines of Soil-Gen2 operate on every compartment separately.Therefore, changes of OC in the profile are modified by different water content, clay content and soil temperature in corresponding compartments.One key soil parameter for converting OC fractions to mass per unit area, bulk density, is calculated by the model by division between of the simulated mass of the solid phase and (fixed) volume per compartment.
Figure 1 shows the OC-cycle process modeled in Soil-Gen2.Vegetation provides dead plant material (leaf and root litter) as model input, which contains Ca 2+ , Mg 2+ , K + , Na + , Al 3+ , Cl − , SO 2− 4 , HCO − 3 and CO 2− 3 previously taken up from the soil solution via the transpiration stream.These ions follow the carbon decomposition pathway described hereunder dependent on vegetation types.Four vegetation types (grass/scrub, conifers, deciduous wood and agriculture/barley) are identified in SoilGen2, each having a unique root distribution pattern, associated water ion-uptake and target ion composition in living biomass (Finke, 2012: Table 2).Decomposition rates are considered invariant with respect to the various ion species, which is a simplification of the true system.
Dead plant material is distributed over root and leaf litter with a vegetation-dependent fraction fr ecto , and the root litter input is distributed over the soil depth such that it reflects the root density distribution.These litter inputs are then split and added to the decomposable (DPM) and resistant plant material (RPM) pools using the (vegetation-dependent) DPM/RPM factor.These plant material pools, existing in both the ectorganic layer and the individual mineral soil layers (endorganic layer), gradually decompose and mineralize while being incorporated and redistributed in the soil by bioturbation, which is modeled as an incomplete mixing process (Finke and Hutson, 2008).Each soil layer that is subject to bioturbation contributes a depth-dependent input mass fraction to a bioturbation pool, which is then mixed vertically.
Further decomposition of OM is modeled according to the concepts of the RothC26.3 model (Coleman and Jenkinson, 2005;Jenkinson and Coleman, 1994) using degradation rates that are modified as a function of soil temperature and moisture conditions.The mineralization process finally produces CO 2 and releases cations and anions (Finke and Hutson, 2008: Fig. 1) into the soil solution.The decomposition rates of DPM, RPM, biomass (BIO) and humified OM (HUM) are considered similar for all vegetation types.This is also assumed for the scaling factor (scalfac) and the fractionation parameter BIO/HUM; scalfac is the constant (1.67) used in the equation to set the CO 2 /(BIO+HUM) ratio as in RothC 26.3 model, and the ratio is influenced by the clay content of the soil.In SoilGen2, the constant is transformed into a parameter that can be calibrated, and the ratio varies per soil compartment with the variation of clay content in the profile.
Figure 1 shows there are eight parameters that describe the carbon decomposition pathways: fr ecto , DPM/RPM (both for deciduous forest), k DPM , k RPM , scalfac, BIO/HUM, k BIO and k HUM .The relative importance of these parameters will be tested by sensitivity analysis and then calibrated in this study.

Belgian loess soil under permanent deciduous forest
This study concerns soils at three topographic positions in the Sonian Forest near Brussels, Belgium (50 • 46 31 N, 4 • 24 9 E), developed in loess deposited during the Weichselian glaciation.The three pedons are located at mutual distances of less than 100 m, but extensive research revealed a clear relation between slope exposition and decalcification depth (Langohr and Sanders, 1985), which was confirmed by model simulations (Finke, 2012).The loess cover is 2-4 m thick and overlies a dissected plateau of pre-Weichselian age in Tertiary clays that locally cause water stagnation, but not at the three plot sites.Langohr and Saunders (1985) proved that the landscape has hardly eroded in the last 20 000 yr. Annals of landowners from the 14th century onwards indicate that this area was never under agriculture, as it was used for hunting by the nobility at least from this time onwards.
Older reports indicate that it was a mixed beech/oak forest previously (Davis et al., 2003;Verbruggen et al., 1996).
Detailed analysis of the mineral soils of these pedons was reported in Finke (2012;Table 3).For this study, samples from the ectorganic litter layers were also taken and analyzed (Table 1) for later comparison with simulation results.Volume and dry weight were measured for bulk density estimation.The weight loss-on-ignition (LOI) method was used to determine OC.Because the Belgian soil has been decalcified and contains fairly low amounts of clay in the top, the bias of LOI method, overestimation of OC by ignoring loss of water from various clay minerals, calcite and gypsum, is negligible.

Chinese loess soil under secondary and artificial deciduous forest
This study also concerns three pedons at plateau position located in Ziwu Mountains (35-36  (Liu, 2007).However, since 1870s population moved out of the region and in 1970s a forest  (Gong et al., 2003;IUSS Working Group WRB, 2006), production forest (Robinia, 20 yr).
Samples from ectorganic litter layers were taken in the same way as in the Sonian forest of Belgium, while in mineral layers samples were taken at 5-10 cm depth intervals.Bulk densities were measured by volume and dry weight method.The potassium dichromate method, which is not sensitive to the high CaCO 3 content, is used for OC analysis (Table 1).

Model input data
Two types of inputs are included in SoilGen2: one is for boundary conditions (e.g., climate, litter input and bioturbation history); another is for initial conditions (e.g., soil properties, typical year weather pattern).A full simulation of soil from incipient stages of soil formation up to today (more than 10 000 yr) is ideal, but it takes a long run time, which would render sensitivity analysis and calibration unfeasible.As former tests showed that amounts of OC in soil pedons would become stable after a 300-yr simulation in case of invariant climate and vegetation conditions (Finke and Hutson, 2008), all the tests in this study have a temporal extent of 1000 yr so that effects of initial values of OC are eliminated and effects of soil properties are minimized.

Inputs for Belgium
Climate and weather data were taken from the nearby weather station of Uccle (near Brussels and at 5 km from the studied site).A typical year of daily rainfall and weekly potential evapotranspiration data was used for the whole simulation period with an annual rainfall sum of 849 mm and potential evapotranspiration of 649 mm.Over this year, the characteristics of precipitation (e.g., total amount of the rainfall, the number of days with rain, frequency of rainfall intensity and monthly distribution of rainfall) were similar to multi-year average conditions.The average January temperature was 3 • C and July temperature 18 • C.An annual litter input of 4.7 Mg C ha −1 yr −1 (A.Verstraeten, personal communication, Research Institute for Nature and Forest, Belgium, 2012) was taken.The bioturbation was assumed to be 8.2 Mg ha −1 yr −1 affecting the upper 30 cm of soil (Gobat et al., 2004).Initial physical and chemical properties of the soil pedons were only partly known from measurements (Finke, 2012: Table 3), and to obtain a complete set of initial soil properties we did the following: 1. Starting from the properties of the C-horizon, we simulated soil formation between 15 000 yr ago (end of loess deposition) to present.See Finke and Hutson (2008) and Finke (2012) for details concerning the modeling approach and inputs.
2. The simulated properties at present were taken as initial inputs for the various scenarios of following tests for the sensitivity analysis and calibration.However, simulated OC was re-initialized to 0.5 % OC throughout the pedon.
As the Belgian climate is a strongly leaching one with practically no dust additions, above reconstruction reflects the effect of leaching on soil properties (e.g., soil texture, calcite and hydraulic properties) over 15 000 yr and shows a more realistic initial condition.
Comparison between simulations and measurements (Finke, 2012: Table 5) showed that simulations could reproduce the A-E-Bt horizon sequence and also the World Reference Base (WRB) soil classifications based on available (non-morphological) measurements.Therefore these simulations were considered as suitable basis for the current study.

Inputs for China
Representative climate data were interpolated from nearby weather stations in China.The average January/July temperatures were −6/22 • C, −5/23 • C and −5/23 • C for LJB, ZW2 and ZW3, respectively, while annual rainfall and potential evapotranspiration were 482/1645 mm, 516/1582 mm and 519/1587 mm, based on inverse distance interpolation of 30-yr (1958-1988) monitoring data of 61 weather stations distributed over the Loess Plateau.A typical year of daily rainfall and weekly potential evapotranspiration was from monitoring data of Xifeng in 1978 through comparing the characteristics of yearly precipitation at three weather stations near these soil pedons.
Annual input of litter (Populus 4.5 Mg C ha −1 yr −1 , Robinia 4.4 Mg C ha −1 yr −1 ) was transformed from measured biomass data (including volumes of growing stock per unit area, net annual increments and removals) in ZW forest station (Zhang and Shangguan, 2005), based on the protocol developed by De Wit et al. (2006).The bioturbation was assumed to be 15.3-17.4Mg ha −1 yr −1 affecting the upper 70 and 100 cm of soil for Populus and Robinia ecosystems, respectively (Gobat et al., 2004).In addition, distribution of monthly litter input and roots were adjusted in the model according to observed data of Populus and Robinia ecosystems in Loess Plateau (Cao et al., 2006;Cui et al., 2003;Hu et al., 2010;Zhang et al., 2001).
The Chinese climate has a large precipitation deficit, and redistributions of calcite and clay towards greater depths are not as significant as in Belgium.Furthermore, dust deposition is not negligible in study region of China, and fertilizes the profile with additional calcite and changes the particle size at the top.Soil properties with uncertain dust addition cannot be simulated accurately to reconstruct the situation of 1000 yr ago, so we decided to assume that the current profile represents this situation fairly well, and initial physical and chemical properties of the soil pedons were from measurements of their parent material layers at the bottom of pedon (Table 1).

Sensitivity analysis method
Sensitivity analysis (SA) determines the response of selected model outputs to variations (within plausible bounds) of uncertain input parameters (Saltelli et al., 2000).Results of SA can be used to select and rank the most important parameters for calibration.Various SA methods have been developed (Saltelli et al., 2005).A choice of a particular method is based on a function of the number of parameters to be evaluated and the CPU-time per run.The number of parameters to be evaluated in the current study is eight, and a typical SoilGen-run for a 1000-yr period takes about 20 h CPU time.Under these circumstances, Saltelli et al. (2005) proposed four methods: Bayesian sensitivity analysis (Oakley and O'Hagan, 2004), fractional factorial designs (Campolongo et al., 2000), automated differentiation techniques (Grievank and Walter, 2008) and the Morris method (Morris, 1991).
Bayesian sensitivity analysis is more efficient than traditional Monte Carlo techniques but still requires substantial amounts of simulations and reprogramming of the Soil-Gen code.This technique is therefore considered beyond the scope of this study.Fractional factorial designs have the disadvantage that assumptions need be made on model behavior.Automated differentiation techniques also require substantial reprogramming of the model code and may lead to results only representing local areas in parameter space (Saltelli et al., 2005).The Morris method is feasible in terms of computing time, because it takes samples from levels rather than from distributions of parameters (which may be a drawback when such distributions are known, but this is not the case here).For the above reasons we chose to apply the latter method.
The Morris method is based on the principle that one factor (model parameter) is varied at a time over a certain number of levels in parameter space.Each variation, comprising two simulations, leads to a so-called elementary effect u i : where x is the parameter value, x i its imposed variation for factor i only ( x i = 0 for the other factors) and Y (x) the model result with parameter set x. Values for x are randomly chosen inside a plausible parameter value range [x i,low x i,high − x i ], and x i is either 0 or a predetermined multiple of 1/(p−1) with p the number of considered parameter levels (rescaled at range [0, 1]) for which p/2 elementary effects are computed.In this study we took p = 8, and fixed x i at 2 × 1/(8 − 1) on the [0, 1] rescaled range.The number 2 is an arbitrary choice in the Morris method (Morris, 1991); it determines what fraction of the plausible parameter range is covered by each elementary effect.The obtained elementary effects u i comprise a simple random sample, of which the mean µ i and standard deviation σ i are used to assess how important a factor is. Hereto, a graph is made displaying the position of a factor i in terms of µ * i (the average of |u i |) and σ i .If the value of µ * i is high, then there is a high linear effect of factor i; large values of σ i indicate either non-linear behavior of the model for factor i or non-additive behavior (relative to other parameter values).In this study we took eight parameter levels, resulting in four elementary effects, for each one of eight model parameters.Thus, 64 simulations (32 pairs) for a period of 1000 yr were done for a typical loess forest soil in Belgium (2.5 m depth with a vertical discretization of 5 cm at the Uccle Plateau location) and 64 more for a loess forest soil in China (1.5 m depth with the same vertical discretization at pedon of LJB).This comprised about 85-CPU-day simulation time on one core (less than 6 days on four quad-core PCs), which was considered feasible.The model output parameters considered were OC (ton ha −1 ) in ectorganic layers and OC (mass ton ha −1 and content %) in the mineral soil.Separate analyses were done for ectorganic layers and endorganic layers, because later calibration would focus on the vertical distribution of OC.

Calibration and validation approach
Calibration is the process of modifying the input parameters to a model until the output from the model matches an observed set of data.Various techniques also have been developed for model calibration, which differ in how parameter combinations are generated and how results are compared.In many cases, the modeler selects parameter combinations and evaluates results by expert judgment, in which case calibration is more or less a skill and may not detect the optimal parameter combinations.An alternative, often used technique is the minimization of an object function describing the deviations between measurements and simulations for various settings of parameters.The minimization process advises on optimal parameter combinations under the assumption that model outputs are differentiable with respect to the model parameters.
A well-known implementation is the PEST software (Doherty, 2004).Model runs are sequential, and the software decides if a new run with changed parameter settings is needed after results of a preceding run have been confronted to measurements by evaluating the object function.Another emerging method is the exploration of parameter space by a Markov chain Monte Carlo method.Model results are evaluated by calculating the posterior probability of the parameter set given the data, using the prior distribution of the parameters and a likelihood that expresses the correspondence between measurements and simulations.Of this Bayesian calibration method, various implementations exist, but even the most efficient ones (Vrugt et al., 2009) require numerous simulations.With time-consuming models such as Soil-Gen, convergence may take very long both in PEST and in Bayesian calibration.Finke (2012) used therefore an alternative approach in which various chosen sets of parameters were run with the model in parallel, confronting the model with measured data to quantify simulation accuracy and fitting a polynomial function predicting simulation accuracy as a function of parameter value.Analyzing the partial derivatives of this function, the position in parameter space with optimal simulation results was predicted.This approach may not find the true optimal parameter set and also depends on the choice of the evaluated parameter values.However, for reasons of runtime it was applied in Finke (2012).
The overall procedure of calibration is given in Fig. 2, following the principle of minimizing the difference between measured and simulated OC step by step.During the steps, parallel tests were firstly done for soil pedons by varying the most sensitive parameters identified in sensitivity analysis.Then the results were evaluated according to the quality indices described below.If there was still possibility to reduce the difference between measured and simulated OC by varying the same parameter, more parallel tests would be done in a sub-range of the parameter; otherwise, the next, less sensitive, parameter would be selected for tests.The process would be repeated, according to the order of parameter sensitivity, until no improvement was identified.In both studied regions, the difference of simulated and measured total carbon of the whole pedon was firstly minimized.In the second stage, the distribution of OC over ectorganic and endorganic layers was calibrated.The calibration aims at minimizing the number of runs by a quick convergence between simulated and measured OC.First, the quantitative relationships between changed OC and calibrated parameters revealed by sensitivity analysis were used to set limits to the variation range of the parameters.Second, a convergence criterion (< 5 %) was applied to a particular quality index during the calibration.Calibration for each parameter was stopped when the best result was less than 5 % better than the second best result (relative to the value of the quality index of the first simulation).
The number of soil pedons per region (3) was low because of the large input requirements of the model.We identified the stability of the calibration result by repeating the calibration three times per region.Each time, two of three soil pedons in each region were used to calibrate the model, while the third one was used for validation.
The indices used for assessing the quality of the C-module of SoilGen2 are the following: 1. root mean square deviation (RMSE) of total simulated and measured OC mass per mineral soil compartment: where f OCM and f OCS are measured and simulated OC mass, respectively, ρ M and ρ S measured and simulated bulk density (kg dm −3 ), and T the thickness of the k soil compartments (all equal to 50 mm).RMSE1 endo,pedon of the two pedons in each group of pedons in each region are averaged to obtain RMSE1 endo ; 2. RMSE of total simulated and measured OC mass in the ectorganic layer: where OC S is the simulated OC mass in ectorganic layers and n is the number of pedons (two per group of pedons in each region); 3. RMSE of total simulated and measured OC mass in the whole soil pedon: 4. mean difference (MD) of total simulated and measured OC mass in mineral soil: 5. MD of total simulated and measured OC mass in ectorganic layers: 8. dissimilarity (DIS) (Gower, 1971) of simulated and measured OC % in mineral soil: where OC % max and OC % min are the maximal and minimal value found in a particular pedon.DIS OC % is calculated by averaging over two pedons and varies between 0 (perfect) and 1 (very poor).DIS gives the absolute difference as a fraction of the observed differences and thus is dimensionless.Therefore, it allows to compare results over different output parameters, which is not possible with other indices.
Of these statistics, the first six express how well the total OC mass in the soil is simulated, while the last two express how well the vertical distribution of OC content over the pedon is simulated.The MD statistics indicate systematic under-or overestimation of simulated values relative to measurements, while the RMSE statistics focus on the absolute values of these differences.Both are ideally close to 0. The improvement of RMSE OCMS (IM RMSE OCMS ) was used as convergence criterion in the first stage of calibration, while IM RMSE ecto and IM RMSE endo were used in the second stage.

Sensitivity analysis
Table 2 gives the model parameters that were considered in the sensitivity analysis and the plausible range of these parameters.In Fig. 3, at a double-logarithmic scale, the µ * and σ of the elementary effects of the eight factors are shown.The µ values for all rate factors were negative, which was expected because these describe decomposition and positive values for x i , leading to higher values for k, and are expected to lead to lower amounts of remaining OC.
The OC mass in ectorganic layers and endorganic layers and the OC % in endorganic layers respond most sensitively to the decomposition rate of humus k HUM and less sensitive to the fraction of dead plant material entering the soil as leaf litter fr ecto , k RPM , and scalfac.The other factors show less sensitivity.Most responses are between the |µ| = SEM (stand error of the mean) and |µ| = 2*SEM lines indicating a fair confidence level.Most certain response of the sensitive factors is that of fr ecto , while the other factors are less certain.This may be caused by non-linear response or non-additive behavior (the model responds to interactions of factors).
The sensitivity order in the ectorganic layer is similar to that of the mineral soil, except for k RPM and fr ecto .k RPM is more important in the ectorganic layer; the rate modification due to moisture deficit is always equal to 1 in ectorganic layers while it can decrease in the mineral soil (Coleman and Jenkinson, 2005).This results in stronger modified k RPM in the ectorganic layer.Comparison (Table 2 and Fig. 3) of results for the Chinese and Belgian loess forest soils shows that the sensitivity orders of the factors follow the same pattern, irrespective of the differences in soil (the Belgian loess soil is strongly leached whereas the Chinese is not) and climate (a large precipitation surplus in Belgium and a large precipitation deficit in China).The values for µ * and σ differ (Fig. 3), which is higher in endorganic layers in Chinese soil pedons but lower in ectorganic layers.Nevertheless, the sensitivity order is the same, which means that the same model parameters could be selected for calibration in both soils: k HUM , fr ecto , k RPM and scalfac.

Calibration and validation
Based on 67 and 57 test simulations, three sets of calibrated parameters were obtained for Belgium and China, respectively.Table 3 gives the model parameters in calibrations tests.Results comparing simulated and measured OC by eight quality indices are given in Figs. 4 and 5 and Table 3, with lower absolute values of them indicating better simulated results.
In the stage of calibration (steps 2 and 3 in Table 3), tests (1bs and 1cs) with default values of parameters in SoilGen2 (step 1 in Table 3) showed that simulated total OC mass was lower than measured ones in both regions with larger deviations in Belgium than in China (Fig. 4d-f and j-l).The differences in total OC mass were minimized by decreasing the most sensitive parameter k HUM and the second sensitive parameter k RPM during the first calibration stage (Table 3).The IM RMSE OCMS became much lower than 5 %, when k HUM values were set as 0.0061-0.0065and 0.014-0.019,and k RPM values were 0.27-0.29 and 0.25-0.29 in Belgium and China, respectively.The RMSE OCMS resulting from these simulations served as target quality for the second stage of calibration.
The following calibrations for both regions were to adjust the distribution of OC in ectorganic (overestimation in stage 1) and endorganic (underestimation) layers.Based on tests in steps 4 in both regions with variation of fr ecto , it was revealed that the decrease of fr ecto could also increase the amount of OC mass in both layers (Table 3).In order to reduce this effect and the values of RMSE ecto and RMSE endo simultaneously, k HUM and k RPM were increased slightly (steps 5 and 6 in Table 3).The best results were finally obtained in Test21b 1, Test22b 2 and Test22b 3 in Belgium and Test18c 1, Test17c 2 and Test19c 3 in China with the following combination of parameters (Belgium, k HUM = 0.0065-0.0074,k RPM = 0.27-0.27and fr ecto = 0.30-0.38;China, k HUM = 0.016-0.019,k RPM = 0.30-0.30and fr ecto = 0.37-0.43),with the lowest values of IM RMSE ecto and IM RMSE endo at the same time while the RMSE OCMS remained close to the value reached in the first stage of the calibration.
Above tests were further confirmed as the best calibration results by comparing eight quality indices synthetically (Table 3 and Fig. 4), because the indices of them all belong to the lower ones in all tests of two regions.Figure 5 further shows that the simulated vertical distributions of OC content are also similar to measured ones visually.
The comparisons between measured and simulated OC of the validation pedons are given in Fig. 6.Two groups of OC mass are evenly distributed along 1 : 1 line (Fig. 6a), and RMSE OCMS values of pedons are all less than 12 % of the total OC mass.Furthermore, the consistency is also evident in vertical distributions of OC content (Fig. 6b-g).
This two-stage calibration approach was repeated three times, always calibrating with two pedons in a region while checking the quality with the third pedon.Results indicate that the calibration results are fairly stable and independent on the subset of two pedons chosen for the calibration.This suggests that the calibrated parameters are applicable to other soil pedons in similar climate conditions.

Comparison with former studies
Our results of sensitivity analysis are in accordance with former studies on RothC model in surface forest soils (Paul and Polglase, 2004;Paul et al., 2003), which indicated that change in soil carbon is particularly sensitive to the decomposition rates of HUM, RPM and BIO pools.Comparing that only the relative importance of the parameters was shown in former analysis (Paul and Polglase, 2004;Paul et al., 2003), a quantitative evaluation of their importance is given in our study and the especially significant sensitivity of k HUM is revealed.
The calibrated values of parameters (k HUM and k RPM ) in our study all fall into the logical range of former calibrations for RothC model (Shirato et al., 2004;Skjemstad et al., 2004;Todorovic et al., 2010), covering various climate conditions  (Jenkinson, 1990).The difference may be attributed to following aspects: firstly, decomposition in agricultural soils is faster than that in forest soils because of its lower lignin content in litter (Lambers et al., 1998) and more favorable micro-climate conditions for decomposition induced by human disturbance (Schlesinger and Andrews, 2000); secondly, carbon at deeper depth (1.5-2.5 m in our study) is older than that near the surface, indicating that it has a greater resistance to decomposition or that the environment at depth is less favorable for decomposition processes (Swift, 2001).
The calibrated fr ecto is lower than the default value (0.58) in SoilGen2 based on measured data (Kononova, 1975).In realistic soil carbon cycle process, part of litter carbon pool in ectorganic layer leaches to endorganic layers in the form of dissolved organic carbon (DOC).However, this process is not simulated in SoilGen2, while only little carbon is being exchanged between two layers by bioturbation in the model.Therefore, fr ecto , as the ratio of carbon pool in ectorganic layer to the total pool, was expected to be lower than the literature value as this decrease mimics the effect of DOCleaching.

Comparison between two regions
Although the orders of sensitivity for parameters are the same in two regions, the responses are less significant in Belgian soil pedons (Table 2).It led to corresponding larger ranges of parameters (k HUM , fr ecto and k RPM ) that varied during calibration in Belgium (Table 3), which shows slower decomposition rate of OC in Belgian soils.The differences may be driven by the following reasons.Firstly, litter chemical composition is one of the most important factors that affect decomposition of litter.Especially in late stage of decomposition for the formation of humus, lignin decomposition exerts the dominant control in soils (Berg and McClaugherty, 2008;Quideau et al., 2001), which is relatively resistant to decomposition (Lambers et al., 1998).Since the study area in Belgium is under beech forest while it is under poplar in China, higher lignin and holocellulose contents in the former ecosystem compared to the later (Coldwell and Delong, 1950) induce slower decomposition rate of OC in Belgium than in China, which is reflected by lower decomposition coefficients (k HUM , fr ecto and k RPM ) of resistant carbon pools.
Secondly, the distribution of temperature, precipitation and evaporation over the year also affects the decomposition rate and the loss of carbon from soil (Raich and Tufekcioglu, 2000;Schimel et al., 1994).High temperatures accompanied by significant rain occur in the summer monsoon climate of China, which could lead to quicker litter decomposition (Raich and Tufekcioglu, 2000;Zhang et al., 2008) without any limit of energy or moisture in this season.The different influences of climate conditions on decomposition process in two regions may be reflected indirectly by setting different  values of these parameters, because just decomposition coefficients (k HUM , fr ecto and k RPM ) of soil carbon pools were calibrated in this study, and not the mechanisms that mimic the effect of temperature and moisture on decomposition.Finally, because the effect of some soil properties (e.g., CaCO 3 content, pH) on OC cycle is not incorporated in Soil-Gen2 (except clay and water content), the difference of these properties in two regions may explain part of the variation in calibrated parameters.This relationship could be quantified via regression analysis based on many simulated and calibrated plots.However, the small number of calibrated plots in our study does not allow for such analysis.In addition, the different calibrated results in Belgium and China indicate that future calibrations are needed for distinct climate conditions, possibly for different soil types as well, and uncertainty bandwidths of calibrated parameters should be given.The current study however does not allow a certain statement on this topic as only a few soil type/climate combinations have been explored.

Conclusions
Sensitivity analysis based on the Morris method shows that k HUM , fr ecto and k RPM are the three most important parameters in SoilGen2 to affect change of OC both in Belgian and Chinese soil pedons.The sensitivity orders of the parameters follow the same pattern in the two regions, but the values of elementary effects differ.
According to the results of sensitivity analysis, SoilGen2 parameters are calibrated by decreasing k HUM , fr ecto and k RPM .The final results are obtained by the following combination of parameters: k HUM = 0.0065-0.0074,k RPM = 0.27-0.27and fr ecto = 0.30-0.38 in Belgium and k HUM = 0.016-0.019,k RPM = 0.30-0.30and fr ecto = 0.37-0.43 in China.The less significant sensitivity of parameters in the sensitivity analysis and larger variation of parameters during calibration in Belgium compared to China may be attributed to their distinct vegetation types and climate conditions.The calibrated parameters follow the law that OC in deeper soil layers is more resistant to decomposition than in surface soil layers, which is induced by the age of carbon and an unfavorable environment for soil micro-organisms in the deeper layers.This indicates that the calibration allows better simulation of carbon storage in the whole soil pedon.The calibrated SoilGen2 will be the base for quantitative soil carbon pool reconstruction based on the application of it to loess-soil sequences deposited in China in future studies, which will offer an opportunity to understand the mechanism of carbon cycle at geological timescale.

Fig. 1 .
Fig. 1.Structure and process parameters of the organic C-module of SoilGen2.indicates a distribution factor; is a rate factor.Process parameters are in italics; grey boxes indicate pools of C and associated ion species.The white square box is added for conceptualization, and white rounded boxes indicate processes.The dotted line indicates the model boundary.

Fig. 2 .
Fig. 2. Calibration procedure for OC cycle in SoilGen2 and the order of events.

Fig. 3 .
Fig. 3.Estimated means (µ * ) and standard deviations (σ ) of the distribution of elementary effects of factors on OC.(a) OC mass in ectorganic layers; (b) OC mass in endorganic layers; (c) OC content in endorganic layers.Closed symbols with names indicate the 4 most important factors.

Fig. 6 .
Fig. 6.Comparison of measured and simulated OC based on soil pedons for validation in Belgium and China.(a) OC mass in Belgium and China; (b) OC content of N pedon in Belgium; (c) OC content of P pedon in Belgium; (d) OC content of S pedon in Belgium; (e) OC content of LJB pedon in China; (f) OC content of ZW2 pedon in China; (g) OC content of ZW3 pedon in China.

Table 2 .
Model parameters in the organic C-module in SoilGen2 and results of sensitivity analysis.

Table 3 .
Model parameters in the organic C-module in SoilGen2 and quality indices of one group of pedons during the calibration in Belgium and China.Bold numbers indicate the best calibration results for each group of pedons.RPM frac ecto IM RMSE OCMS IM RMSE ecto IM RMSE1 endo RMSE2 endo DIS pedon