Soil organic matter (SOM) simulation models describe carbon (C) and nutrient - mainly nitrogen (N) - cycling in soil at different spatial (molecular to continental) and temporal (hours to decades) scales (Manzoni and Porporato, 2009). The use of SOM models (alone or integrated into more complex ecological, cropping and farming system models) spans many fields of application: e.g., the interpretation of SOM dynamics at micro-scale in laboratory incubation experiments (Van Veen et al., 1985; Blagodatsky et al., 1998), the estimation of longterm C dynamics in soils (Gijsman et al., 2002; Grace et al., 2006; Izaurralde et al., 2006), the quantification of N availability for plants (Coucheney et al., 2015) and N losses in the environment (Acutis et al., 2000), and the study of how SOM mineralization changes in response to a change of environmental variables (Davidson and Janssens, 2006; Conant et al., 2011; Liang et al., 2015; Sierra et al., 2015b).
Many SOM models are compartmental because they represent SOM heterogeneity through discrete pools, defined as homogeneous-in terms of composition and accessibility-fractions of SOM (Manzoni and Porporato, 2009; Sierra et al., 2011; Sierra and Müller, 2015). The number of pools varies from few to many. For example, in the ICBM model (Kätterer and Andrén, 2001) there are two pools, one of resistant (humified) and one of fresh organic matter, while in the BACWAVEWEB model (Zelenev et al., 2006) there are many pools of plant residues, microbial biomass, microfauna and native soil organic matter. Model pools are connected through a network of transformation/ decomposition fluxes that describe C and N turnover in soil. In the majority of the models (Manzoni and Porporato, 2009), SOM decomposition is simulated with first-order decomposition kinetics with respect to the decomposing pool (substrate). Instead, fewer models take simultaneously into account of the size of both the substrate and the microbial biomass (as a whole pool, an active fraction or as an exoenzymes pool) in simulating the SOM decomposition rate (Manzoni and Porporato, 2009). These models often assume that SOM decomposition is saturated with respect to substrate availability (Blagodatsky and Richter, 1998; Moorhead and Sinsabaugh, 2006; Ingwersen et al., 2008; Wutzler and Reichstein, 2008; Allison et al., 2010) or with respect to decomposers (Hadas et al., 1998; Schimel and Weintraub, 2003; Wutzler and Reichstein, 2008; Lawrence et al., 2009; Wutzler and Reichstein, 2013) using classical and reverse Monod (or Michaelis-Menten) equations, respectively.
Models in which decomposition is co-limited by substrate and decomposers have been used in simulating C priming effects in soils (Blagodatsky et al., 1998; Neill and Gignoux, 2006; Blagodatsky et al., 2010; Wutzler and Reichstein, 2013), i.e., a change in microbial biomass activity (apparent priming effects) or in the turnover of SOM (real priming effects) due to addition of fresh (exogenous) organic matter or nutrient to soil (Kuzyakov et al., 2000; Blagodatskaya and Kuzyakov, 2008; Kuzyakov, 2010).
Another important aspect concerning the simulation of SOM turnover is related to the effect that shortage of N has on C decomposition. In many models, each pool is supposed to have a constant C to N ratio (i.e., pools are strictly homeostatic), and N decomposition follows C decomposition stoichiometrically (Manzoni and Porporato, 2009). Thus, organic N can be mineralized to, or additional N can be immobilised from, soil mineral N (SMN), depending on the C to N ratios of the substrate and the microbial biomass, and to the substrate-C use efficiency (the fraction of decomposed C assimilated by decomposers) (Verberne et al., 1990). According to Manzoni and Porporato (2009), when SMN requirement (the potential N immobilisation) exceeds N availability, two different approaches can be employed to guarantee homeostasis of decomposers: the N inhibition hypothesis and the Coverflow mechanism. In the first case, a N-limitation factor is calculated such that C decomposition rates (and thus the overall SMN demand) are reduced, while in the second case, extra available-C is respired as CO2 (Schimel and Weintraub, 2003; Neill and Gignoux, 2006), or temporarily stored as extracellular polysaccharides (Hadas et al., 1998), without any reduction of SOM decomposition.
The implementation of Monod-type kinetics and the C-overflow mechanism increases model complexity and thus the number of parameters that must be initialized. Moreover, it is almost unavoidable that some parameters are subjected to calibration to achieve a good agreement between predictions and measurements (Müller et al., 2007). However, increasing model complexity can lead to undetermined model parameterisation and model overfitting (Sierra et al., 2015a). Thus, it is necessary to accurately select which model parameters have to be calibrated, according to available datasets. Parameters to be calibrated are normally those that have a major impact on the output variables of the model. To identify such parameters, sensitivity analysis (SA) can be used.
Sensitivity analysis studies how variation (uncertainty) in the output of a model can be apportioned to different sources of variation in model inputs (Saltelli et al., 2010). Sensitivity analysis is a necessary step of model development and a useful preliminary step to parameter calibration (Saltelli et al., 2004). Usually, SA requires running a model several times (the number of executions depending on the chosen method) by continuously changing the values of model parameters; these values are obtained by sampling parameters from pre-defined statistical distributions. Model outputs (resulting from the different model runs) are then analysed to assess the contribution of each parameter (and of their interactions) on model output variation. Among the different techniques for executing SA, global methods have been strongly recommended (Chan et al., 1997; Saltelli, 2002; Saltelli et al., 2004) because they are model-free (i.e., do not make any assumption about model structure and formulation) and evaluate model behaviour and output sensitivity in a wide range of parameters combinations (they are global). Global methods are suitable to rank parameters according to their sensitivity measures (Cariboni et al., 2007), and to select for calibration only the parameters to which model output is actually most sensitive (Cariboni et al., 2007; Saltelli et al., 2004). Examples of global SA of SOM models are not frequent (Cavalli and Bechini, 2011; Wang et al., 2005a, 2005b; Wang et al., 2013). Most published studies used simplified local methods, which change one parameter at a time in a specific point of the input space, while keeping the others constant (Schimel and Weintraub, 2003; Petersen et al., 2005a, 2005b).
Decomposition of organic input applied to soil is often studied with aerobic incubations conducted in the laboratory under constant temperature and soil water content. With such an approach it is possible to isolate unwanted environmental effects that affect SOM turnover and SMN dynamics. Results of laboratory aerobic incubations show that application of fresh organic matter to soil causes a temporary increase in CO2 emissions (Jensen et al., 2005; Morvan et al., 2006; Peters and Jensen, 2011); respiration of C usually follows a two-phase dynamic characterised by high emission rates in the first weeks and slower rates thereafter, when the more labile C in soil is exhausted. In soils amended with animal manures (Morvan et al., 2006; Bechini and Marino, 2009; Peters and Jensen, 2011) and low-N crop residues (Trinsoutrot et al., 2000; Jensen et al., 2005), net N immobilisation usually occurs immediately after amendment. The extent and duration of immobilisation depends mostly on the quality of added organic matter (in term of chemical recalcitrance, nutrient content and availability to microorganisms).
It is therefore very important to evaluate sensitivity of SOM decomposition at different time periods after soil amendment, and to consider organic inputs of different quality; this allows to identify the effect of model parameters on model output variation in relation to input quality, both in the short-term (about one-four weeks referring to laboratory conditions), medium-term (about four-eight weeks) and long-term (more than eight weeks).
The objective of this work was to modify a simple four-pool SOM model (ICBM/2BN, Kätterer and Andrén, 2001), replacing first-order with Monod-type decomposition kinetics, and introducing the C-overflow mechanism or the N inhibition hypothesis to limit C decomposition. Combining the modifications introduced in the original model, we had six variants of ICBM/2BN: a base model (MODFirst), two variants that implement classical and reverse Monod kinetics instead of first-order kinetics (MODMonod and MODRevMonod), and three variants equal to the three previous models but that implement the C-overflow mechanism instead of N inhibition hypothesis (MODFirst+Over, MODMonod+Over and MODRevMonod+Over). The six models were subjected to global SA, highlighting the contribution of model parameters on the variation of CO2 emissions and SMN concentration at different stages of SOM decomposition, from 3 to 189 days after soil amendment with animal manures or crop residues, under constant soil temperature and soil water content.
Materials and methods
The six models
The base model (MODFirst) builds on the ICBM/2BN model (Kätterer and Andrén, 2001), which represents soil organic matter with four pools: two pools of young organic matter (labile, YL, and resistant, YR), one of microbial biomass (B) and one of old organic matter (0). The structure of the base model is presented in Figure 1A, while model parameters and pools are listed in Table 1 and Table 2, respectively.
The base model represents young organic matter differently than Katterer and Andrén (2001); in fact, each input (e.g., animal manures and crop residues) is allocated to three young pools: a labile and a resistant one (YL and YR) as in ICBM/2BN, plus a very resistant pool (YHL). This additional pool represents the lignified fraction of crop residues (Henriksen and Breland, 1999; Hadas et al., 2004) and the partially humified organic matter of animal manures (Petersen et al., 2005b). Thus, a fraction of input-C is assigned to the labile pool (fCYL), another to the humus-lignin like pool (YHL), and the remaining to the resistant pool (fCYR = 1 – fCYL – fCO). Similarly, input organic N is partitioned among the three young pools with the parameters fNYL, fNYR and fNYHL. In the base model, soil (without added inputs: grey boxes in Figure 1A) is supposed to have some organic matter not only in the old pool (O) but also in the young resistant pool (YR soil), due to soil disturbance (Hadas et al., 2004).
Connections (decomposition and transformation fluxes) among model pools are presented in Figure 1A. Microbial biomass feeds on both young and resistant pools; a fraction (fB) of microbial residues is recycled through the living microbial biomass, while the remaining fraction (1 – fB) is assimilated to the pool O. During organic matter decomposition, microbial biomass assimilates only a fraction (eff) of the incoming C, and the other part (1 – eff) is released as CO2.
We modified the base model for the decomposition kinetics and the effects of N shortage on pool decomposition (Table 3).
Decomposition of organic pools
In the base model, all pools decay with first-order kinetics with respect to the decomposing pool. The outflow of C from each pool (pool decomposition rate, DEC) depends on the size of the substrate pool (CS) and on its decomposition constant (k):
A first model modification (MODMonod) assumes that decomposition of substrates CS depends also on the size of the microbial biomass (CB) pool, according to a classical formulation of the Monod equation (Blagodatsky et al., 1998):
Alternatively, a second model modification (MODRevMonod) is based on a reverse Monod equation (Eq. 3), which assumes that CS is sufficiently large, and that microbial biomass-rather than the substrate-saturates decomposition (Hadas et al., 1998):
Nitrogen limitation on decomposition
In all models, all pools are constrained to have a defined and constant C to N ratio; thus, microorganisms are considered homeostatic (Manzoni and Porporato, 2009). With this assumption, N decomposition follows C decomposition stoichiometrically, and N in each flux is equal to DEC divided by the C to N ratio of the decomposing pool. Moreover, the entire N in the flux (organic N) is considered available for microorganisms, according to the direct hypothesis for microbial N immobilisation (Manzoni and Porporato, 2009).
Considering that the C to N ratio of the receiving pool (e.g., microbial biomass, C/NB) can be different from the C to N ratio of the decomposing pool (C/NS), and that the C entering the receiving pool is different than the C leaving the decomposing pool, the N balance (ΔN) of each flux is calculated as (Verberne et al., 1990):
The C to N ratio of microbial biomass and eff allow calculating the critical C to N ratio (C/Ncrit) for microbial N immobilisation, which represents the substrate C to N ratio at which N immobilisation is equal to N mineralisation (ΔN = 0):
If ΔN is positive (Eq. 4), excess organic N is mineralised to ammonium, otherwise, if it is negative, the N required by the receiving pool is immobilised from the soil mineral N pool (SMN = ammonium-N + nitrate-N). Potential N immobilisation (IMMPOT) from all n-fluxes requiring extra N is then calculated as:
We first implemented the effect of SMN shortage (i.e., when IMMPOT is higher than available SMN) according to the N inhibition hypothesis (Manzoni and Porporato, 2007): decomposition fluxes having a negative ΔN are limited by a reduction factor (Nlim), ensuring that IMMPOT never exceeds SMN availability. In the definition of Nlim we assume that the maximum amount of N that can be immobilised is equal to a fraction (fmax = 0.95) of SMN, in order to avoid the complete depletion of SMN within one time step (Petersen et al., 2005b). The factor Nlim is calculated as:
Thereafter, DEC of all immobilising fluxes are multiplied by Nlim, thus reducing both C mineralization and SMN demand if Nlim < 1.
Alternatively to the N inhibition hypothesis, three variants of the model (MODFirst+Over, MODMonod+Over and MODRevMonod+Over; Figure 1B) implement the C-overflow mechanism (Manzoni and Porporato, 2009), only on the fluxes from young pools to B (Hadas et al., 1998). According to this hypothesis, when shortage of SMN occurs, DEC of immobilising fluxes is not limited by mineral N availability; instead, excess C is directed to polysaccharides (P) C-pool (Hadas et al., 1998). Therefore, excess of C (DECoverflow, i.e., C that cannot be assimilated in B under the current SMN availability) in each flux is calculated as (Manzoni and Porporato, 2009):Table 3).
Treatments simulated and model outputs recorded
The six models underwent SA in order to assess the sensitivity of C and N decomposition to a variation of model parameters. For each model, we performed separate SAs for the following treatments: i) unamended (bare) control soil (CON); ii) liquid manure (LM)-amended soil; iii) solid manure (SM)-amended soil; iv) low-N crop residue (CR)- amended soil (Table 4).
The selected model outputs of interest were CO2 and SMN concentration simulated in seven dates (after 3, 7, 14, 29, 44, 73, and 189 days after amendment). For amended treatments (LM, SM and CR), model outputs were calculated subtracting the values simulated in CON, to obtain the net effect.
Soil organic C and organic N were set to 11,065 mg C kg–1 and 1086 mg N kg–1, respectively, while inorganic soil N content was set to 7 mg N kg–1; these values were measured during an incubation experiment (unpublished data). The amount of C and N applied with LM, SM, and CR are reported in Table 4. Manures were applied at a rate of about 340 kg N ha–1, according to national regulations, while their chemical composition was set according to Cavalli et al. (2015). Crop residue was applied considering a maize grain production of 12 t dry matter (DM) ha–1, a grain DM to stalk DM ratio of 0.8 and a C concentration of maize stalk of 45%. For both manures and crop residue, we considered an incorporation depth of 0.3 m and a soil bulk density of 1.2 g cm–3.
We have run the simulations of SA by assuming constant soil temperature and water content, and, no response function (Kätterer and Andrén, 2001) was implemented.
Therefore, as a combination of six models, four treatments, two model outputs, and seven dates, we performed 336 SAs.
The method of Sobol’
As mentioned in the introduction, SA requires running the model many times by changing the values of parameters (which are sampled from a pre-defined statistical distribution).
The method of Sobol’ (Sobol’, 1993; Saltelli et al., 2004; Sobol’ and Kucherenko, 2005) is based on the decomposition of total model output variance (V) into different sources of variation: i.e., the partial variances due to single parameters (Vi), and all possible interactions of order w among model parameters (Vij, Vijm, etc.):
Sensitivity indices (S) are then calculated as the ratios of conditional variances (Vi, Vij, etc.) to total (unconditional) variance (Si = Vi/V; Sij = Vi/V, etc.), and represent the contribution to the total output variance of single parameters (Si) or combinations of parameters (Sij, etc.). Total-order indices (ST) represent the whole contribution of one parameter (alone and in combination with all other parameters) to total output variance (i.e., STi is the sum of all terms in Eq. 9 involving parameter i, divided by V). Given the rather high number of model parameters (from 13 to 20, Table 1), we calculated first-order (Si) and total-order (STi) sensitivity indices only. Thus, the difference between STi and Si, estimates the fraction of total output variance due to interactions (from order 2 up to order w) between parameter i and all other w-1 parameters.
Monte Carlo estimation of sensitivity indices was done according to Saltelli et al. (2010) using a sample size equal to 216.
Parameters distributions and models initialisation
During SAs, model parameters were sampled from uniform distributions, due to the lack of detailed a priori information about parameters variability. Lower and upper boundaries of each distribution were derived from values of parameters of similar significance used in other models, experimental measurements, or both (Table 1). When more than one parameter value was available in the literature, we set the limits of the distribution according to the minimum and the maximum values found. Conversely, when only one value was found, we set the lower and upper boundary of the distribution as ±50% of the reference value (Cavalli and Bechini, 2011).
Ranges of the decomposition and half-saturation constants of models MODMonod and MODMonod+Over were calculated from those imposed for models MODRevMonod and MODRevMonod+Over according to Eq. 10 and Eq. 11, respectively. In fact, Moorhead and Sinsabaugh (2006) gave evidence that models based on classical and reverse Monod kinetics would provide similar DECs under the following two constraints:
Limits for the resistant fraction of soil C (fCYR soil) were set by assuming that CYR soil represents 1-10% of soil C (Hadas et al., 2004). Microbial biomass C was initialised at 2% of soil C (Blagodatskaya and Kuzyakov, 2008; Wutzler and Reichstein, 2013). The microbial biomass C to N ratio was let to vary between 4 and 12 (Hassink, 1994; Friedel and Gabiel, 2001; Cleveland and Liptzin, 2007; Griffiths et al., 2012), by changing the N content of the pool B.
The obtained critical C to N ratio (Eq. 5) averaged 18 and was in agreement, with values commonly determined experimentally (in the range 15-35) for soils amended with crop residues (Trinsoutrot et al., 2000; Jensen et al., 2005), manure components (Van Kessel et al., 2000), and animal manure and organic fertilisers (Sørensen et al., 2003; Peters and Jensen, 2011; Delin et al., 2012).
The fraction of input organic N allocated to the pool YHL (fNYHL) was calculated in order to keep the C to N ratio of the pool YHL equal to 15. This value permitted to explore different partitioning of input C and N among the YL, YR and YHL pools, avoiding that the sum of fNYL, fNYR and fNYHL was higher than one, and ensuring that the C to N ratio of the pool YHL was close to that of O.
Variability of model outputs
Besides reporting sensitivity indices, we decided to show also the variation of model outputs over time. For each model output (CO2 or SMN, at one of the seven dates) we selected four reference points corresponding to the average, median, 10th and 90th percentile of the distribution of model output obtained after SA. Four reference points by seven dates yielded 28 model simulations. These 28 model simulations (different for CO2 and for SMN), for each treatment and model, were selected as representative of the variation in model outputs, and thus plotted for the entire incubation period.
Results of CO2 emissions (mg CO2-C kg–1) were divided by the amount of C applied with manures or crop residue (mg C kg–1), thus yielding the percentage of applied C that was mineralized (mg CO2-C 100 mg–1 of applied C). Results of net SMN (mg N kg–1) are presented as net N mineralisation (NNM), expressed as percentage of applied C. The NNM (mg N 100 mg–1 of applied C) represents the variation of net SMN compared to Day 0 (i.e., NNMday = i = (SMNDay = i – SMNDay = 0)/applied C). Thus, an increase of NNM indicates net mineralisation of organic N, while a decrease of NNM indicates net immobilisation of SMN. We have chosen to standardise NNM per unit of added C in order to be consistent with a way many experimental measurements of NNM are reported.
Variability of CO2 emissions and soil mineral nitrogen
The trend over time and the variability of simulated net CO2 and NNM obtained during SAs is reported in Figures 2-4.
Liquid manure amended soil
In the LM treatment, the implementation of the C-overflow mechanism produced no appreciable differences in net CO2 (Figure 2A) and NNM (Figure 2B). Conversely, the replacement of first-order with Monod kinetics slightly reduced the amount of mineralised C (4-8% of manure-C) in the first two weeks (Figure 2A). As incubation proceeded, net CO2 emissions raised slightly faster in MODMonod and MODMonod+Over than in other models. At Day 189, less C was mineralised in models with first-order and reverse Monod kinetics than in models with classical Monod kinetics (on average 6% and 11% of manure-C, respectively) in all simulations with the exception of those corresponding to 10th percentile of CO2 distributions. In such case, MODFirst and MODMonod decomposed the same amount of manure-C (on average 28%), and slightly more than MODRevMonod (23% of manure-C).
Until Day 14, with Monod kinetics respect to first-order kinetics NNM was higher (by 0.5-0.8% of manure-C) in simulations immobilising N (those identified with the 10th percentile points in Figure 2B), or lower (by 0.1-0.8% of manure C) in simulations mineralising N (average, median and 90th percentile points). Thereafter, models with firstorder and reverse Monod kinetics gave similar NNMs, while those implementing classical Monod equations mineralised more N (0.5-0.8% of manure-C at Day 189) in simulations corresponding to average, median and 90th percentile points, while no differences among models occurred in simulations corresponding to the 10th percentile points (Figure 2B).
Solid manure amended soil
Opposite of LM, replacing the N inhibition hypothesis with the Coverflow mechanism slightly enhanced manure-C and N mineralisation in SM (Figure 3). This occurred especially for models implementing first-order and classical Monod decomposition kinetics, with the exceptions of simulations corresponding to the 90th percentile of output distributions (Figure 3A and B). At the end of incubation (Day 189), extra C mineralised in MODFirst+Over, MODMonod+Over, and MODRevMonod+Over, compared to respective model variants implementing the N inhibition hypothesis (MODFirst, MODMonod, and MODRevMonod), corresponded to 3-12%, 3-8% and 2-5% of manure-C, respectively (Figure 3A). Similarly, results of NNM (Figure 3B) showed that extra N mineralised in models implementing the C-overflow mechanism increased as net N immobilisation rose: from 0.1-0.4% of manure-C at Day 189 for simulations corresponding to average and median, up to 0.5-0.8% of manure-C for simulations corresponding to 10th percentile of NNM distribution.
Replacing first-order with Monod decomposition kinetics produced different results over time for both net CO2 emissions and NNM. Until Day 14 (Figure 3A), mineralised C in MODFirst and MODFirst+Over was slightly higher compared to that of models with Monod (by 2-5% of manure-C) and reverse Monod kinetics (by 2-7% of manure-C). After Day 14, as previously observed for LM, CO2 emissions occurred at faster rates in MODMonod and MODMonod+Over than in other models, and at Day 189, net CO2 was higher, on average by 5 and 7% of manure-C, compared to models with first-order and reverse Monod kinetics, respectively (Figure 3A).
Differently from net CO2, NNM was higher in models with Monod kinetics than in those with first-order kinetics (by 0.1-1.0% of manure-C), and, in models implementing classical Monod compared to reverse Monod equations (by 0.1-0.4% of manure-C); moreover, differences among models amplified as incubation proceeded (Figure 3B).
Crop residue amended soil
In CR, the implementation of C-overflow mechanism always enhanced C mineralisation, even if the effect was narrow (+1-2% of residue-C at Day 189) for those corresponding to 90th percentile of CO2 distribution (Figure 4A). Conversely, for other points (average, median and 10th percentile), the increase at Day 189 corresponded to 7-19%, 4% and 2-3% of residue-C for models implementing first-order, classical, and reverse Monod kinetics, respectively. Conversely, the effect of Coverflow on NNM (Figure 4B) was relevant only for models with firstorder and classical Monod kinetics in simulations corresponding to 10th percentile points after Day 73. In fact, NNM at Day 189 was higher in MODFirst+Over, MODMonod+Over, and MODRevMonod+Over (by 0.6%, 0.3%, and 0.1% of residue-C, respectively), compared to respective model variants implementing the N inhibition hypothesis.
As observed for manured treatments, substituting first-order with Monod decomposition kinetics had an appreciable effect on both net CO2 emissions and NNM also in CR. At the end of incubation, in simulations corresponding to average, median and 90th percentile of net CO2 distribution, more C was mineralised, corresponding to 1-9% and 6-16% of residue-C, in models with first-order kinetics compared to those with classical and reverse Monod equation, respectively (Figure 4A). The opposite happened in simulations corresponding to 10th percentile of net CO2 distributions, when more CO2 accumulated in both MODMonod and MODRevMonod (on average 10% of residue-C) compared to MODFirst. During all the incubation, NNM was higher (+0.6-1.2% of residue-C at Day 189) in models with Monod kinetics than in those with first-order kinetics, while differences between models implementing classical and reverse Monod equations were extremely small (very frequently lower than 0.1% of residue-C).
After SAs, parameters affecting C and N turnover were considered important if they had a first-order sensitivity index higher than 10% for at least one of the output variables (CO2 or SMN) in at least one date. The mean and the standard deviation of first-order sensitivity indices were calculated across the seven dates, separately for each treatment and output variable (Table 5). Moreover, a qualitative trend of the indices over time was given, separately for the short-medium term (Day 3-44) and medium-long term (Day 44-189). The trend was assigned three categorical values: increasing, decreasing, and constant. A value, within each time period, was considered constant if its absolute variation (max – min) was less than 5%.
In CON, parameters that mostly affected CO2 emissions were those related to the turnover of the young resistant pool (fCYR soil and kYR soil) and the substrate use efficiency (eff); these three parameters account together for a relevant fraction of CO2 variance (from 45 to 82%, depending on the model). In addition, in models implementing reverse Monod kinetics, the decomposition constant of microbial biomass (kB) and the half-saturation constant (krm) accounted for a further 29% of CO2 variance, while in those implementing classical Monod kinetics, the half-saturation constant (kmYRsoil) alone explained 26% of CO2 variation. The importance of microbial biomass parameters eff and kB usually decreased with time, especially after Day 44, while that of the half-saturation constants increased. Conversely, the effect of fCYR soil increased after Day 44 in all models, while the trend of kYR soil depended on the type of decomposition kinetics: Si decreased after Day 44 in MODFirst and MODFirst+Over, while it remained constant until Day 189 in all other models.
The three microbial biomass parameters eff, fNB and, with the exception of MODMonod and MODMonod+over, kB, explained a relevant fraction also of the variation of SMN (from 56 to 69% depending on the model). As observed for CO2, first-order indices of these parameters usually tended to decrease after Day 44.
As expected, in the amended treatments parameters related to the turnover of pools CYR soil and O had no appreciable effect in any model (first-order indices always less than 10%) on both net CO2 emissions and net SMN. Microbial biomass parameters had an effect similar to what was observed in CON. In addition, in the amended treatments some input-related parameters markedly contributed to net CO2 and net SMN variances. With few exceptions, average first-order indices of the labile fraction of input C (fCYL) were high for both net CO2 (5-25%) and net SMN (9-23%) in all models. Moreover, the importance of fCYL usually diminished as time passed, and as the C to N ratio of the input raised (especially from SM to CR).
Other input parameters that had an effect on net CO2 variation were the decomposition constants of YL and YR pools (kYL and kYR) that accounted together for 9-40% of net CO2 variance. Similarly to fCYL, the effect of first-order indices of kYL decreased with increasing the C to N ratio of the input; the opposite happened for kYR, especially for manures and CR, with few exceptions. Also the trend over time of Si differed between kYL and kYR: the former had a decreasing effect, while the latter had an increasing effect over time. Differently from what was observed for net CO2, sensitivity of net SMN to kYL was relevant in SM (for models implementing classical and reverse Monod kinetics), while sensitivity to kYR was relevant in CR (for all models). The effect of kYR usually increased until Day 44 and then stabilised thereafter.
In models with Monod kinetics, also half-saturation constants were important (krm, kmYL and, with the exception of SM, kmYR). These parameters gave an average contribution of 8-25% to net CO2 variance and, solely for CR, of 2-8% to net SMN variance. For both output variables, sensitivity indices of krm and kmYL decreased until Day 44, and then remained constant until Day 189, while those of kmYR usually increased during time.
The sum of first-order sensitivity indices (ΣSi) accounted for 70-94% and 57-88% of total output variance for CO2 and SMN, respectively (Table 5), showing that interactions among model parameters (1 – ΣSi) considerably contributed to total output variance in some model x input combinations, while in other combinations the interactions were not important. Moreover, interactions among parameters increased with increasing the C to N ratio of the input, especially in models implementing the N inhibition hypothesis (Table 5).
In most cases, parameters with high interactions (Table 6) were the same that had high first-order indices. In CON the most important interactions affecting CO2 variance were for parameters fCYR soil, kYR soil and the half-saturation constant kmYR soil. Interactions of microbial biomass parameters eff and fNB also markedly contributed to SMN variance. Besides eff and fNB, in the amended treatments net CO2 and net SMN were also sensitive to interactions of input parameters with high Si (fCYL, kYL, kYR, krm, kmYL, and kmYR) in many model × input combinations.
In this section, we first compare the simulated dynamics of CO2 emissions and NNM with published experiments. Even if we did not calibrate the six models to reproduce a specific set of measurements, it is important to verify if the simulated statistical distributions and the trends over time are in agreement with real data. Thereafter, we discuss the parameters to which model outputs are most sensitive. We then compare models with respect to N limitation effects on C decomposition, and with respect to decomposition kinetics. Finally, we emphasize that some models were able to simulate priming effects.
Variability of simulated CO2 emissions and net nitrogen mineralisation
Dynamics of simulated net CO2 emissions in amended treatments (Figures 2A, 3A, and 4A) were in general agreement with those observed in aerobic incubations experiments of soil amended with liquid animal manures (Sørensen and Fernández, 2003; Sørensen et al., 2003; Morvan et al., 2006; Bechini and Marino, 2009; Cavalli et al., 2014), solid animal manures (Thomsen and Olesen, 2000; Calderón et al., 2005; Morvan et al., 2006; Peters and Jensen, 2011) and low-N crop residues (Henriksen and Breland, 1999; Trinsoutrot et al., 2000; Nicolardot et al., 2001; Hadas et al., 2004; Jensen et al., 2005; Redin et al., 2014).
Carbon mineralisation dynamics in SL, SM and CR were described by two phases of progressively decreasing decomposition rates. The first phase, lasting for about two months, was characterised by high and exponentially decreasing CO2 emissions, and was accompanied by the complete depletion of the more labile input organic matter in 14-73 days (YL); thereafter, decomposition of the more resistant input organic matter (YR) mostly contributed total CO2 emission (Figures 2A, 3A and 4A). During the whole incubation, only a fraction of initial YR (36-70%, 42-67%, and 33-61% in LM, SM, and CR, respectively) was mineralised (data not shown). Conversely mineralisation of YHL was insignificant (1-3% of input C) for all inputs.
In our simulations, the application of exogenous organic matter to soil often induced temporary net immobilisation of N (up to several weeks) by soil microbial biomass, especially when the C to N ratio of the added organic input was far higher than that of the soil microbial biomass, as demonstrated experimentally for solid manures (Calderón et al., 2005; Morvan et al., 2006; Peters and Jensen, 2011) and low-N crop residues (Henriksen and Breland, 1999; Trinsoutrot et al., 2000; Nicolardot et al., 2001; Hadas et al., 2004; Jensen et al., 2005; Redin et al., 2014). This pattern was well represented by all models after application to soil of both SM (Figure 3B) and CR (Figure 4B).
Conversely, the simulated NNM by all models was excessive, already from the beginning of the incubation, when compared to measurements carried out in incubation studies with LM. In fact, even in the case of liquid manures, with low C to organic N ratios, considerable net N immobilisation is often experimentally observed in the first weeks of manure decomposition (Sørensen and Fernández, 2003; Sørensen et al., 2003), frequently lasting for months (Bechini and Marino, 2009). On the contrary, in our study we observe that only model simulations corresponding to the 10th percentile of net SMN distributions provided manure-N mineralisation similar to what observed in reality, while the rest of the distributions were far too high (Figure 2B). To explain this fact, we point our attention to the ranges of model parameters used in SAs. These ranges were chosen in order to simulate both C and N decomposition. However, Cavalli and Bechini (2012) have shown that, with models characterised by fixed microbial C to N ratio and C use efficiency, a trade-off exists between an accurate simulation of C and N mineralisation. In particular, to simulate the immobilisation of N that is commonly observed in experiments, it is necessary to simulate a very low C respiration, a situation that, with the ranges of parameter values that we have chosen, occurred rarely.
Parameters common to all models affecting CO2 emissions and net nitrogen mineralisation
Despite the rather high number of parameters needed in the six models (Table 1), respired C and SMN were sensitive only to a variation of few of them (Table 5) and their interactions (Table 6). In all models, substrate use efficiency (eff) accounted (on average for the entire period) for 14-54% of CO2 variance. The high sensitivity of CO2 to eff was due to its direct effect on C turnover (Figure 1), regulating the amount of C respired as CO2, and its feedback on C decomposition by defining N required by microbial biomass (Eq. 4). Thus, when shortage of N occurs, CO2 emissions are reduced proportionally to the value of eff by diminishing decomposition rates (N inhibition hypothesis, Eq. 7) or by temporary storing excess-C in the polysaccharides pool (C-overflow mechanism, Eq. 8). In addition to eff, concentration of SMN (Tables 5 and 6) was also sensitive to a variation of fNB (and thus to the C to N ratio of soil microbial biomass). Together, eff and fNB explained from 26 to 63% of SMN variance; this is justified by their effect in determining N requirements by microbes (Eq. 4).
In the amended treatments, neither the most recalcitrant fraction of input C (fCYHL) nor its decomposition constant (kYHL) gave effects on net CO2 emissions and SMN in all models. This can be attributed to range that we have chosen for kYHL (Table 1) and to the reduced simulation time span that made the YHL pool almost inert (data not shown). Conversely, C decomposition and SMN were very sensitive to variations of the size (fCYL) and decomposition constant (kYL) of the labile young pool YL (Tables 5 and 6). Both parameters influence the mineralisation of labile C; moreover, fCYL is fundamental in regulating input-N availability by switching from N-rich to N-poor YL pools (Tables 1 and 2), and thus strongly impacts SMN dynamics.
Effects of carbon-overflow on CO2 emissions and net nitrogen mineralisation
Differences between models implementing N inhibition hypothesis and C-overflow mechanism occurred in SM and CR (Figures 3A and 4A), but not LM (because for LM N was never limiting C decomposition). In SM and CR, models implementing C-overflow decomposed more input-C (on average 8% and 10% of added C at Day 189 in SM and CR, respectively) and immobilised less N (on average –0.7% and –0.4% of added C at Day 189 in SM and CR, respectively) than models implementing the N inhibition hypothesis. The simulated dynamics of the polysaccharides pool were in agreement with that of Hadas et al. (1998). For a soil amended with a low-N straw (C to N ratio of 91), they simulated a net growth of the polysaccharides pool until N limitation persisted (Day 50-80). Thereafter, due to ceased accumulation, the polysaccharides declined, with a velocity depending on kP. In our work net C accumulation in the polysaccharides pool lasted until Day 44-73 in SM and Day 189 in CR (Figure 5), when the maximum size of P was reached. The maximum size of the polysaccharides pool (Figure 5) was similar for MODFirst+over and MODMonod+Over(15-16% and 18-23% of inputC in SM and CR, respectively) and was lower for MODRevMonod+Over (11% and 9% of input-C in SM and CR, respectively). Results of MODMODFirst+over and MODMonod+Over are in agreement with those of Hadas et al. (1998) that reported a maximum accumulation of C in the pool P of about 2225% of input-C.
However, sensitivity of CO2 emissions to the parameter describing polysaccharides decomposition rate (kP) was nil (Tables 5 and 6) because the amount of polysaccharides-C decomposed after Day 73 (Figure 5) was relatively low in SM (2-6% of manure-C) and almost zero in CR. Thus, the effect of previously mentioned parameters (eff, fCYL and kYL) masked that of kp in determining CO2 variation. Different results would have been obtained if all C in excess was respired as CO2, without a temporary storage in P (Schimel and Weintraub, 2003; Neill and Gignuoux, 2006). However, the polysaccharides pool was specifically added to avoid strong increase of C respiration under N-limited conditions (Hadas et al., 1998), even if both mechanisms (i.e., exudation of organic molecules and extra CO2 production) were proved to occur (Russel and Cook, 1995).
Effects of decomposition kinetics on CO2 emissions and nitrogen mineralisation
In models implementing classical and reverse Monod kinetics, C mineralisation, and to a less extent SMN, were sensitive to a variation of the half-saturation constants kmYR soil, kmYL, and kmYR (MODMonod and MODMonod+Over), and krm (MODRevMonod and MODRevMonod+Over) (Tables 5 and 6) due to their direct effect on C decomposition rates (Eqs. 2 and 3).
In the first days of incubation (until Day 7, 14 and 29 in LM, SM and CR, respectively) the general impact of decomposition kinetics on CO2 emissions was similar for all inputs, and consisted in a reduction (on average by 0.5) of C mineralisation with Monod and reverse Monod compared to first-order kinetics (Figures 2A, 3A, and 4A). Reduced C decomposition induced less N mineralisation or immobilisation (Figures 2B, 3B, and 4B).
After Day 73, differences among models depended on input type and on the amount of decomposed C and N. When microbial biomass grows without long-term N limitations, models with classical Monod kinetics decomposed more C and N than the other models. This was evident in LM and SM, excluding those simulations where CO2 emissions were very low (10th percentile), when C mineralisation was on average +6% and +10% of manure-C in models with classical Monod kinetics compared to first-order and reverse Monod kinetics, respectively (Figures 2A and 3A). The increase in CO2 emissions in MODMonod and MODMonod+Over was due to the growth of soil microbial biomass that in turn enhanced C decomposition rates (k’ × CR, Eq. 2), while the effect of substrate exhaustion (Cs/(km+Cs), Eq. 2) was lower. Therefore, after Day 73, NNM was higher (on average +0.6% of manure-C) in models implementing classical Monod kinetics than in others, in those cases when input decomposition did not induce long-term net N immobilisation (Figure 2B and 3B).
Conversely, when substrate decomposed slowly, and microbial biomass was C- or N- limited (as for CR), models with first-order kinetics decomposed more C (on average 8% of manure-C at Day 189) and immobilised more N (on average 1.0% of manure-C at Day 189) than models with Monod and reverse Monod kinetics.
An extreme case occurred for MODFirst that showed negative net CO2 emissions after Day 29. This result was an artefact of model parametrization. In these simulations native soil pools (YR soil and O) had a C to N ratio higher than C/Ncrit (Eq. 5) and thus induced net N immobilisation during decomposition. The application of CR to soil caused additional N immobilisation and generated the onset of N limitation. This did not happen in the unamended soil, which therefore released more CO2 than the treated soil.
Simulation of priming effects
A very important finding is that some of the models were able to simulate negative or positive priming effects (PE). Indeed differences among models for net CO2 emissions and net SMN were not only due to a different turnover of input pools but also of native soil pools.
In SM and CR, when the C to N ratio of YR soil was higher than C/Ncrit and SM and CR exhausted the available mineral N causing N limitation, the models implementing the N inhibition hypothesis decomposed less YR soil than in CON. This negative PE was small for models implementing Monod kinetics (–1 and –3% of added C at Day 189 in SM and CR, respectively), while it was higher for MODFirst (–10% and – 21%). In addition, models implementing Monod kinetics (irrespective of N control mechanism on C decomposition), always showed positive PEs on YR soil when decomposition of native SOM was not N-limited. This is represented in Figure 6 (only for models with C-overflow) with negative net CYR soil, i.e., less CYR soil at the end of the incubation in the amended treatments than in CON.
Increasing the amount of added-C, more YR soil was decomposed (Figure 6A), in particular for LM and SM. Decomposition of YR soil increased as the input C to organic N ratio decreased (Figure 6B); this was in accordance with the assumption that as more input-C was decomposed (in the order LM > SM > CR), more microbial biomass grew and thus more YR soil was mineralised. Thus, in this case, the input C to organic N ratio is a proxy of input decomposability.
Moreover, distribution statistics of net YR soil C showed that PEs varied from few percent points to –27% and –14% of added C, for models implementing classic and reverse Monod kinetics, respectively (Figure 6B). Variability of PE within each added input was also due to differences in decomposability, resulting from the partitioning of added C and N among young pools: high C decomposition induced high microbial growth that decomposed more YR soil, while low, often N-limited, microbial growth did not stimulate YR soil decomposition, and thus did not produce appreciable PEs.
Models implementing Monod kinetics predicted PE according to cometabolism of fresh (YL, YR) and native (YR soil) organic matter (Kuzyakov et al., 2000; Wang et al., 2015), under constraints of stoichiometric decomposition theory (Hessen et al., 2004; Craine et al., 2007; Chen et al., 2014). According to this theory, high-quality substrates (high accessibility, fast decomposition, and with C to N ratio close to C/Ncrit) induce balanced microbial growth that in turn enhances decomposition of native SOM. Our results were in agreement with this theory. The size of PE was lower in CR compared to LM and SM, due to lower decomposition rates and higher C to N ratio of CR. Moreover, PE was lower in simulations when microbial biomass was C or N-limited following LM, SM or CR additions (Figure 6B).
Our models do not implement the mechanisms of priming effect proposed by two published theories, and thus our results are not in agreement with these theories. First, predicted PEs (mg C kg–1) increased linearly with C addition (Figure 6A). While this is in accordance with experimental results for low C additions, a negative relationship was demonstrated between applied C and primed C when the ratio of added-C to microbial biomass C exceeds about 50 (high C additions), possibly due to preferential substrate utilisation (Blagodatskaya and Kuzyakov, 2008). The preferential substrate decomposition theory (Wang et al., 2015) states that a negative PE might occur when microorganisms preferentially decompose the added organic materials, rather than the native soil organic matter. In our analysis, models did not predict this negative PE because no preferential decomposition mechanisms are implemented in our models.
Second, predicted PE was not in agreement with the microbial Nmining theory (Moorhead and Sinsabaugh, 2006; Craine et al., 2007; Chen et al., 2014; Wang et al., 2015), stating that soil microbial biomass feeds on added fresh organic matter to get energy, and decomposes native SOM to release limiting nutrients (N). Based on this theory, we should have obtained higher PE when the shortage of N increased (i.e., in CR compared to manures, and in simulations with N shortage). However, our models do not implement differential mechanisms for C and N mineralisation, and thus this effect could not be obtained.
Model modifications (introduction of C-overflow and Monod decomposition kinetics) had an effect on simulated CO2 emissions, soil mineral nitrogen concentration, and turnover of added (manures and crop residues) and native soil organic matter pools. In particular, models implementing Monod decomposition kinetics were able to simulate positive priming effects after application of all inputs. Moreover, sensitivity analysis indicated that simulated CO2 and SMN were sensitive to a variation of only few model parameters; these are the parameters that should be selected for model calibration using experimental measurements.