## Introduction

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 CO_{2} (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 CO_{2} 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 (MOD_{First}), two variants that implement classical and reverse Monod kinetics instead of first-order kinetics (MOD_{Monod} and MOD_{RevMonod}), and three variants equal to the three previous models but that implement the C-overflow mechanism instead of N inhibition hypothesis (MOD_{First+Over}, MOD_{Monod+Over} and MOD_{RevMonod+Over}). The six models were subjected to global SA, highlighting the contribution of model parameters on the variation of CO_{2} 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

#### Model structure

The base model (MOD_{First}) 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, *Y _{L}*, and resistant,

*Y*), one of microbial biomass

_{R}*(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 (*Y _{L}* and

*Y*as in ICBM/2BN, plus a very resistant pool (

_{R})*Y*). This additional pool represents the lignified fraction of crop residues (Henriksen and Breland, 1999; Hadas

_{HL}*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 (

*fC*), another to the humus-lignin like pool (

_{YL}*Y*), and the remaining to the resistant pool (

_{HL}*fC*= 1 –

_{YR}*fC*–

_{YL}*fC*Similarly, input organic N is partitioned among the three young pools with the parameters

_{O}).*fN*and

_{YL}, fN_{YR}*fN*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 (

_{YHL}.*Y*), due to soil disturbance (Hadas

_{R}soil*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 (*f _{B}*) of microbial residues is recycled through the living microbial biomass, while the remaining fraction (1 –

*f*is assimilated to the pool O. During organic matter decomposition, microbial biomass assimilates only a fraction

_{B})*(eff)*of the incoming C, and the other part (1 –

*eff)*is released as CO

_{2}.

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 (*C _{S}*) and on its decomposition constant (k):

A first model modification (MOD_{Monod}) assumes that decomposition of substrates *C _{S}* depends also on the size of the microbial biomass (

*C*) pool, according to a classical formulation of the Monod equation (Blagodatsky

_{B}*et al.*, 1998):

*k’*is the decomposition constant,

*k*m is the half-saturation constant which indicates the size of

*C*at which

_{S}*DEC*proceeds at half rate compared to first-order.

Alternatively, a second model modification (MOD_{RevMonod}) is based on a reverse Monod equation (Eq. 3), which assumes that *C _{S}* is sufficiently large, and that microbial biomass-rather than the substrate-saturates decomposition (Hadas

*et al.*, 1998):

*k*is the half-saturation constant for the reversed Monod models.

_{rm}#### 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/N _{B})* can be different from the C to N ratio of the decomposing pool

*(C/N*, and that the C entering the receiving pool is different than the C leaving the decomposing pool, the N balance (Δ

_{S})*N*) of each flux is calculated as (Verberne

*et al.*, 1990):

*eff*is the substrate assimilation efficiency (equal to one for the flux from

*B*to

*O).*

The C to N ratio of microbial biomass and *eff* allow calculating the critical C to N ratio *(C/N _{crit})* 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 *(IMM _{POT})* from all

*n*-fluxes requiring extra N is then calculated as:

We first implemented the effect of SMN shortage *(i.e.*, when *IMM _{POT}* 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 (

*N*), ensuring that

_{lim}*IMM*never exceeds SMN availability. In the definition of

_{POT}*N*we assume that the maximum amount of N that can be immobilised is equal to a fraction (

_{lim}*f*= 0.95) of SMN, in order to avoid the complete depletion of SMN within one time step (Petersen

_{max}*et al.*, 2005b). The factor

*N*is calculated as:

_{lim}Thereafter, *DEC* of all immobilising fluxes are multiplied by *N _{lim}*, thus reducing both C mineralization and SMN demand if

*N*< 1.

_{lim}Alternatively to the N inhibition hypothesis, three variants of the model (MOD_{First+Over}, MOD_{Monod+Over} and MOD_{RevMonod+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 *(DEC _{overflow}, i.e.*, C that cannot be assimilated in

*B*under the current SMN availability) in each flux is calculated as (Manzoni and Porporato, 2009):

*IMM*(equal to –

_{MAX}*ΔN*×

*N*is the maximum amount of SMN that can be immobilised. Therefore, decomposed C assimilated by

_{lim})*B*is higher compared to that assimilated under the N inhibition hypothesis. Decomposition of the polysaccharides pool is limited by N according to the N inhibition hypothesis (Table 3).

### Sensitivity analysis

#### 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 CO_{2} 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 (*V _{i}*), and all possible interactions of order

*w*among model parameters (

*V*,

_{ij}*V*, etc.):

_{ijm}*m*, ..., and

*w*represent model parameters.

Sensitivity indices (S) are then calculated as the ratios of conditional variances (*V _{i}, V_{ij}, etc.)* to total (unconditional) variance

*(S*=

_{i}*V*;

_{i}/V*S*=

_{ij}*V*, etc.), and represent the contribution to the total output variance of single parameters

_{i}/V*(Si)*or combinations of parameters (

*S*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.,

_{ij},*ST*is the sum of all terms in Eq. 9 involving parameter i, divided by

_{i}*V).*Given the rather high number of model parameters (from 13 to 20, Table 1), we calculated first-order

*(S*and total-order

_{i})*(ST*sensitivity indices only. Thus, the difference between

_{i})*ST*and

_{i}*S*, estimates the fraction of total output variance due to interactions (from order 2 up to order

_{i}*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 2^{16}.

#### 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 MOD_{Monod} and MOD_{Monod+Over} were calculated from those imposed for models MOD_{RevMonod} and MOD_{RevMonod+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 *(fC _{YR}*

_{soil}) were set by assuming that

*C*

_{YR}_{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 *Y _{HL} (fN_{YHL})* was calculated in order to keep the C to N ratio of the pool

*Y*equal to 15. This value permitted to explore different partitioning of input C and N among the

_{HL}*Y*,

_{L}*Y*and

_{R}*Y*pools, avoiding that the sum of

_{HL}*fN*and

_{YL}, fN_{YR}*fN*was higher than one, and ensuring that the C to N ratio of the pool

_{YHL}*Y*was close to that of O.

_{HL}#### Variability of model outputs

Besides reporting sensitivity indices, we decided to show also the variation of model outputs over time. For each model output (CO_{2} or SMN, at one of the seven dates) we selected four reference points corresponding to the average, median, 10^{th} and 90^{th} 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 CO_{2} 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 CO_{2} emissions (mg CO_{2}-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 CO_{2}-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.*, NNM_{day = i} = (SMN_{Day = i} – SMN_{Day = 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.

## Results

### Variability of CO_{2} emissions and soil mineral nitrogen

The trend over time and the variability of simulated net CO_{2} 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 CO_{2} (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 CO_{2} emissions raised slightly faster in MOD_{Monod} and MOD_{Monod+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 10^{th} percentile of CO_{2} distributions. In such case, MOD_{First} and MOD_{Monod} decomposed the same amount of manure-C (on average 28%), and slightly more than MOD_{RevMonod} (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 10^{th} percentile points in Figure 2B), or lower (by 0.1-0.8% of manure C) in simulations mineralising N (average, median and 90^{th} 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 90^{th} percentile points, while no differences among models occurred in simulations corresponding to the 10^{th} 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 90^{th} percentile of output distributions (Figure 3A and B). At the end of incubation (Day 189), extra C mineralised in MOD_{First+Over}, MOD_{Monod+Over}, and MOD_{RevMonod+Over}, compared to respective model variants implementing the N inhibition hypothesis (MOD_{First}, MOD_{Monod}, and MOD_{RevMonod}), 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 10^{th} percentile of NNM distribution.

Replacing first-order with Monod decomposition kinetics produced different results over time for both net CO_{2} emissions and NNM. Until Day 14 (Figure 3A), mineralised C in MOD_{First} and MOD_{First+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, CO_{2} emissions occurred at faster rates in MOD_{Monod} and MOD_{Monod+Over} than in other models, and at Day 189, net CO_{2} 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 CO_{2}, 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 90^{th} percentile of CO_{2} distribution (Figure 4A). Conversely, for other points (average, median and 10^{th} 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 10^{th} percentile points after Day 73. In fact, NNM at Day 189 was higher in MOD_{First+Over}, MOD_{Monod+Over}, and MOD_{RevMonod+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 CO_{2} emissions and NNM also in CR. At the end of incubation, in simulations corresponding to average, median and 90^{th} percentile of net CO_{2} 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 10^{th} percentile of net CO_{2} distributions, when more CO_{2} accumulated in both MOD_{Monod} and MOD_{RevMonod} (on average 10% of residue-C) compared to MOD_{First}. 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).

### Sensitivity indices

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 (CO_{2} 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%.

#### Unamended soil

In CON, parameters that mostly affected CO_{2} emissions were those related to the turnover of the young resistant pool *(fC _{YR} soil* and

*k*and the substrate use efficiency

_{YR}soil)*(eff);*these three parameters account together for a relevant fraction of CO

_{2}variance (from 45 to 82%, depending on the model). In addition, in models implementing reverse Monod kinetics, the decomposition constant of microbial biomass (

*k*) and the half-saturation constant

_{B}*(k*accounted for a further 29% of CO

_{rm})_{2}variance, while in those implementing classical Monod kinetics, the half-saturation constant

*(k*alone explained 26% of CO

_{mYR}soil)_{2}variation. The importance of microbial biomass parameters

*eff*and

*k*usually decreased with time, especially after Day 44, while that of the half-saturation constants increased. Conversely, the effect of

_{B}*fC*increased after Day 44 in all models, while the trend of

_{YR}soil*kYR soil*depended on the type of decomposition kinetics:

*S*decreased after Day 44 in MOD

_{i}_{First}and MOD

_{First+Over}, while it remained constant until Day 189 in all other models.

The three microbial biomass parameters eff, *fN _{B}* and, with the exception of MOD

_{Monod}and MOD

_{Monod+over},

*k*, explained a relevant fraction also of the variation of SMN (from 56 to 69% depending on the model). As observed for CO

_{B}_{2}, first-order indices of these parameters usually tended to decrease after Day 44.

#### Amended soil

As expected, in the amended treatments parameters related to the turnover of pools *C _{YR} soil* and

*O*had no appreciable effect in any model (first-order indices always less than 10%) on both net CO

_{2}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 CO

_{2}and net SMN variances. With few exceptions, average first-order indices of the labile fraction of input C (

*fC*) were high for both net CO

_{YL}_{2}(5-25%) and net SMN (9-23%) in all models. Moreover, the importance of

*fC*usually diminished as time passed, and as the C to N ratio of the input raised (especially from SM to CR).

_{YL}Other input parameters that had an effect on net CO_{2} variation were the decomposition constants of *Y _{L}* and

*Y*pools

_{R}*(k*and

_{YL}*k*that accounted together for 9-40% of net CO

_{YR})_{2}variance. Similarly to

*fC*, the effect of first-order indices of

_{YL}*k*decreased with increasing the C to N ratio of the input; the opposite happened for

_{YL}*kYR*, especially for manures and CR, with few exceptions. Also the trend over time of

*S*differed between

_{i}*k*and

_{YL}*k*the former had a decreasing effect, while the latter had an increasing effect over time. Differently from what was observed for net CO

_{YR}:_{2}, sensitivity of net SMN to

*kYL*was relevant in SM (for models implementing classical and reverse Monod kinetics), while sensitivity to

*k*was relevant in CR (for all models). The effect of

_{YR}*k*usually increased until Day 44 and then stabilised thereafter.

_{YR}In models with Monod kinetics, also half-saturation constants were important *(k _{rm}*,

*k*and, with the exception of SM,

_{mYL}*k*These parameters gave an average contribution of 8-25% to net CO

_{mYR})._{2}variance and, solely for CR, of 2-8% to net SMN variance. For both output variables, sensitivity indices of

*k*and

_{rm}*k*decreased until Day 44, and then remained constant until Day 189, while those of

_{mYL}*k*usually increased during time.

_{mYR}#### Parameter interactions

The sum of first-order sensitivity indices (Σ*S _{i}*) accounted for 70-94% and 57-88% of total output variance for CO

_{2}and SMN, respectively (Table 5), showing that interactions among model parameters (1 – Σ

*S*) 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).

_{i}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 CO_{2} variance were for parameters *fC _{YR} soil, k_{YR} soil* and the half-saturation constant

*km*Interactions of microbial biomass parameters

_{YR}soil.*eff*and

*fN*also markedly contributed to SMN variance. Besides

_{B}*eff*and

*fN*, in the amended treatments net CO

_{B}_{2}and net SMN were also sensitive to interactions of input parameters with high

*S*, and

_{i}(fC_{YL}, k_{YL}, k_{YR}, k_{rm}, k_{mYL}*k*in many model × input combinations.

_{mYR})## Discussion

In this section, we first compare the simulated dynamics of CO_{2} 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 CO_{2} emissions and net nitrogen mineralisation

Dynamics of simulated net CO_{2} 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 CO_{2} emissions, and was accompanied by the complete depletion of the more labile input organic matter in 14-73 days (*Y _{L}*); thereafter, decomposition of the more resistant input organic matter (

*Y*) mostly contributed total CO

_{R}_{2}emission (Figures 2A, 3A and 4A). During the whole incubation, only a fraction of initial

*Y*(36-70%, 42-67%, and 33-61% in LM, SM, and CR, respectively) was mineralised (data not shown). Conversely mineralisation of

_{R}*Y*was insignificant (1-3% of input C) for all inputs.

_{HL}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 10^{th} 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 CO_{2} 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 CO_{2} variance. The high sensitivity of CO_{2} to *eff* was due to its direct effect on C turnover (Figure 1), regulating the amount of C respired as CO_{2}, and its feedback on C decomposition by defining N required by microbial biomass (Eq. 4). Thus, when shortage of N occurs, CO_{2} 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 *fN _{B}* (and thus to the C to N ratio of soil microbial biomass). Together,

*eff*and

*f*explained from 26 to 63% of SMN variance; this is justified by their effect in determining N requirements by microbes (Eq. 4).

_{NB}In the amended treatments, neither the most recalcitrant fraction of input C *(fC _{YHL})* nor its decomposition constant

*(k*gave effects on net CO

_{YHL})_{2}emissions and SMN in all models. This can be attributed to range that we have chosen for

*k*(Table 1) and to the reduced simulation time span that made the

_{YHL}*Y*pool almost inert (data not shown). Conversely, C decomposition and SMN were very sensitive to variations of the size

_{HL}*(fC*and decomposition constant

_{YL})*(k*of the labile young pool

_{YL})*Y*(Tables 5 and 6). Both parameters influence the mineralisation of labile C; moreover,

_{L}*fC*is fundamental in regulating input-N availability by switching from N-rich to N-poor

_{YL}*Y*pools (Tables 1 and 2), and thus strongly impacts SMN dynamics.

_{L}### Effects of carbon-overflow on CO_{2} 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 *k _{P}.* 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 MOD

_{First+over}and MOD

_{Monod+Over}(15-16% and 18-23% of inputC in SM and CR, respectively) and was lower for MOD

_{RevMonod+Over}(11% and 9% of input-C in SM and CR, respectively). Results of MODMOD

_{First+over}and MOD

_{Monod+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 CO_{2} emissions to the parameter describing polysaccharides decomposition rate *(k _{P})* 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,

*fC*and

_{YL}*kYL)*masked that of

*k*in determining CO

_{p}_{2}variation. Different results would have been obtained if all C in excess was respired as CO

_{2}, 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 CO

_{2}production) were proved to occur (Russel and Cook, 1995).

### Effects of decomposition kinetics on CO_{2} 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 *k _{mYR} soil, k_{mYL}*, and

*k*(MOD

_{mYR}_{Monod}and MOD

_{Monod+Over}), and

*krm*(MOD

_{RevMonod}and MOD

_{RevMonod+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 CO_{2} 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 CO_{2} emissions were very low (10^{th} 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 CO_{2} emissions in MOD_{Monod} and MOD_{Monod+Over} was due to the growth of soil microbial biomass that in turn enhanced C decomposition rates *(k’ × C _{R}*, Eq. 2), while the effect of substrate exhaustion

*(C*, 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).

_{s}/(k_{m}+C_{s})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 MOD_{First} that showed negative net CO_{2} emissions after Day 29. This result was an artefact of model parametrization. In these simulations native soil pools (Y_{R} *soil* and *O)* had a C to N ratio higher than *C/N _{crit}* (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 CO

_{2}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 CO_{2} 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 *Y _{R} soil* was higher than

*C/N*and SM and CR exhausted the available mineral N causing N limitation, the models implementing the N inhibition hypothesis decomposed less

_{crit}*Y*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 MOD

_{R}soil_{First}(–10% and – 21%). In addition, models implementing Monod kinetics (irrespective of N control mechanism on C decomposition), always showed positive PEs on

*Y*when decomposition of native SOM was not N-limited. This is represented in Figure 6 (only for models with C-overflow) with negative net

_{R}soil*C*, less

_{YR}soil, i.e.*C*at the end of the incubation in the amended treatments than in CON.

_{YR}soilIncreasing the amount of added-C, more *Y _{R} soil* was decomposed (Figure 6A), in particular for LM and SM. Decomposition of

*Y*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

_{R}soil*Y*was mineralised. Thus, in this case, the input C to organic N ratio is a proxy of input decomposability.

_{R}soilMoreover, distribution statistics of net *Y _{R} 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

*Y*, while low, often N-limited, microbial growth did not stimulate

_{R}soil*Y*decomposition, and thus did not produce appreciable PEs.

_{R}soilIn addition, higher C decomposition in models implementing classical Monod than reverse Monod kinetics (Figures 2A, 3A and 4A) gave higher PEs (Figure 6).

Models implementing Monod kinetics predicted PE according to cometabolism of fresh (Y_{L}, Y_{R}) and native (Y_{R} *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/N _{crit})* 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.

## Conclusions

Model modifications (introduction of C-overflow and Monod decomposition kinetics) had an effect on simulated CO_{2} 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 CO_{2} 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.