Information content of incubation experiments for inverse estimation of pools in the Rothamsted carbon model: a Bayesian perspective

. A major drawback of current soil organic carbon (SOC) models is that their conceptually deﬁned pools do not necessarily correspond to measurable SOC fractions in real practice. This not only impairs our ability to rigorously evaluate SOC models but also makes it difﬁcult to derive accurate initial states of the individual carbon pools. In this study, we tested the feasibility of inverse modelling for estimating pools in the Rothamsted carbon model (ROTHC) using mineralization rates observed during incubation experiments. This inverse approach may provide an alternative to existing SOC fractionation methods. To illustrate our approach, we used a time series of synthetically generated mineralization rates using the ROTHC model. We adopted a Bayesian approach using the recently developed DiffeRential Evolution Adaptive Metropolis (DREAM) algorithm to infer probability density functions of the various carbon pools at the start of incubation. The Kullback-Leibler divergence was used to quantify the information content of the mineralization rate data. Our results indicate that measured mineralization rates generally provided sufﬁcient information to reliably estimate all carbon pools in the ROTHC model. The incubation time necessary to appropriately constrain all pools was about 900 days. The use of prior information on microbial biomass carbon signiﬁcantly reduced the uncertainty of the initial carbon pools, decreasing the required incubation time to about 600 days. Simultaneous estimation of initial carbon pools and decomposition rate constants signiﬁcantly increased the of This effect was most for and Altogether, our results demonstrate that it is particularly to derive reasonable estimates of the pool and the inert organic matter from inverse modelling of mineralization rates observed during incubation experiments.


Introduction
A substantial part of the global carbon cycle occurs in the pedosphere. Jobbágy and Jackson (2000) estimated that the upper three meters of the soils of the world contain about 2340 Pg organic carbon. More recently, Tarnocai et al. (2009) reported that the amount of soil organic carbon (SOC) in the northern circumpolar permafrost region was previously underestimated. They concluded that the upper three meters of the soils in this region contain about 1020 Pg organic carbon. These numbers suggest that at the global scale total SOC amounts to about 3070 Pg in the upper three meters of the soil profile. This is about fourfold the amount of carbon in the atmosphere and sixfold that stored in living plant material (Prentice et al., 2001). The link between terrestrial carbon cycling and climate change has recently received much attention because the pedosphere has the potential to either amplify or dampen global warming (e.g. Friedlingstein et al., 2006;Heimann and Reichstein, 2008). Moreover, it has long been recognized that SOC levels exert strong influence on the physical and chemical properties of soils, and consequently, determine soil productivity and thus the functioning Published by Copernicus Publications on behalf of the European Geosciences Union. of terrestrial ecosystems. Both fields of science demand for reliable modelling of SOC dynamics.
The turnover of soil organic matter (SOM) is usually described with multi-compartment models. The model compartments (or pools) are generally defined by functional properties such as decomposition rate constants, input and output pathways, and partitioning of decomposition products. These models have been successfully applied for simulating total SOC levels observed in long-term field experiments (Smith et al., 1997). A major drawback of these models, however, is that their functionally defined compartments do not necessarily correspond to measurable SOC fractions in real practice (e.g. Christensen, 1996;Elliott et al., 1996;Cambardella, 1998;Arah and Gaunt, 2001). The ongoing difficulty with direct measurement of the various carbon pools not only impairs our ability to rigorously evaluate SOC models, but also downplays the utility of these models for predictive purposes. Accurate estimates of the initial states of the various carbon pools are a prerequisite for studying and predicting organic carbon dynamics in the soil.
In recent years, the Rothamsted carbon model (ROTHC, Coleman and Jenkinson, 1999) has found widespread application and use to study SOC dynamics. Several studies have developed fractionation procedures that yield SOC fractions that match as closely and consistently as possible the various pools considered in ROTHC (Jenkinson et al., 1992;Balesdent, 1996;Skjemstad et al., 2001;Ludwig et al., 2003;Skjemstad et al., 2004;Zimmermann et al., 2007). Notwithstanding this progress made, we think that this approach of "measuring the modelable" (Elliott et al., 1996) is subject to considerable uncertainty and debate. For instance, Smith et al. (2002) pointed out that to be equivalent to a model pool a measured fraction must be both unique (isolating all SOM pertaining to a conceptual pool) and non-composite (isolating SOM from one conceptual pool only). SOM fractionation is generally based on physical or chemical properties, or both. But it seems questionable if these properties alone or a combination thereof can be used to isolate unique and noncomposite SOC fractions that exhibit similar behaviour and fate as the individual carbon pools defined in ROTHC. In a recent review on SOM fractionation methods and their relevance for modelling, von Lützow et al. (2007) concluded that only the microbial biomass and the light fraction (SOM not firmly associated with soil minerals and consisting mostly of plant residues) can be reliably isolated. They argue that the currently available fractionation techniques are not specific enough with regard to SOM stabilization mechanisms, and therefore, do not yield functional pools with homogeneous turnover rates. This might explain why the opposite strategy of "modelling the measurable" (Elliott et al., 1996) -that is, building mathematical models based on measurable SOM fractions -has not yet found its way into modelling practice.
In this study, we explored the feasibility of inverse modelling to estimate the initial carbon pools in the ROTHC model. This approach may provide an alternative to exist-ing SOM fractionation methods. To test our approach, we used a synthetic time series of mineralization rate data. This has the advantage that the true pool sizes are known, which enables us to benchmark the reliability and robustness of our inverse modelling approach. In our analyses we considered measurements made during incubation because such experiments are admirably suited to maximize the information content of the measured mineralization rate data by controlling the boundary conditions. Incubation experiments also reveal the dynamics of the intermediate and slow pools after depletion of the labile compounds. To appropriately infer uncertainty of the initial states of the carbon pools, we posed our inverse problem in a Bayesian framework using the recently developed DiffeRential Evolution Adaptive Metropolis (DREAM) algorithm (Vrugt et al., 2008(Vrugt et al., , 2009. We explicitly addressed the following questions: (i) do incubation experiments provide sufficient information to reliably estimate all carbon pools in the ROTHC model, especially those pools with intermediate and slow turnover? (ii) What length of incubation is required to appropriately constrain all carbon pools? (iii) Would prior information on microbial biomass carbon reduce uncertainty of the initial pool size estimates? And (iv), how does uncertainty in some of the ROTHC model parameters affect the final carbon pool size estimates?

Model description
A schematic overview of the Rothamsted carbon model (ROTHC, Coleman and Jenkinson, 1999) appears in Fig. 1. The ROTHC model has four active carbon pools: decomposable plant material (DPM), resistant plant material (RPM), microbial biomass (BIO) and humified organic matter (HUM). The mass balance of the ith active pool, c i (g C kg −1 soil), is described with the following ordinary differential equation: where t (year) is time, k * i (year −1 ) denotes the actual decomposition rate, and s i (g C kg −1 soil year −1 ) represents inputs. The actual decomposition rate k * i is modelled as a multiplicative product of the decomposition rate constant, k i (year −1 ), and three rate modifying functions, f (dimensionless), accounting for the effects of soil temperature, soil moisture, and soil cover: Decomposition rate constants (Table 1) were derived from field experiments with 14 C labelled plant material and 15 N labelled microbial biomass, as well as measurements of total SOC from long-term field experiments at Rothamsted (Jenkinson et al., 1992). These constants are assumed to be time  (Coleman et al. , 1997). The incoming plant material is partitioned between DPM and RPM in a ratio that depends on vegetation type, p 1 (vegetation). The decomposition process releases part of the carbon as CO 2 to the atmosphere and distributes the remaining carbon between the BIO and HUM pools. The proportion of carbon that is mineralized, p 2 (clay), depends on the clay content of the soil whereas a fixed ratio, p 3 , is used to partition the carbon flux between the BIO and HUM pools. A detailed description of the rate modifying functions and partitioning coefficients appears in Coleman and Jenkinson (1999) and so will not be repeated here. In addition to the four active pools (DPM, RPM, BIO and HUM), the ROTHC model also contains a fifth pool that is made up of inert organic matter (IOM). This pool is resistant to microbial decomposition and does not receive any carbon from the four active pools or other sources. It is assumed that the carbon in the IOM pool is of geological rather than pedological age, implying that it virtually contains no 14 C. The IOM pool is necessary to be able to match the 14 C signature of soils, which usually indicates the presence of some small amount of very old carbon (Jenkinson et al., 1992). In all our calculations, we used a MATLAB implementation of the ROTHC model with an explicit Runge-Kutta numerical solution scheme to solve for the system of ordinary differential equations. This scheme is especially designed to maintain mass balance, a requirement for efficient inverse modelling (Schoups et al., 2010). Precision of mineralization rate measurements using an automated respirometer. The circles represent the actual measurements, error bars depict the width of the corresponding 95% confidence intervals. The solid line represents the fit of the linear error model, dashed lines mark the associated 95% confidence envelope.

Generation of synthetic incubation data
We generated a time series of synthetic mineralization rates using the ROTHC model with values of the initial pool sizes and model parameters listed in Table 1. The sizes of the initial carbon pools were obtained from a simulation of the Broadbalk Continuous Wheat experiment following the procedure and input data given in Jenkinson et al. (1992). In all our simulations, incubation temperature was set to 20 • C, soil moisture was assumed to be at an "optimal" level and soil cover was set to "bare". Following Ellert and Bettany (1988) and Hess and Schmidt (1995), we used rates of mineralization rather than cumulative amounts during inverse modelling to obtain realistic estimates of pool and parameter uncertainty.
We simulated an incubation experiment of 1200 days. From this data set, we selected n = 80 different measurements with equidistant intervals on a transformed time scale t 0.6 . This approach ensures that more measurements are available at the start of incubation when mineralization rates typically show the steepest decline. To include the effect of measurement error, simulated mineralization rates were corrupted with a normally distributed (Gaussian) error. The size of this random error was determined using a data set of actual mineralization rate measurements conducted with an automated respirometer. Figure 2 plots the standard deviation of the observed mineralization rates against their mean value. We used a total of 8 different samples with 12 replicate measurements on each sample. This was the maximum number of repetitions possible within a 24 h interval using our respirometer setup. The scatter plot clearly illustrates the presence of heteroscedasticity: the measurement error increases with increasing mineralization rate. To use this information in our Bayesian analysis, we fitted a straight line through the data with intercept of 0.00014 and slope of 0.019 g C kg −1 soil day −1 . This error model was used to calculate the measurement error for each observation of mineralization rate. We also used 12 realizations to generate the synthetic data. From these realizations, the mean and standard deviation of each observation was calculated. The synthetic data set is depicted in Fig. 3a. The error bars, which represent 95% confidence intervals of the synthetic observations, are hardly visible and mostly hidden behind the circles that represent the mean values. The residuals between (true) model simulation and generated measurements are shown in Fig. 3b. Both the tight confidence intervals and the small residuals illustrate the high precision of the synthetically generated mineralization rate data. This synthetic data set was used in our subsequent analyses. To investigate the information content of incubation experiments for estimating the initial carbon pools, we varied the duration of incubation from 25 to 1200 days. In the remainder of this paper, we focus on incubation experiments with length of 300, 600, and 900 days with n=34, 52, and 67 observations of mineralization rate, respectively. Consistent conclusions were drawn for the other durations.

Inverse modelling
To help describe the inverse method utilized herein, let us denote the ROTHC model with the symbol M that simulates n mineralization rates R={R 1 ,...,R n } from the initial states S, with measured forcing variablesF , and model parameters P : where S is the subject of inference in the rest of this paper, as it contains the sizes of the individual carbon pools at the start of the ROTHC simulation. The forcing variablesF include observed soil temperature and soil moisture, whereas P consists of decomposition rate constants (k), rate modifying factors (f ) and partitioning coefficients (p). Inverse modelling is a general method for estimating some of the arguments of M by minimizing a predefined measure of the residual vector ε={ε 1 ,...,ε n }. The residuals are defined by the differences between simulated variables R and corresponding measurementsR={R 1 ,...,R n }: The closer the residuals are to zero, the better the model represents the observational data. However, because of errors in observed forcing variablesF , structural inadequacies in the model (the model is not perfect), errors in the actual response variable measurementsR, and uncertainty in the model parameters P , the residuals are not expected to go to zero. Note that our use of inverse modelling in this paper is rather different than typical implementations of this methodology. Here, we estimated the initial states whereas most inverse modelling approaches focus on estimating the model parameters.
A common measure that is minimized during inverse modelling is the weighted least squares objective function: where σ 2 i is the variance of the ith measurement. For models that are nonlinear in their estimated arguments the minimum of (S|F ,P ,R) cannot be found analytically. Various optimization methods have therefore been developed during the past decades to efficiently minimize Eq. (5) for multidimensional search spaces. Unfortunately, most of these algorithms only provide an estimate of the best values of the estimated arguments. This would suffice in many practical applications, yet it would also be desirable to have an estimate of the uncertainty of the final optimized arguments. A probability density function (pdf) will help assess the information content of the measured variables and help generate predictive distributions of R.
Bayesian statistics provide a generic framework to estimate uncertainty of state variables, parameters and model predictions. In this approach, the estimated arguments are treated as probabilistic variables having a joint posterior pdf, p(S|F ,P ,R). The posterior pdf summarizes what is known about the estimands given the measurements and prior information. If we assume that the measurement errors are independent and normally distributed with individual variance σ 2 i , the posterior pdf takes the following form: where p(S) signifies the prior pdf of S. It summarizes the information on S before any measurements are available. In many applications of Bayesian statistics prior knowledge about the estimands is typically vague. In that case a noninformative (uniform) prior distribution is usually imposed. The Bayesian framework of statistical inference and prediction offers several important advantages over the classical frequentist approach. A detailed description and overview of both approaches can be found, for example, in Reichert and Omlin (1997) and Omlin and Reichert (1999). Unfortunately, for most practical problems the posterior pdf in Eq. (6) cannot be obtained by analytical means nor by analytical approximation. We therefore resort to iterative approximation methods such as Markov chain Monte Carlo (MCMC) simulation to generate a sample from the posterior pdf. The basis of the MCMC method is a Markov chain that generates a random walk through the search space with stable frequency stemming from a fixed probability distribution, p(S|F ,P ,R). Here, we use the DiffeRential Evolution Adaptive Metropolis (DREAM) algorithm that was recently introduced by Vrugt et al. (2008Vrugt et al. ( , 2009). The DREAM sampling scheme is an efficient MCMC sampler that runs multiple chains simultaneously for global exploration of the search space and automatically adapts the scale and orientation of the proposal distribution during the evolution to the posterior distribution. This scheme is an adaptation of the Shuffled Complex Evolution Metropolis algorithm (Vrugt et al., 2003) and has the advantage of maintaining detailed balance and ergodicity while showing good efficiencies on complex, highly nonlinear and multi-modal target distributions (Vrugt et al., 2009).
In all our calculations reported herein, we assumed that soil temperature and soil moisture are constant and known during the entire incubation experiment. In the first case, we fixed all model parameters at their standard values (Table 1) and inversely estimated the sizes of the active pools (DPM, RPM, BIO and HUM) at the beginning of the incubation experiment. We further assumed that no prior information was available for any of these (Sect. 3.1). The inert pool (IOM) is by definition not decomposed, and hence, observed mineralization rates do not contain any information about the size of this particular carbon pool. Consequently, its size cannot be estimated from measured mineralization rates. The size of the IOM pool was therefore calculated from the difference between total SOC measured at the beginning of the experiment and the combined size of the other carbon pools. The amount of total initial SOC (28.6 g C kg −1 soil) was assumed to be known. Our second case study investigated whether prior information of microbial biomass carbon at the beginning of incubation will reduce the uncertainty associated with the various carbon pools. We used this information to formulate an informative prior for the BIO pool (Sect. 3.2). Finally, our third case study jointly estimated decomposition rate constants and initial pool sizes to address ROTHC model parameter uncertainty explicitly (Sect. 3.3). Our first two studies assume that the ROTHC model parameters are perfectly known, whereas in reality these are subject to considerable uncertainty. The findings from this last study were used to make recommendations about the usefulness of inverse modelling for determining the initial carbon pool sizes in the ROTHC model.

Kullback-Leibler divergence
We used the Kullback-Leibler divergence (Kullback and Leibler , 1951) to quantify the information content of the synthetic time series of mineralization rates. Kullback-Leibler divergence, D KL (dimensionless), is a non-symmetric measure of the dissimilarity between two probability distributions. In the Bayesian context of this study, it was used to quantify the gain of information that results from moving from the marginal prior distribution of the ith state variable, p(S i ), to its posterior, p(S i |F ,P ,R). For continuous variables Kullback-Leibler divergence can be defined as: = p(S i |F ,P ,R)ln where the integral is evaluated over the feasible space of the state variable S i . We solved Eq. (7) using numerical integration.

Informative prior for microbial biomass carbon
In this study, we investigated whether an informative prior for c BIO would improve the identifiability of this and possibly other carbon pools. For that purpose, we assumed that microbial biomass carbon was measured at the beginning of the experiment using the fumigation-extraction method (Vance et al., 1987). Jenkinson et al. (1992) demonstrated that microbial biomass carbon measured using fumigationextraction is closely related to simulated c BIO in the ROTHC model. The basic principle of the fumigation-extraction method is that soil microorganisms die after chloroform fumigation and part of the microbial constituents is degraded by enzymatic autolysis and transformed into extractable components (Joergensen, 1996). Microbial biomass carbon, C BIO (g C kg −1 soil), is calculated as: where E C (g C kg −1 soil) is the organic carbon extracted from fumigated soil minus that extracted from nonfumigated soil and K EC (dimensionless) is the extractable part of microbial biomass carbon after fumigation. Throughout this paper we used the upper case C to denote measurable fractions and lower case c to denote the modelled pools. The K EC value can be obtained from direct or indirect calibration but most often a fixed value of 0.45 (Wu et al., 1990) is used. As Joergensen (1996) pointed out, the measurement of E C is very precise (having a coefficient of variation <5%), but the estimate of C BIO is less certain because of variation of K EC between different soils.
To derive a probability density function of C BIO that can be used as an informative prior for c BIO , we used the K EC values presented in Joergensen (1996) for arable soils. A normal distribution with mean 0.42 and standard deviation 0.08 was shown to describe this data set well. This probability distribution was then used in conjunction with Eq. (8) to obtain the pdf of C BIO . It is noteworthy that this pdf exhibited significant skew to the left. That is, even if the measurement of E C and the mean of K EC were unbiased, the mode of the probability distribution of C BIO did not exactly match its true value. We did not impose a measurement error on E C because its magnitude was considered negligible compared to the uncertainty that stems from the natural variability of K EC .

Non-informative priors
Figure 4 displays marginal probability density distributions of the initial pool sizes using measured mineralization rates from 300, 600 and 900 days of incubation. To quantify the information content of the individual experiments, the Kullback-Leibler divergence is reported for each experiment and active carbon pool. Incubation experiments of 300 days provided insufficient information to warrant the inverse estimation of all carbon pools in the ROTHC model. Only c DPM and c RPM were well defined by calibration against measured mineralization data. The other carbon pools could not be constrained with posterior distributions extending over the entire prior defined ranges. Incubation experiments of 600 and 900 days improved the identifiability of the c BIO , c HUM and c IOM pools. The posterior pdfs of these pools became increasingly well identified with increasing length of the incubation experiment. Unfortunately, even after processing all mineralization data some residual uncertainty remained in the estimates of the initial carbon pools. Also, the estimates of the carbon pools deviated somewhat from their actual values used to generate the synthetic mineralization data. This is because of measurement error and the limited size of the data set. The posterior pdf of the IOM pool closely followed that of the HUM pool. This is not surprising since HUM makes up the largest part of total SOC. Note that the IOM pool was not estimated from the mineralization data but calculated from the difference between total SOC and carbon stored in the four active pools. Any uncertainty in the estimate of c HUM is therefore directly propagated into the estimate of c IOM . This finding highlights the need for an accurate estimate of the HUM pool.
To better understand the previous results consider Table 2 that reports the Pearson correlation coefficients of the posterior samples using 900 days of incubation. In general, strong correlation between individual carbon pools impairs their identifiability because a change in one pool can be compensated for by a change in another pool. A perfect (linear) correlation was found between c HUM and c IOM , which reflects how IOM is calculated. Most active carbon pools exhibited considerable correlation with Pearson coefficients up to −0.94. The only carbon pool that was relatively uncorrelated with the others was c DPM . This explains why this particular pool was best identified in all our synthetic experiments. The correlation structure induced in the posterior samples is determined by the pathways of carbon flow within ROTHC (Fig. 1) and the kinetic properties of the various pools (Table 1).
In Fig. 5a we plotted the evolution of the Kullback-Leibler divergence with increasing length of incubation. Obviously, there was only little improvement in the estimate of c DPM after about 200 days of incubation. At this time, DPM was almost completely mineralized (with less than 1% of the initial amount remaining). On the contrary, the first 200 days contained very little information about c HUM . Extending the incubation experiment beyond 900 days only marginally reduced the uncertainty of the various carbon pools. This raises the question when to stop an incubation experiment because the gain from more experimental data would be marginal. Whereas traditional inverse modeling endeavors start when data collection has terminated, arguably there is significant advantage of simultaneous data collection and inverse modeling in this case. Such an approach continuously updates the posterior distributions of the various carbon pools when new mineralization data become available. A description of such sequential MCMC methods can be found, for example, in Del Moral et al. (2006). The Kullback-Leibler divergence can be especially helpful in this context to judge when to stop the incubation experiment.

Informative prior for the BIO pool
The results in the previous section have shown that relatively long incubation times were required to appropriately constrain all carbon pools in the ROTHC model. Such long experiments are not only impractical, but also downplay the utility of inverse modelling for assessing SOC dynamics. This limitation inspires thinking into alternative measurements to improve the identifiability of the various carbon pools. Let us assume that microbial biomass carbon was measured separately at the beginning of the incubation experiment using the fumigation-extraction method. Let us also assume that this measurement was relatively unbiased. Figure 6 shows probability distributions of the individual carbon pools when explicitly using prior information on c BIO in our analysis with DREAM. Each of the three horizontal panels again consider a different length of incubation. Compared to the previous results with a non-informative prior for c BIO (Fig. 4), the identifiability of the various carbon pools had clearly improved. The various pools were better defined by calibration against measured mineralization rates with posterior pdfs that are considerably tighter around their true values used to generate the synthetic data. This is not surprising. Explicit prior information about microbial carbon reduces the acceptable range of c BIO , which in turn reduces the uncertainty of the other carbon pools due to their correlation structure (Table 2). It is interesting to observe that when an informative prior for c BIO was used, measured mineral- ization data could not further reduce the uncertainty of this pool. The posterior distribution was virtually identical to the prior distribution for all lengths of incubation. Figure 5b shows the corresponding evolution of the Kullback-Leibler divergence criteria with increasing length of the incubation experiment. This diagnostic measure showed that c RPM benefited the most from the use of an informative prior for c BIO . This improvement was most pronounced at shorter incubation times. Altogether, these results demonstrate that an informative prior for the BIO pool improved the identifiability of the various carbon pools, thereby reducing the required length of incubation to about 600 days.
In many practical applications, however the measurement of C BIO might be biased due to an erroneous estimate of K CE , the extractable part of microbial biomass carbon after fumigation. Here we examined the effect of a biased prior of the BIO pool on the inference of the ROTHC carbon pools. As illustration, we assume that the true value of K CE was overestimated by one standard deviation, which results in an underestimation of C BIO (see Eq. 8). Figure 7 shows the resulting marginal posterior probability distributions. Indeed, the biased prior had a dominating effect on the posterior distribution of the various pools. Its influence slightly decreased with increasing length of incubation. This is the effect of the increasing number of mineralization data points, and hence, information becoming available. In general, underestimation of c BIO resulted in an overestimation of c RPM and an underestimation of c HUM . This is consistent with the correlation structure between these various carbon pools as found in the posterior samples. Notwithstanding, the estimates are still reasonable, with their true values lying within the 95% confidence intervals of the respective posterior distributions. Similar results were obtained when the prior of the BIO pool overestimated the actual value of c BIO (not shown).
It is important to realize that the use of an informative prior of the BIO pool does not simply fix the value of initial c BIO to measured microbial biomass carbon. The use of an informative prior is advantageous for two reasons. First, the informative prior appropriately treats the uncertainty in measured C BIO , which is then propagated into the posterior pdfs of the remaining carbon pools. This preserves the correlation structure of the carbon pools and should therefore result in more realistic estimates of uncertainty. Second, if the C BIO measurement is biased, the additional information contained in the time series of mineralization data may partially correct for this bias. In the present study, however, this mitigating effect was relatively weak. To ensure that the measurement of C BIO is unbiased, we therefore recommend using direct calibration methods such as in situ labelling of microorganisms (Bremer and van Kessel, 1990;Dictor et al., 1998). This is important because the bias in the prior of c BIO was not fully compensated for by the information contained in the time series of observed mineralization rates. To further improve our results, replicate in situ labelling measurements can be used to calculate the standard deviation of the sample distribution of K CE (Dictor et al., 1998). This information can be used to formulate an informative prior of c BIO that is -based on the data given in Dictor et al. (1998) -more accurate than the prior we used in the present study. Future work should also focus on evaluating ROTHC against a time series of observed microbial biomass carbon. This will help establish whether ROTHC can describe the decline in microbial biomass that is commonly observed during long-term incubation experiments (e.g. Nicolardot et al., 1994;Follett et al., 2007).
One possible way forward to reduce uncertainty in the initial states of the various carbon pools would be to make use of explicit prior information about the SOC light fraction. This respective fraction (also termed particulate organic carbon) is conceptually equivalent to the sum of c DPM and c RPM in the ROTHC model. Zimmermann et al. (2007) used density fractionation to isolate particulate organic carbon (<1.8 g cm 3 ) and found good agreement between their measurements and modelled c DPM + c RPM . Shirato and Yokozawa (2006) proposed to use acid hydrolysis to partition plant material into the decomposable and resistant pools considered in ROTHC. The combination of both methods seems particularly promising to derive informative priors for these two different pools. Our previous results suggest that prior information on c RPM could greatly improve the estimate of c HUM (and consequently c IOM ). This would further reduce the required length of incubation. Moreover, we think that the Bayesian approach utilized herein could prove useful to independently test the proposed fractionation methods in a reliable and comparatively quick manner.

Informative prior of the BIO pool and model parameter uncertainty
In the previous two sections, we focused on estimating the uncertainty in the initial carbon pool sizes of the ROTHC model using measurements of mineralization rate data collected during incubation experiments. This approach, however, neglects uncertainty in the model parameters. In this section, we explore how uncertainty in the model parameters affects the final posterior distributions of the various carbon pools. To this end, we considered the decomposition rate constants to be unknown and added these parameters to our inverse estimation with DREAM. We used informative priors for all decomposition rate constants, k, and jointly estimated these model parameters and the various carbon pools, c. The priors of the decomposition rate constants followed a normal distribution with mean corresponding to the standard values of k used in ROTHC and standard deviation equal to a coefficient of variation of 25%. These informative priors thus reflected our belief that the true values of the decomposition rate constants remain in close vicinity of their standard values used in ROTHC. The bounds of the feasible parameter space were set to ± three times the standard deviation of the prior pdfs (Table 1). Figure 8 shows the corresponding marginal posterior pdfs of initial pool sizes, again using an (unbiased) informative prior for the BIO pool. Obviously, the uncertainty in the k estimates has propagated into that of c. The estimates of initial pool sizes now showed substantial more uncertainty compared to the case where model parameter uncertainty was ignored (Fig. 6). This is especially true for c HUM (and consequently c IOM ). Even after 900 days of incubation the identifiability of these two pools was rather poor. Also, the uncertainty of the RPM pool had clearly increased, now requiring at least 600 days of incubation to reasonably constrain this pool. The posterior of c BIO closely followed its prior distribution, demonstrating that measured mineralization rates did not contain any additional information about this pool. The DPM pool required only 300 days of incubation to be well defined. The evolution of the Kullback-Leibler divergence for the different initial pool sizes is depicted in Fig. 5c.
The marginal posterior pdfs of decomposition rate constants are shown in Fig. 9. The posterior pdfs of k BIO and k HUM closely followed their prior distributions, indicating that observed mineralization rates contained very little additional information to further constrain the rate constants of these two pools. The posterior distribution of k RPM was better constrained but time evolution of its pdf showed inconsistent behaviour, underestimating the true decomposition rate after 300 days of incubation. Most likely, this was due to the particular realization of the measurement error in this early phase of incubation. The pdf of k DPM was the only marginal distribution that showed a distinct reduction in uncertainty after processing the mineralization rate data with DREAM. This was true even when using only the first 300 days of incubation.
Note that our approach of estimating decomposition rate constants, k, is equivalent to estimating actual decomposition rates, k * . This is because of the multiplicative formulation used in Eq. (2). This approach implicitly accounts for the uncertainty in the rate modifiers. We could also explicitly account for this model parameter uncertainty, but this requires that soil samples are exposed to varying temperature and moisture regimes during incubation. This would further  increase the uncertainty of the estimates of the initial carbon pool sizes. Treating the partitioning coefficients, p, as additional variables would have a similar effect.
Another point that deserves special consideration is model structural uncertainty. Inherent to the inverse approach is the assumption that the model is a valid description of reality. This was certainly true for our synthetic case study because we used the same model to generate the data and to infer the pool sizes. We should have in mind, however, that "environmental models are inherently imperfect because they abstract and simplify real patterns and processes that are themselves imperfectly known and understood" (Brown and Heuvelink, 2005). These imperfections will likely be revealed in systematic deviations between observed and simulated mineralization rates. In real practice, we may therefore expect that model structural uncertainty will play a significant role and further increase the uncertainty in the carbon pool size estimates. Our results demonstrate that -even with a perfect model -the estimation of c HUM (and consequently c IOM ) is particularly difficult, when model parameter uncertainty is explicitly accounted for.
Our results comply with related studies in which MCMC methods were used to constrain parameters and initial conditions in terrestrial ecosystem models. For instance, Xu et al. (2006) determined posterior distributions of transfer coefficients from observations of biomass growth, litterfall, carbon content of the litter layer and mineral soil, and soil respiration. Only transfer coefficients linked to fast processes were warranted by calibration against the experimental data. Contrary to our study, the initial state of the system was assumed to be known. Fox et al. (2009) used real and synthetically generated time series of net ecosystem exchange of CO 2 and leaf area index to simultaneously estimate model parameters and initial states, but without recourse to estimating the SOM pool. Also in this work, the parameters describing fast processes were best identified. Moreover, the turnover rate of SOM was well constrained by available data. Fox et al. (2009) attributed this, among others, to the known initial state of the SOM pool. Both studies clearly illustrated an inability of current measurement technologies and data sources to appropriately constrain all important model parameters and states of the system under study.

Summary and conclusions
Results presented in this paper illustrate that, in principle, all carbon pools considered in the ROTHC model can be estimated from mineralization rates observed during incubation experiments. Yet, about 900 days of incubation were required to appropriately constrain all pools, especially those with intermediate and slow turnover. Such long experiments are required because of significant correlation between the various carbon pools. The use of prior information on microbial biomass substantially reduced the uncertainty in the pool size estimates. This in turn reduced the required incubation time to about 600 days. Care must be taken, however, that the measurement of microbial biomass carbon is unbiased because measured mineralization rates did not contain sufficient information to fully compensate for this error. This led to bias in the estimates of all pool sizes. When part of model parameter uncertainty was explicitly accounted for, the uncertainty in all initial carbon pool estimates increased considerably. This effect was most pronounced for the intermediate and slow pools. Altogether, our results demonstrate that measured mineralization rates do not warrant the identification of the HUM and IOM pools, especially in the face of model parameter uncertainty. Inverse estimation of the initial states of these two carbon pools might therefore not be feasible in real practice. Other data sources or measurement techniques are required to appropriately constrain the various carbon pools of the ROTHC model.