Evaluation of a soil greenhouse gas emission model based on Bayesian inference and MCMC: Model uncertainty

Ecological Modelling 253 (2013) 97–106

Contents lists available at SciVerse ScienceDirect

Ecological Modelling journal homepage: www.elsevier.com/locate/ecolmodel

Evaluation of a soil greenhouse gas emission model based on Bayesian inference and MCMC: Model uncertainty Gangsheng Wang a,b,∗ , Shulin Chen a,∗ a b

Department of Biological Systems Engineering, Washington State University, Pullman, WA 99164-6120, USA Environmental Sciences Division, Oak Ridge National Laboratory, Oak Ridge, TN 37831-6301, USA

a r t i c l e

i n f o

Article history: Available online 6 October 2012 Key words: Bayesian inference Greenhouse gas (GHG) Markov Chain Monte Carlo (MCMC) Metropolis–Hastings algorithm Model uncertainty

a b s t r a c t We combined the Bayesian inference and the Markov Chain Monte Carlo (MCMC) technique to quantify uncertainties in the process-based soil greenhouse gas (GHG) emission models. The Metropolis–Hastings sampling was examined by comparing four univariate proposal distributions (UPDs: symmetric/asymmetric uniform and symmetric/asymmetric normal) and one multinormal proposal distribution (MPD). Almost all the posterior parameter ranges from the MPD could be reduced to 1 order of magnitude. The simulation errors in CO2 fluxes were much greater than those in N2 O fluxes, which resulted in a greater importance in model structure than in model parameters for CO2 simulations. We suggested deriving the covariance matrix of parameters for MPD from the sampling results of a UPD; and generating a Markov chain by updating a single parameter rather than updating all parameters at each time. The method addressed in this paper can be used to evaluate uncertainties in other GHG emission models. © 2012 Elsevier B.V. All rights reserved.

1. Introduction Process-based mathematical models have been developed to simulate the greenhouse gas (GHG) emissions as an important part of the carbon–nitrogen dynamics in soils (Chen et al., 2008; Ma and Shaffer, 2001; Smith et al., 1997; Wu and McGechan, 1998). However, studies on the uncertainties in these models and model applications are limited (Wang and Chen, 2012). A subjective interpretation of uncertainty is “the degree of confidence that a decision maker has about possible outcomes and/or probabilities of these outcomes” (Refsgaard et al., 2007). Uncertainty assessment is therefore important when models are used for decision-making (Yohe and Oppenheimer, 2011). Uncertainty analysis not only gives the uncertainty from different sources (i.e., model parameters, model structure, model inputs and outputs), but also gives an evaluation of model performance and limitations. Some studies on GHG models are concentrated on sensitivity analysis, which is one of the methods described in Refsgaard et al. (2007). The uncertainty of the PnET-N-DNDC model was evaluated by examining the sensitivity of the model outputs to such environmental factors as temperature, precipitation, photosynthetically active radiation (PAR) and model input variables, e.g., N-concentration in precipitation, litter mass, soil organic carbon

∗ Corresponding authors at: Department of Biological Systems Engineering, Washington State University, Pullman, WA 99164-6120, USA. E-mail addresses: [email protected] (G. Wang), [email protected] (S. Chen). 0304-3800/$ – see front matter © 2012 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.ecolmodel.2012.09.010

(SOC), pH, and soil texture (Stange et al., 2000). The sensitivity analysis was conducted by changing one factor at a time while keeping all others constant. Monte Carlo is the most widely used method in uncertainty anaylsis. Thorsen et al. (2001) used Monte Carlo to analyze the propagation of uncertainty from input data to model output, and found that the magnitude of uncertainty was closely associated with the investigated spatial scale, i.e., smaller output uncertainty on catchment scale than on grid level. The Monte Carlo technique was also applied to assess the uncertainty of DenNit model output with regard to parameterization (Reth et al., 2005). A Gaussian distribution within one standard error of mean was used for parameter sampling. In the comparison of carbon and nitrogen dynamics under conditions of conventional and diversified rotations, 64 parameter combinations were identified to test their impacts on model outcomes in Tonitto et al. (2007). An uncertainty analysis tool for the DNDC model allows for the selection of either the Monte Carlo or the Most Sensitive Factor (MSF) method (Li et al., 1992a). The MSF method involves running the model twice for each spatial unit (e.g., grid cell or polygon) with the maximum and minimum parameter values. These two runs generate two gas fluxes to form an interval, which is assumed to cover the real gas flux with a high probability. Using this uncertainty model, several highly sensitive factors influencing DNDC model were identified (Li et al., 1992a,b). The Monte Carlo analysis was also adapted to assess uncertainties in soil N2 O simulations from model input and structure of DAYCENT (Del Grosso et al., 2010). Probability distribution functions (PDFs) of major model inputs (weather, soil texture, and N applications)

98

G. Wang, S. Chen / Ecological Modelling 253 (2013) 97–106

were assigned to quantify the uncertainty due to model inputs. The structural uncertainty was a synthesis of parametric uncertainty and model residual errors derived by an empirically based linear mixed effect model (Ogle et al., 2007). The aforementioned methods usually did not take into account the probability distribution of parameters. Even if a distribution was considered, e.g., a Gaussian distribution in Reth et al. (2005), it was a priori. The posterior distribution of a parameter is more informative than a priori for modelers to use and evaluate a data-driven model (Hartig et al., 2011). The Bayesian inference based on Markov Chain Monte Carlo (MCMC) has been used to calibrate parameters and quantify parametric uncertainty in the N2 O submodel of CERES-EGC (Lehuger et al., 2009), where the uniform probability distributions were assigned as a priori for 11 global parameters. The model structure in Del Grosso et al. (2010) referred to both model parameters and structure (quantified by simulation errors); whereas model parameters and structure were usually regarded as two distinct sources for total uncertainty (Refsgaard, 1997; Wang et al., 2009). Lehuger et al. (2009) focused on parametric uncertainty and used the observation errors (standard deviations) of N2 O fluxes to represent the simulation errors in the likelihood estimator. However, the simulation errors may be treated as a latent variable and incorporated into the Bayesian framework with model parameters for quantifying the model structural uncertainty (Wang and Chen, 2012). The objective of this paper was to evaluate model uncertainties due to model parameters and structure by coupling the Bayesian theory with MCMC method. Based on Bayes’ theorem, the posterior distribution of both model parameters and model output variance can be derived from the prior distribution and observed outputs, and the 95% confidence intervals (CIs) of any output variables due to parameter uncertainty and model structure uncertainty can be estimated. The uncertainties from these two sources were also compared with that due to observation errors in the GHG fluxes. A soil GHG emission model (Appendix A) was used to test the proposed uncertainty analysis method. 2. Material and methods 2.1. Bayesian inference According to Bayesian inference (Hartig et al., 2011), the posterior distribution (|yt ), i.e., the likelihood function L(|yt ), of parameter set  can be generated from the prior distribution f() conditioned on observations yt :



  |yt





=



   



f yt | · f 

 



f yt | · f  · d



 

∝ f yt | · f 

(1)

where (,  y ) is a vector including the model parameter set () and the standard deviation ( y ) depicting simulation errors; f(yt |) is the distribution function of model output variable yt conditioned on ; and t is a time index. Generally, yt is a transformation of the model output Yt to obtain a homoscedastic variance for the simulation errors (Kuehl, 1999).The square-root transformation Yt , is adopted in this study. (Engeland et al., 2005), i.e., yt = As for the prior distribution, a common approach is to assume uniform priors (Iskrev, 2007), which means f() is a constant. In addition, the model output (yt ) is assumed to follow a normal distribution (Congdon, 2001):





f yt | = 1/



Thus

  |yt



√



 

2 

2y · exp − yt − yt ()





 

/ 2y2



2 

∝ f yt | ∝ 1/y · exp − yt − yt ()

/ 2y2

(2)



(3)

 

 ∗ |yt  k |yt

   2  2  yt − yt ( ∗ ) yt − yt ( k ) yk  = ∗ exp − − (4)  2  2 y

2 y∗

2 yk

When all the simulation errors are assumed to be independent, the likelihood of the model outcome can be expressed as the product of the likelihood of each individual outcome at each time step



  |y



=



  |yt



(5)

t

2.2. Metropolis–Hastings algorithm for MCMC The Metropolis–Hastings (MH) algorithm is a typical Markov Chain Monte Carlo (MCMC) sampling method to randomly sample from the posterior distribution described by Eq. (5). The procedure of MH can be found in many reports (Chumbley et al., 2007; Hastings, 1970; Link and Barker, 2008; Mathe and Novak, 2007; Tiana et al., 2007). An important criterion in MH is the acceptance probability:



∗

a x |x

k





 = min 1,   

 (x∗ ) J xk |x∗

 

 xk J x∗ |xk

(6)

where a(x* |xk ) denotes the acceptance probability; xk is the current state of the chain; x* is the new state of the chain generated from xk using a specified irreducible proposal distribution J(x* |xk ); and (• ) is the posterior distribution function defined by Eq. (5). If we draw a random number (Z) from the uniform distribution U(0,1), then a new state (i.e., k + 1) of the chain can be determined by



xk+1 =









x∗ ,

if a x∗ |xk ≥ Z

xk ,

if a x∗ |xk < Z

(7)

Four univariate (symmetric and asymmetric uniform, symmetric and asymmetric normal) and one multivariate (symmetric multinormal) proposal distributions were examined (see Appendix B). In the transition from the current state to a new state, three strategies may be used (Hastings, 1970). (i) All elements in the parameter vector are changed. In this case, x* and xk are no longer the one-dimensional parameters as in the previous four distributions, they are vectors containing d parameters. (ii) Only one of the elements is randomly selected and changed. (iii) Only one element is changed, and this element is selected in a fixed, rather than a random, sequence. It is often inefficient to sample small values in MCMC if a parameter ranges several orders of magnitude. Thus we conducted MCMC pertaining to the logarithmic transformation of parameter values. Hereafter the log-transformed parameter space is called the logarithmic parameter space. 2.3. Generating random samples from multivariate normal distribution It is not as simple to implement a multivariate normal distribution as a univariate distribution. It is required to know the covariance matrix in advance, and to randomly generate a vector, not a single parameter value, from the distribution. Although the true covariance matrix is unknown, it can be approximately estimated from the parameter samples generated by MCMC using any of the four univariate proposal distributions (UPDs). The algorithm described in Hernadvolgyi (1998) was followed to generate random vectors from a multinormal proposal distribution (MPD). In summary, the application of MPD in the MH algorithm may include six steps:

G. Wang, S. Chen / Ecological Modelling 253 (2013) 97–106

99

Table 1 Model parameters and their prior ranges. ID

Parameter name

Lower bound

Upper bound

Units

Remark

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20

rR KMB kN KMDN RDNmax kA kH kB kNO3 kNO2 kNO kN2 O KEM NO3 KEM NO2 KEM NO KEM N2 O Ks PE BE SDm

0.5 0.1 0.0001 0.1 0.1 0.0001 0.0001 0.0001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 1E−08 0.25 1 0.1

0.75 1000 0.1 100 100 0.01 0.01 0.01 10 10 10 10 10 10 10 10 1E−05 1.5 5 10

– g C m–3 h–1 g N m–3 g N m–3 h–1 g C g–1 CB h–1 g C g–1 CB h–1 g C g–1 CB h–1 h–1 h–1 h–1 h–1 g N m–3 g N m–3 g N m–3 g N m–3 m s–1 J kg–1 – –

Eq. (A.1c), Fraction of decomposed organic carbon that produces CO2 Eq. (A.2c), M–M half-saturation constant Eq. (A.5b), Rate constant of nitrification Eq. (A.7c), M–M constant for calculation of overall denitrification intensity Eq. (A.7c), Maximum overall denitrification intensity Eq. (A.2a), Specific decomposition rate in added organic matter (AOM) pool Eq. (A.2b), Specific decomposition rate in humus pool Eq. (A.2c), Specific death rate for biomass Eq. (A.7a), Rate constant for denitrification of nitrate Eq. (A.7a), Rate constant for denitrification of nitrite Eq. (A.7a), Rate constant for denitrification of nitric oxide Eq. (A.7a), Rate constant for denitrification of nitrous oxide Eq. (A.7b), M–M constant to calculate relative enzyme activity of nitrate Eq. (A.7b), M–M constant to calculate relative enzyme activity of nitrite Eq. (A.7b), M–M constant to calculate relative enzyme activity of nitric oxide Eq. (A.7b), M–M constant to calculate relative enzyme activity of nitrous oxide Saturated hydraulic conductivity Air entry potential Eq. (A.11), Exponent of the response function of soil water for decomposition Eq. (11), Multiplier for standard deviation of model simulation errors

(i) Generate parameter samples by MCMC using any of the four UPDs. (ii) Calculate the covariance matrix  from the parameter samples, and compute the eigenvalues and the corresponding eigenvectors. Let  be the eigenvector matrix whose columns are the normalized eigenvectors and let be the diagonal matrix whose diagonal entries are the eigenvalues in the order corresponding to the columns of . Calculate Q =



˚

(8)

(iii) Generate a vector x with all of its elements following the standard normal distribution N(0,1). (iv) Use a linear transformation to get x∗ = Qx + xk

(9)

where x* is the vector from the multivariate normal distribution Nd [xk , ]. (v) As mentioned before, the entire vector or a single element in the vector may be updated in each state transition. 2.4. Uncertainty assessment due to model parameters and structure

max{0,



y = SDm × 0.1 × y¯ obs

(11)

where y¯ obs is the mean value of yt , which is the square-roottransformed CO2 or N2 O flux (Yt ) in this study.

Soil GHG emissions after the application of liquid manure were measured on a dairy field in Pullman, WA (46◦ 45 N, 117◦ 11 W). The fluxes of CO2 and N2 O were calculated for the sampling period from August 28 to September 1 of 2006 (the period length is about 100 h). In order to verify the model simulations, an equally spaced (1 h) time series of CO2 and N2 O fluxes were generated from 13 observed data points by the Unit Response (UR) curve method (Wang et al., 2012a). Henceforth, the 101 data points (including one background flux before manure application) were used as observations to evaluate the model performance. Model initialization is presented in Table 2. 3. Results

2

Yt,sim () + N(0, y )}

A soil greenhouse gas emissions model (SoilGHG) was used to test the uncertainty analysis method. The major transformation processes in SoilGHG include decomposition, mineralization/immobilization, nitrification, and denitrification. A description of model equations is presented in Appendix A. In summary, 19 model parameters (parameters 1–19) are undetermined (Table 1). The prior parameter ranges were suggested through literature research (D’Odorico et al., 2003; Gu et al., 2004; Muller, 1999; Schimel et al., 1994). It is noted that the specific decomposition/death rates for the three pools (kA , kH , and kB ) have units of (g C g−1 CB h−1 ), which are different from the commonly used turnover rates in h−1 . The 20th “parameter” (SDm) in Table 1 is a multiplier for the standard deviation of model simulation errors. The standard deviation ( y ) of the simulation errors in Eq. (2) for CO2 or N2 O can be calculated by

2.6. Study site and model initialization

The 95% CIs for any output of interest due to parameter uncertainty can be computed by the output samples from a certain number of parameter sets generated by the MH algorithm. The 95% CIs for both parameter uncertainty and model structure uncertainty (PS uncertainty) were calculated by adding the model simulation errors to each of the output values at each time step, where the simulation errors due to structural uncertainty followed a normal distribution N(0, y ) with the standard deviation ( y ) determined in MCMC  (Wang et al., 2009). Since the square-root transformation Yt ) is used in this study, the GHG (CO2 or N2 O) flux with (yt = simulation errors, Yt,PS , is calculated by Yt,PS =

2.5. A test model for uncertainty analysis

(10)

where Yt,sim () is the simulated GHG flux at time t in terms of parameter . It is noted that the values of  or  y could change while the Markov chain is updated.

Comparing with the results from 400,000 iterations of the Markov chain, we found that the model runs with 40,000 iterations could generate almost the same Markov chain. This means that the model runs with 40,000 iterations were enough to obtain the posterior space of the parameters. There were totally 20 parameters

100

G. Wang, S. Chen / Ecological Modelling 253 (2013) 97–106

Table 2 Model initialization. Parameter

Value

Units

Remark

 CA CH CB CNA CNH CNB NAdd CNAdd CO2 NH4 NO3 NO2 NO N2 O N2 h

0.13 3800 34,000 700 19 13 10 28.25 19 0.0001 1 1 1 0.0001 0.0001 0.0001 41.3

m3 m–3 g C m–3 g C m–3 g C m–3 – – – g N m–2 – g C m–3 g N m–3 g N m–3 g N m–3 g N m–3 g N m–3 g N m–3 mm

Volumetric soil water content Carbon content in AOM pool Carbon content in humus pool Carbon content in biomass pool C:N ratio in AOM pool C:N ratio in humus pool C:N ratio in biomass pool Nitrogen applied to AOM pool C:N ratio in manures CO2 concentration in soil NH4+ concentration in soil NO− concentration in soil 3 NO− concentration in soil 2 NO concentration in soil N2 O concentration in soil N2 concentration in soil Liquid manure depth on soil surface

(Table 1) involved in the simulations and averagely each parameter had a chain with 2000 states, where the first 5% states were used for burn-in, and the other 95% were used for uncertainty analysis. 3.1. Posterior distributions of parameters by UPDs As shown in Fig. 1, the Bayesian inference allowed the identification of a narrower parameter range for each parameter than the prior space. For any UPD, all the posterior parameter ranges were reduced to be within 1 order of magnitude except four parameters (KMB , kN2 O, KEM NO2 , and KEM N2 O), which had acceptance rates above 50%. The acceptance rate is the fraction of proposed samples that is accepted (Hastings, 1970). The posterior range of Parameter 12 (kN2 O) was also reduced to be within a narrower range except for the asymmetric normal. Compared to the prior ranges of parameters, the posterior parameter ranges were averagely reduced by 67–82% for the four UPDs. In addition, the acceptance rate of symmetric uniform was lower than that of asymmetric uniform. As shown in Fig. 2, the average acceptance rate for univariate normal distributions (27%) was higher than that for univariate uniform distributions (22%).

To synthesize the results from the UPDs, the back-transformed mean values from the four UPDs was adopted for each parameter (Table 3). The average median and the 95% CI for each parameter were listed in Table 3. Parameter 20 (SDm) with an interval of (1.18–1.36) implied that the standard deviation was 11.8–13.6% (with a median of 12.4%) of the mean values of the square-root-transformed fluxes (y¯ obs ). In this study, the mean CO2 and N2 O fluxes were 2.12 kg C ha−1 h−1 and 0.0041 kg N ha−1 h−1 , respectively, thus the mean values (y¯ obs ) for the square-root-transformed CO2 and N2 O fluxes were 1.31 and 0.052. Therefore, the standard deviations of the simulation errors for y¯ obs of CO2 and N2 O fluxes had intervals of (0.15–0.18) and (0.006–0.007), respectively. 3.2. Posterior distributions of parameters by MPD The parameter samples from the asymmetric univariate normal distribution were used to estimate the covariance matrix. The matrix Q was computed using Eq. (8). Note that the covariance matrix is symmetric, but the matrix Q is asymmetric. For comparison, the results for MPD are also shown in Table 3. All the parameters’ posterior distributions were within one order of magnitude except KMB (Parameter 2 in Fig. 1). The 95% CI for Parameter 20 (SDm) was (1.33–1.78), corresponding to a range of standard deviation being 13.3–17.8% (with median of 15.2%) of the mean values (y¯ obs ). The total acceptance rate for MPD was 43% (Fig. 2), which was significantly higher than those of UPDs. The posterior frequency distributions of the logarithmic parameter values indicated that 14 parameters (3, 5–7, 9–12, 14–19) had unimodal-like distributions, three parameters (2, 4, and 20) showed a uniform-like distributions, and the other three parameters (1, 8, and 13) exhibited higher probabilities with the increase of parameter values within the parameter ranges. Overall, the posterior parameter spaces were much smaller than the priors (Fig. 1). 3.3. Uncertainty due to model parameters and structure The parameter samples generated by MCMC using MPD were used for uncertainty assessment. The 19 model parameters and the standard deviation of simulation errors for the model structure

Table 3 The posterior distributions of model parameters from univariate and multinormal proposal distributions. ID

Parameter name

Univariate proposal distributiona Median

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 a

rR KMB kN KMDN RDNmax kA kH kB kNO3 kNO2 kNO kN2 O KEM NO3 KEM NO2 KEM NO KEM N2 O Ks PE BE SDm

6.03E−01 9.33E+00 3.25E−03 6.71E+01 7.27E+01 3.49E−03 7.95E−04 9.09E−04 5.13E+00 2.98E−02 2.33E−02 9.21E−01 8.73E+00 4.46E−01 1.17E−02 3.83E−02 5.64E−06 9.67E−01 3.46E+00 1.24E+00

Multinormal proposal distribution

95% CI

Median

Lower bound

Upper bound

5.10E−01 1.16E+00 1.91E−03 6.17E+01 6.96E+01 2.77E−03 3.16E−04 3.52E−04 4.40E+00 1.84E−02 1.93E−02 7.33E−03 7.22E+00 1.20E−02 1.01E−02 3.47E−03 5.36E−06 8.48E−01 3.06E+00 1.18E+00

6.91E−01 1.53E+02 5.92E−03 7.04E+01 7.79E+01 4.49E−03 1.39E−03 1.18E−03 8.21E+00 3.94E−02 2.87E−02 2.46E+00 9.45E+00 1.34E+00 5.45E−02 2.80E−01 6.14E−06 1.09E+00 3.89E+00 1.36E+00

Median and 95% confidence interval (CI) are back-transformed mean values of four univariate distributions.

7.50E−01 6.42E+00 4.19E−03 7.09E+01 7.02E+01 3.34E−03 3.28E−04 7.65E−04 5.33E+00 4.35E−02 3.12E−02 1.29E−02 9.45E+00 4.74E−01 6.74E−02 1.57E−01 4.78E−06 1.27E+00 3.37E+00 1.52E+00

95% CI Lower bound

Upper bound

6.30E−01 1.02E+00 3.19E−03 6.49E+01 6.21E+01 2.72E−03 1.05E−04 1.20E−04 3.59E+00 2.71E−02 2.23E−02 8.09E−03 7.97E+00 3.30E−01 4.67E−02 1.14E−01 4.49E−06 1.10E+00 3.13E+00 1.33E+00

7.50E−01 1.41E+02 6.90E−03 8.22E+01 7.87E+01 4.10E−03 9.61E−04 9.99E−04 8.03E+00 5.96E−02 4.50E−02 1.75E−02 1.00E+01 5.96E−01 9.38E−02 2.81E−01 5.38E−06 1.43E+00 3.56E+00 1.78E+00

G. Wang, S. Chen / Ecological Modelling 253 (2013) 97–106

101

Fig. 1. Orders of magnitude of parameter ranges (each parameter is denoted by its ID as indicated in Table 1; Prior: prior ranges of parameter; Multinormal, Symmetric Normal, Asymmetric Normal, Symmetric Uniform, and Asymmetric Uniform denote the posterior ranges derived by MCMC using corresponding proposal distributions).

102

G. Wang, S. Chen / Ecological Modelling 253 (2013) 97–106

Acceptance rate (%)

50

4. Discussion

40

4.1. Symmetric vs. asymmetric proposal distribution

30 20 10 0

Symmetric Uniform

Asymmetric Uniform

Symmetric Normal

Asymmetric Normal

Multinorma l

Proposal distribution Fig. 2. Acceptance rates for different proposal distributions used in the Metropolis–Hastings algorithm.

(a)

8

CO2-C (kg C ha−1 h−1)

uncertainty were included in the analysis. The simulation errors were quantified by Eq. (11) with the parameter (SDm) shown in Table 1. As mentioned before, the model output (yt ) in Eq. (2) was assumed to follow a normal distribution. The 95% CIs of the simulated CO2 and N2 O fluxes due to parameter uncertainty and model structure uncertainty are shown in Fig. 3. The model parameter uncertainty had minimal impact on the CO2 simulations. However, the 95% CIs for CO2 simulations became much wider when the model structure uncertainty was included (Fig. 3a), which indicated that the uncertainty in CO2 emissions due to parameters was less important than due to model structure. In addition, Fig. 3a shows that almost all the observed data points were within the 95% CIs of the simulated CO2 fluxes when PS uncertainties were considered. The model simulations of N2 O underestimated the peak fluxes (Fig. 3b). The uncertainty in N2 O simulations due to model structure was not as significant as the structural uncertainty in CO2 simulations, which can be further verified by the analysis of the CI widths. For CO2 fluxes, the average CI width due to parameter uncertainty was 12% (ranging 7–23%) of that due to PS uncertainty, whereas this percentage increased to 30% (ranging 8–64%) for N2 O fluxes. Therefore, the parametric uncertainty in N2 O simulations was more significant than that in CO2 simulations.

6

Simulation (PS) Simulation (P) Observation

4

2

0 9

19

29

39

49 59 Time (h)

69

79

89

99

N2O-N (kg N ha−1 h−1)

(b) 0.03

The symmetric uniform distribution (J(xk |x* ) = J(x* |xk )) is commonly used in order to simplify the calculation of the acceptance probabilities. In addition, the asymmetric distribution and nonuniform distribution can also be employed in the MH procedure to investigate the performance and convergence of the algorithm. As for the posterior distributions of the parameters, the four UPDs led / J(x* |xk ), to similar results. Asymmetric distribution with J(xk |x* ) = on one hand, implies that more information (i.e., the proposal distribution) can be used for calculating the acceptance probability; on the other hand, it can improve the acceptance rate. However, the asymmetric distributions introduce additional complexity into the computations of the acceptance probability. Take the asymmetric uniform distribution as an example, as long as xk → 0, then x* → 0, eventually it would result in x* = 0, and the acceptance probability (Eq. (6)) could become incalculable. The same problem will occur for the asymmetric normal distribution. An alternative solution to this problem is to generate a new chain state from the parameter space instead of the current state. 4.2. Univariate vs. multivariate proposal distribution Compared with the UPD, the MPD is complex and requires more prior information about the parameters. In this study, although the MPD is a symmetric one, which may simplify the computation in the same manner as the symmetric UPDs, the covariance matrix of the parameters must be known in advance (Fang et al., 2004), and specific procedure is required for the sampling. Although either the entire vector or a single element may be updated during each step, the trials to simultaneously update the entire vector including 20 elements did yield an extremely low acceptance rate – nearly 0% in this study. It was then difficult to generate a new state with different values for the Markov chain. Therefore, we recommended to select a single element at random from the parameter-vector and update it at each step, which was the second method proposed by Hastings (1970). The total acceptance rate using the MPD was 43%, which was about twice of that from the UPDs. In most cases, the acceptance rate for a single parameter was higher with the MPD than with the UPD. This outcome might be attributed to that the covariance matrix from the UPDs was used, which could result in higher probability for the new points to be accepted. As shown in Fig. 1, the posterior parameter ranges from the MPD were usually narrower than those from the UPDs. The same priors were employed under conditions of UPDs and MPD. However, the posterior ranges of the three parameters (12, 14 and 16) were also reduced to be within one order of magnitude in case of MPD, whereas they were within 2–4 orders of magnitude in case of UPDs. A parameter vector generated by MPD in the sub-space near the current state has a higher probability to be chosen and to yield a higher likelihood.

Simulation (PS) Simulation (P)

0.02

4.3. Parametric and structural uncertainties vs. observation uncertainty

Observation

0.01

0 9

19

29

39

49 59 Time (h)

69

79

89

99

Fig. 3. The 95% confidence intervals for simulated GHG fluxes with uncertainty due to parameter (P) and both parameter and structure (PS): (a) CO2 and (b) N2 O.

The outputs typically include those data used for model calibration or model performance evaluation, e.g., the CO2 and N2 O fluxes in this study. The measured CO2 and N2 O fluxes in this study were averaged over 12 replicates (Wang et al., 2012a). The coefficient of variations (CVs) for the observed CO2 and N2 O fluxes were 0.45 and 0.55, respectively. When comparing the average CI widths due to parametric or PS uncertainty with the one-standarddeviation interval widths of measured fluxes, we found that the

G. Wang, S. Chen / Ecological Modelling 253 (2013) 97–106

parametric uncertainties were 16% and 23% of the observational uncertainties for CO2 and N2 O, respectively; and the PS uncertainties were 118% and 60% of the observational uncertainties for CO2 and N2 O, respectively. Therefore, the PS uncertainties in CO2 simulations were comparative to the observation errors, but the PS uncertainties in N2 O simulations were less than the observation errors. It suggests that reliable sampling and gas analysis procedures should be followed to control the measurement errors in GHG fluxes (Rochette and Eriksen-Hamel, 2008). 5. Conclusion The Bayesian inference and the MCMC simulations were applied to evaluate the uncertainties due to parameters and model structure of a soil GHG emission model. The normal UPDs could yield higher acceptance rates than the uniform UPDs. The sampling results from UPDs could also provide such essential information as covariance for the MCMC simulations with MPD. Nearly all the posterior parameter ranges (95% CIs) from MPD were reduced to 1 order of magnitude. We recommended updating a single parameter rather than to update the entire parameter vector at each time because the former is more efficient to update the chain. The simulation errors in CO2 fluxes were much greater than those in N2 O fluxes, which implied a greater importance of model structure than of model parameters for CO2 simulations. In conclusion, the uncertainties in parameters and structure of the SoilGHG model can be controlled within a reasonable range compared to the observation errors in GHG fluxes. In future studies, we may use more effective diagnosis methods (Chumbley et al., 2007; Marshall et al., 2005) to determine the convergence of the Markov chain. The method presented in this paper can be applied to assess the uncertainty of other ecosystem models. Acknowledgments The authors thank the Paul Allen Family Foundation and Climate Friendly Farm project for providing funding for this research, and the excellent editorial comments of Dr. Joan Wu. Thanks also go to the two anonymous reviewers for their constructive comments. Appendix A. SoilGHG model description The major transformation processes in SoilGHG include decomposition, mineralization/immobilization, nitrification, and denitrification. A description of model equations is presented as follows. Decomposition

dNB = dt

1 − rH

dCA = CAdd + CBD − CdecA dt

(A.1a)

dCH = rH CdecA − CdecH dt

(A.1b)

dCB = (1 − rR − rH )CdecA + (1 − rR )CdecH − CBD dt

(A.1c)

C CBD CdecA dNA = Add + − CNAdd CNB CNA dt

(A.1d)

CdecH CdecA dNH = rH − CNH CNH dt

(A.1e)

 Cdec

CNA CNH

A

CNA

+

CdecH CBD − − NnetM CNH CNB

(A.1f)

where t is time; the subscripts A, H, and B refer to AOM pool, humus pool, and biomass pool, respectively; CA , CH , and CB denote organic carbon contents in each pool (g C m−3 ); NA , NH , and NB are nitrogen concentration in each pool (g N m−3 ); CAdd is the external input to the system (g C m−3 h−1 ); CN is the C/N ratio; CdecA and CdecH are the decomposition rates of carbon (g C m−3 h−1 ); CBD is the carbon transformed back to the AOM pool due to the death of microbial biomass (g C m−3 h−1 ); rR is the fraction of decomposed organic carbon that produces CO2 ; rH is the fraction of CdecA that goes into the humus pool, and determined by min(0.25, CNH /CNA ) (D’Odorico et al., 2003); NnetM is the net mineralization (NnetM > 0) or net immobilization (NnetM < 0) (g C m−3 h−1 ). The decomposition of carbon and the death of microbial biomass is modeled using the Michaelis–Menten (M–M) kinetics (Wang et al., 2012b):













CdecA = kA fD (, T )

CB CA KMA + CA

(A.2a)

CB CH KMH + CH

(A.2b)

CB CB KMB + CB

(A.2c)

CdecH = kH fD (, T ) CdecB = kB fD (, T )

where k is the specific rate constant for decomposition or biomass death (g C g–1 CB h–1 ); fD (, T) is a dimensionless factor that defines soil water () and temperature (T) effects; KMB , KMA , and KMH are the half-saturation constants (g C m−3 ); and it is assumed that KMA = 5KMB and KMH = 50KMB . The mass balance of CO2 can be written as dCO2 = rR (CdecA + CdecH ) − ECO2 /Z dt

(A.3)

where CO2 is the concentration of CO2 in soil (g C m–3 ); ECO2 is the emission rate of CO2 from soil to air (g C m–2 h–1 ); Z is depth of the soil layer (m). Mineralization/immobilization The net mineralization/immobilization (NnetM ) is the difference between the rate of mineralization and the rate of immobilization in terms of the transformation between organic and inorganic nitrogen. The switch between the two states (net mineralization and net immobilization) and the calculation of NnetM depends on the assumption that the C/N ratio in the biomass pool remains constant in time (Manzoni et al., 2004; Porporato et al., 2003). NnetM = CdecH

Three functional pools are considered in SoilGHG, i.e., added organic matter (AOM, i.e., litter, fertilizer, or manure), humus, and microbial biomass (Porporato et al., 2003). The carbon and nitrogen balance equations for the three pools can be expressed as (Manzoni et al., 2004; Porporato et al., 2003):

103

 1

+ CdecA

CNH

−

 1 CNA

1 − rR CNB

−



rH 1 − rR − rH − CNH CNB

 (A.4)

where NnetM is controlled by the availability of NB and NH4+ . Nitrification The nitrification process is modeled by the first-order kinetics (Porporato et al., 2003): dNH4 = NnetM − NN dt

(A.5a)

NN = [kN fN (, T )] · NH4

(A.5b)

where NH4 is the concentration of NH4+ in soil (g N m–3 ); NN is nitrification rate (g N m–3 h–1 ); kN is the rate constant of nitrification (h–1 ); fN (, T) is a dimensionless factor that defines soil water and temperature effects on nitrification.

104

G. Wang, S. Chen / Ecological Modelling 253 (2013) 97–106

Denitrification Denitrification is the primary pathway for N2 O emission under wet conditions (Mathieu et al., 2006), which was the case in our study (after liquid manure application). The mass balance of nitrogen during the denitrification sequence can be expressed by dNO3 = NN − DNO3 dt

(A.6a)

dNO2 = DNO3 − DNO2 dt

(A.6b)

dNO = DNO2 − DNO − ENO/Z dt

(A.6c)

dN2 O = DNO − DN2 O − EN2 O/Z dt

(A.6d)

dN2 = DN2 O − EN2 /Z dt

(A.6e)

where NO3 , NO2 , NO, N2 O, and N2 represent the concentration of − nitrate (NO− 3 ), nitrite (NO2 ), nitric oxide (NO), nitrous oxide (N2 O), and nitrogen (N2 ), respectively (g N m–3 ); DNOx denotes the deni− trification rate (g N m–3 h–1 ) of a N oxide (NO− 3 , NO2 , NO, and N2 O), which is modeled using an enhanced competitive M–M kinetics; ENO, EN2 O, and EN2 are the emission rates of NO, N2 O, and N2 (g N m–2 h–1 ), respectively. The competitive M–M kinetics (Wang and Chen, 2012) assumes that the activity of the reductase changes with the availability of NOx following a M–M notation:





EMNOx =

MEM · NOx KEM + NOx

DNOx = kNOx · fDN , T

RDN = RDNmax · =

Nsum =

· EMNOx ·

Vden

max · Nsum KMDN + Nsum

RDNmax · Nsum KMDN + Nsum



Denom =



NOx



kNOx · NOx

RDN · NOx Denom

(A.7a) (A.7b)

of the thickness of a soil layer; GHGa and GHGs are gas concentrations in air and soil, respectively (g C m–3 or g N m–3 ), GHGa can be estimated by the method described in Muller (1999). Gas diffusivity in bimodal (aggregated) porous media can be expressed in terms of internal and external porosities (Jones et al., 2003; Potter et al., 1996):

 Ds = Da



εi / i

εi / i

+

2 

ε 2 e

e

εi / i − e

εi / i − e

2x0

2

2x1

1 − e 2x1

1 − e





2x

(εe − εe 2 ) 2x

+ (εe − εe 2 )

2 ε2x e

(A.9a)

where Ds denotes the gas diffusivity in soil (m2 h–1 ); Da is gas diffusivity in air (m2 h–1 ), for NO, N2 O, and N2 , Da = 0.05148 m2 h–1 at standard temperature and pressure (STP, i.e., 0 ◦ C, 101.3 kPa), for CO2 , Da = 0.05652 m2 h–1 at normal temperature and pressure (NTP, i.e., 20 ◦ C, 101.3 kPa), Da can be adjusted by temperature and pressure (Campbell and Norman, 1998); i and e are internal and external porosity, respectively, and = i + e ; εi and εe are internal and external air-filled porosity, and ε = εi + εe ; the exponents (x0 , x1 , and x2 ), representing the probability functions for tortuosity, are obtained by numerically solving the equation (Jones et al., 2003):

 2xj

j



= 1 − 1 − j

xj

, j = 0, 1, 2

(A.9b)

i is estimated as soil water content (SWC) at field capacity ( FC ), and e = −  FC (Potter et al., 1996). and  FC can be calculated from Saxton and Rawls (2006). Since micropores are mostly water-filled and macropores are mostly air-filled at  FC (Potter et al., 1996), i , e , εi , εe , and j can be estimated as follows: i = FC , e = − FC and s = 1 −



 /Vdenmax

εi =

FC − ,  < FC 0,  ≥ FC



(A.7c)

0 =



and εe =

(A.9c)

e ,  < FC



i −  / ( i + s ) ,  < i 0,  ≥ i

(A.9d)

e − ( − FC ),  ≥ FC

 , 1 = e ,

and 2 =

e + i − ,  > i e ,  ≤ i

(A.9e)

(A.7d) (A.7e)

where kNOx denotes the denitrification rate constant of NOx (h–1 ); fDN (, T) is a dimensionless factor that defines soil water and temperature effects on denitrification; EMNOx is a dimensionless parameter representing the relative enzyme activity; MEM (usually set to 1, the maximum relative reductase activity) and KEM (g N m–3 ) are the M–M parameters of NOx ; RDN is the overall denitrification intensity (g N m–3 h–1 ); RDNmax is the maximum RDN under optimum conditions (g N m–3 h–1 ); Vdenmax is the maximum growth rate for denitrifiers (h–1 ); KMDN is the M–M parameter for RDN calculation (g N m–3 ); Nsum is the total N oxide concentration (g N m–3 ); and Denom is the calculated denitrification intensity (g N m–3 h–1 ).

where  is the volumetric SWC. Response functions of soil temperature and water The effects of soil water and temperature can be expressed as the product of the response function of soil water and the response function of soil temperature (Wu and McGechan, 1998): f (, T ) = f () · f (T )

(A.10)

The influence of soil temperature for each transformation process was modeled by the Arrhenius equation (Gonc¸alves et al., 2007). The response functions of soil water for different transformation processes are introduced as follows. (i) Decomposition: The soil water response function decreases on either side of an optimum level (Johnsson et al., 1987):

Gas emissions Gas emissions were described by a diffusion process governed by Fick’s law (Campbell and Norman, 1998; Muller, 1999): EGHG = Ds · (GHGs − GHGa ) /z

2 

(A.8)

where EGHG is the gas (CO2 , NO, N2 O, or N2 ) emission rate (g C m–2 h–1 or g N m–2 h–1 ); Ds denotes the gas diffusivity in soil (m2 h–1 ); z is the diffusion distance (m), and is assumed as half

fD () =

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨



 − w l − w

BE

,

w ≤  < l

1, l ≤  ≤ h

⎪ ⎪  BE ⎪ ⎪ s −  ⎪ ⎪ , h <  ≤ s ⎩ es + (1 − es ) · s − h

(A.11)

G. Wang, S. Chen / Ecological Modelling 253 (2013) 97–106

where  s and  w are the saturated SWC and permanent wilting point, respectively;  h (= s − 0.08) and  l (= w + 0.1) are the high and low SWC, respectively; es denotes the relative effect of SWC at  s and BE is an exponent. (ii) Nitrification/denitrification: The influence of soil water can be elucidated by piecewise linear equations (Muller, 1999):

f (WFP) =

⎧ 0 WFP ≤ WFP1 ⎪ ⎪ ⎪ ⎪ ⎪ a + b · WFP WFP1 < WFP ≤ WFP2 ⎪ ⎨ + + 1

WFP2 < WFP ≤ WFP3

⎪ ⎪ ⎪ ⎪ ⎪ a− + b− · WFP WFP3 < WFP ≤ WFP4 ⎪ ⎩ 0

(A.12)

105

(ii) Asymmetric uniform distribution x∗





← U xk /ı, xk ı , ı > 1,



if xk > 0



(B.2a)

or x∗ ← U xk .ı, xk /ı , ı > 1 if xk < 0





  1 1  and J x∗ |xk =   x∗ · (ı−1/ı) xk · (ı−1/ı)       J xk |x∗ /J x∗ |xk = xk /x∗  , J xk |x∗ = 







    1 (x∗ − xk ) J xk |x∗ = J x∗ |xk = √ exp − 2 2 2

where WFP = / is water–filled pore space; a+ and b+ are intercept and slope of linear regression for increasing activity

(B.3a)



(B.3b)

(iv) Asymmetric normal distribution



∗



(B.2c)

(iii) Symmetric normal distribution ∗

x ← N xk , 

WFP > WFP4

(B.2b)

x ← N xk ,  k

  1 (x∗ − xk ) J xk |x∗ = √ exp − 2 2 k 2( k )

2







  1 (x∗ − xk ) J xk |x∗ = √ exp − 2 ∗ 2( ∗ )2

and

2



(B.4a) (B.4b)

If (i.e., b+ > 0); and a− and b− are intercept and slope of linear regression for decreasing activity (i.e., b– < 0). The values for WFPi , (i = 1,2,3,4), and a+ , b+ , a− , b− are adapted from Muller (1999). Basically, higher soil water content will inhibit nitrification and stimulate denitrification. A piecewise function is used for nitrification (e.g., Muller, 1999; Wu and McGechan, 1998), where f(WFP) increases with the increasing of SWC at lower WFP, but decreases with the increasing of SWC at higher WFP. f(WFP) for the denitrification of NO− 2 and NO follows a similar way to the f(WFP) for nitrification. As for the denitrification of NO− 3 and N2 O, f(WFP) increases with the increasing of WFP with a starting WFP greater than 0.35. Measured soil temperature at a 1–h interval was used as model input. One dimensional soil water movement was modeled by the Richards equation solved by the standard matric potential with Newton–Raphson method (Campbell, 1985). The soil profile was divided into 10 layers with unequal thickness, i.e., geometrically increasing node spacing as described in Campbell (1985). SoilGHG is a single-layer model, and this single layer is shallower than the whole soil profile. In this case, the weighted average water content over these m layers was used as input to the SoilGHG model, with the thickness of each layer used as weighting coefficient. At the upper boundary, the surface water infiltration and evapotranspiration were considered as a source term added to the first node. The lower boundary condition was set to a constant potential (Campbell, 1985). Appendix B. Proposal distribution functions in the Metropolis–Hastings algorithm Four univariate (symmetric and asymmetric uniform, symmetric and asymmetric normal) and one multivariate (symmetric multinormal) proposal distributions were examined. The ratios of the proposal distributions between the current state and the new state were derived for the asymmetric distributions. (i) Symmetric uniform distribution x∗ ← U[xk − ı, xk + ı]









J xk |x∗ = J x∗ |xk =

(B.1a) 1 2ı

(B.1b)

 k = xk

and  ∗ = x∗

then J(xk |x∗ )

k = ∗ · exp  J(x∗ |xk ) xk = ∗ · exp x





(B.4c)

(x∗ − xk ) 2( k )

 1−

2

(x∗ − xk ) 2(xk )

2

2

2

 1−

( k )

2



( ∗ )2 (xk )

2



(B.4d)

(x∗ )2

(v) Symmetric multinormal distribution



x∗ ← Nd xk ,



J(x∗ |xk ) = J(xk |x∗ ) = (2)−d/2 |

 exp

(B.5a)



|−1/2



−1 T 1 − (x∗ − xk ) (x∗ − xk ) 2

(B.5b)

where d is the number of elements in x* ; || is the determinant of ; and xk , the current state in MCMC, is the vector of means. References Campbell, G.S., 1985. Soil Physics with BASIC: Transport Models for Soil-Plant Systems. Elsevier Science Publishers B.V., Amsterdam, Netherlands. Campbell, G.S., Norman, J.M., 1998. An Introduction to Environmental Biophysics, 2nd ed. Springer-Verlag, New York, 286pp. Chen, D.L., Li, Y., Grace, P., Mosier, A.R., 2008. N2 O emissions from agricultural lands: a synthesis of simulation approaches. Plant and Soil 309, 169–189. Chumbley, J.R., Friston, K.J., Fearn, T., Kiebel, S.J., 2007. A Metropolis–Hastings algorithm for dynamic causal models. Neuroimage 38, 478–487. Congdon, P., 2001. Bayesian Statistical Modelling. Wiley, Chichester, NY, 531pp. D’Odorico, P., Laio, F., Porporato, A., Rodriguez-Iturbe, I., 2003. Hydrologic controls on soil carbon and nitrogen cycles. II. A case study. Advances in Water Resources 26, 59–70. Del Grosso, S., Ogle, S., Parton, W., Breidt, F., 2010. Estimating uncertainty in N2 O emissions from US cropland soils. Global Biogeochemical Cycles 24, GB1009. Engeland, K., Xu, C.-Y., Gottschalk, L., 2005. Assessing uncertainties in a conceptual water balance model using Bayesian methodology. Hydrological Sciences Journal 50, 45–63. Fang, S.F., Gertner, G.Z., Anderson, A.A., 2004. Estimation of sensitivity coefficients of nonlinear model input parameters which have a multinormal distribution. Computer Physics Communications 157, 9–16. Gonc¸alves, E.M., Pinheiro, J., Abreu, M., Brandão, T.R.S., Silva, C.L.M., 2007. Modelling the kinetics of peroxidase inactivation, colour and texture changes of pumpkin (Cucurbita maxima L.) during blanching. Journal of Food Engineering 81, 693–701.

106

G. Wang, S. Chen / Ecological Modelling 253 (2013) 97–106

Gu, L., Post, W.M., King, A.W., 2004. Fast labile carbon turnover obscures sensitivity of heterotrophic respiration from soil to temperature: a model analysis. Global Biogeochemical Cycles 18, GB1022. Hartig, F., Calabrese, J.M., Reineking, B., Wiegand, T., Huth, A., 2011. Statistical inference for stochastic simulation models – theory and application. Ecology Letters 14, 816–827. Hastings, W.K., 1970. Monte Carlo sampling methods using Markov Chains and their applications. Biometrika 57, 97–109. Hernadvolgyi, I.T., 1998. Generating Random Vectors from the Multivariate Normal Distribution. University of Ottawa, Canada, Ottawa, ON, pp. 1–14. Iskrev, N., 2007. How Much Do We Learn from the Estimation of DSGE Models? A Case Study of Identification Issues in a New Keynesian Business Cycle Model. Department of Economics, University of Michigan, pp. 1–49. Johnsson, H., Bergstrom, L., Jansson, P.E., Paustian, K., 1987. Simulated nitrogen dynamics and losses in a layered agricultural soil. Agriculture, Ecosystems & Environment 18, 333–356. Jones, S.B., Or, D., Bingham, G.E., 2003. Gas diffusion measurement and modeling in coarse-textured porous media. Vadose Zone Journal 2, 602–610. Kuehl, R.O., 1999. Design of Experiments: Statistical Principles of Research Design and Analysis, 2nd ed. Duxbury Press, Pacific Grove, CA, 666pp. Lehuger, S., Gabrielle, B., Oijen, M., Makowski, D., Germon, J.C., Morvan, T., Hénault, C., 2009. Bayesian calibration of the nitrous oxide emission module of an agroecosystem model. Agriculture, Ecosystems & Environment 133, 208–222. Li, C.S., Frolking, S., Frolking, T.A., 1992a. A model of nitrous oxide evolution from soil driven by rainfall events. 1. Model structure and sensitivity. Journal of Geophysical Research – Atmospheres 97, 9759–9776. Li, C.S., Frolking, S., Frolking, T.A., 1992b. A model of nitrous oxide evolution from soil driven by rainfall events. 2. Model applications. Journal of Geophysical Research – Atmospheres 97, 9777–9783. Link, W.A., Barker, R.J., 2008. Efficient implementation of the Metropolis–Hastings algorithm, with application to the Cormack–Jolly–Seber model. Environmental and Ecological Statistics 15, 79–87. Ma, L., Shaffer, M.J., 2001. A review of carbon and nitrogen processes in nine US soil nitrogen dynamics models. In: Shaffer, M.J., Ma, L. (Eds.), Modeling Carbon and Nitrogen Dynamics for Soil Management. CRC Press, Boca Raton, FL, pp. 55–102. Manzoni, S., Porporato, A., D’Odorico, P., Laio, F., Rodriguez-Iturbe, I., 2004. Soil nutrient cycles as a nonlinear dynamical system. Nonlinear Processes in Geophysics 11, 589–598. Marshall, L., Nott, D., Sharma, A., 2005. Hydrological model selection: a Bayesian alternative. Water Resources Research 41, W10422. Mathe, P., Novak, E., 2007. Simple Monte Carlo and the metropolis algorithm. Journal of Complexity 23, 673–696. Mathieu, O., Henault, C., Leveque, J., Baujard, E., Milloux, M.J., Andreux, F., 2006. Quantifying the contribution of nitrification and denitrification to the nitrous oxide flux using N-15 tracers. Environmental Pollution 144, 933–940. Muller, C., 1999. Modeling Soil–Biosphere Interactions. CABI Publishing, New York, NY, 354pp. Ogle, S.M., Breidt, F.J., Easter, M., Williams, S., Paustian, K., 2007. An empirically based approach for estimating uncertainty associated with modelling carbon sequestration in soils. Ecological Modelling 205, 453–463. Porporato, A., D’Odorico, P., Laio, F., Rodriguez-Iturbe, I., 2003. Hydrologic controls on soil carbon and nitrogen cycles. I. Modeling scheme. Advances in Water Resources 26, 45–58. Potter, C.S., Davidson, E.A., Verchot, L.V., 1996. Estimation of global biogeochemical controls and seasonality in soil methane consumption. Chemosphere 32, 2219–2246.

Refsgaard, J.C., 1997. Parameterisation, calibration and validation of distributed hydrological models. Journal of Hydrology 198, 69–97. Refsgaard, J.C., van der Sluijs, J.P., Hojberg, A.L., Vanrolleghem, P.A., 2007. Uncertainty in the environmental modelling process – A framework and guidance. Environmental Modelling & Software 22, 1543–1556. Reth, S., Hentschel, K., Drosler, M., Falge, E., 2005. DenNit – Experimental analysis and modelling of soil N2 O efflux in response on changes of soil water content, soil temperature, soil pH, nutrient availability and the time after rain event. Plant and Soil 272, 349–363. Rochette, P., Eriksen-Hamel, N.S., 2008. Chamber measurements of soil nitrous oxide flux: are absolute values reliable? Soil Science Society of America Journal 72, 331–342. Saxton, K.E., Rawls, W.J., 2006. Soil water characteristic estimates by texture and organic matter for hydrologic solutions. Soil Science Society of America Journal 70, 1569–1578. Schimel, D.S., Braswell, B., Holland, E.A., McKeown, R., Ojima, D., Painter, T.H., Parton, W.J., Townsend, A.R., 1994. Climatic, edaphic, and biotic controls over storage and turnover of carbon in soils. Global Biogeochemical Cycles 8, 279–293. Smith, P., Powlson, D.S., Smith, J.U., Elliott, E.T., 1997. Special issue – Evaluation and comparison of soil organic matter models using datasets from seven long-term experiments – Preface. Geoderma 81, 1–3. Stange, F., Butterbach-Bahl, K., Papen, H., Zechmeister-Boltenstern, S., Li, C.S., Aber, J., 2000. A process-oriented model of N2 O and NO emissions from forest soils. 2. Sensitivity analysis and validation. Journal of Geophysical Research – Atmospheres 105, 4385–4398. Thorsen, M., Refsgaard, J.C., Hansen, S., Pebesma, E., Jensen, J.B., Kleeschulte, S., 2001. Assessment of uncertainty in simulation of nitrate leaching to aquifers at catchment scale. Journal of Hydrology 242, 210–227. Tiana, G., Sutto, L., Broglia, R.A., 2007. Use of the Metropolis algorithm to simulate the dynamics of protein chains. Physica A (Statistical Mechanics and its Applications) 380, 241–249. Tonitto, C., David, M.B., Li, C.S., Drinkwater, L.E., 2007. Application of the DNDC model to tile-drained Illinois agroecosystems: model comparison of conventional and diversified rotations. Nutrient Cycling in Agroecosystems 78, 65–81. Wang, G., Chen, S., 2012. A review on parameterization and uncertainty in modeling greenhouse gas emissions from soil. Geoderma 170, 206–216. Wang, G., Chen, S., Frear, C., 2012a. Estimating greenhouse gas emissions from soil following liquid manure applications using a unit response curve method. Geoderma 170, 295–304. Wang, G., Post, W.M., Mayes, M.A., Frerichs, J.T., Jagadamma, S., 2012b. Parameter estimation for models of ligninolytic and cellulolytic enzyme kinetics. Soil Biology and Biochemistry 48, 28–38. Wang, G., Xia, J., Chen, J., 2009. Quantification of effects of climate variations and human activities on runoff by a monthly water balance model: a case study of the Chaobai River basin in northern China. Water Resources Research 45, W00A11. Wu, L., McGechan, M.B., 1998. A review of carbon and nitrogen processes in four soil nitrogen dynamics models. Journal of Agricultural Engineering Research 69, 279–305. Yohe, G., Oppenheimer, M., 2011. Evaluation, characterization, and communication of uncertainty by the intergovernmental panel on climate change – an introductory essay. Climatic Change 108, 629–639.