Evaluation and optimisation of the soil carbon turnover routine in the MONICA model (version 3.3.1)

Simulation models are tools commonly used to predict changes in soil carbon stocks. Prior validation is essential, however, for determining the reliability and applicability of model results. In this study, the process-based biogeochemical model MONICA (Model of Nitrogen and Carbon dynamics on Agro-ecosystems) was evaluated with respect to soil organic carbon (SOC), using long-term monitoring data from 46 German agricultural sites. A revision and parameterisation of equations, encompassing crop- and fertiliser-specific C contents and the abiotic factors of soil temperature, soil water and clay content, were undertaken and included in the model. The modified version was also used for a Morris elementary effects screening method, which confirmed the importance of environmental and management factors to the model's performance. The model was then calibrated by means of Bayesian inference, using the Metropolis–Hastings algorithm. The performance of the MONICA model was compared with that of five established carbon turnover models (CCB, CENTURY, C-TOOL, ICBM and RothC). The original MONICA model systematically overestimated SOC decomposition rates and produced on average a ∼17 % greater mean absolute error (MAE) than the other models. The modification and calibration significantly improved its performance, reducing the MAE by ∼30 %. Consequently, MONICA outperformed CENTURY, CCB and C-TOOL, and produced results comparable with ICBM and RothC. Use of the modified model allowed mostly adequate reproduction of site-specific SOC stocks, while the availability of a nitrogen, plant growth and water submodel enhanced its applicability when compared with models that only describe carbon dynamics.

1 Introduction

The historic conversion from native to agriculturally managed soils released 116 Pg of soil organic carbon (SOC) as CO₂ to the atmosphere (Sanderman et al., 2018). Current projections estimate that unless no measures are taken, the decline in SOC stocks will continue, further intensifying the release of CO₂ from soils (Riggers et al., 2021; Zhao et al., 2021). However, the sequestration of atmospheric CO₂ in the soil is a climate change mitigation strategy that has substantial potential (Poeplau and Don, 2015; Fuss et al., 2018; Amelung et al., 2020). In particular, croplands are promising landscapes, since they are depleted in SOC and generally accessible for cultivation practices (Poeplau and Don, 2015; Sanderman et al., 2018). Researchers have proposed various agronomic systems, including agroforestry, reduced tilling and perennial and cover cropping, to optimise agricultural opportunities for sequestering carbon and thus partly restore the historically lost SOC (Smith et al., 2005). Estimates of the global capacity of agricultural soils for sequestering carbon range from 2 to 5 Pg C a⁻¹ (Fuss et al., 2018). However, identifying and promoting site-specific measures in a highly diverse agricultural landscape encompassing more than 570 million farms globally and affecting ∼36 % of the Earth's terrestrial surface is a challenge (Lowder et al., 2016; Amundson and Biardeau, 2018). Despite the uncertainty in their practicability, international initiatives, such as the “4 per 1000” initiative, and scientific working groups advise the rapid modification of agricultural systems so that they act as carbon sinks (Minasny et al., 2017; Amelung et al., 2020). In contrast, some scientists focus on the potential for underestimated detrimental effects and advocate a well-thought-out approach to sequestering SOC on a large scale (Sommer and Bossio, 2014; Lugato et al., 2014; Amundson and Biardeau, 2018; Murphy, 2020).

Process-based models are tools that can be used to assess mitigation capabilities of agricultural systems for different environmental and management conditions and to analyse disagreements and accuracy of alternative estimation approaches. They are designed to represent the behaviour of natural systems, based on a selection of mathematical formulations of generally accepted biophysical principles. Carbon turnover models predict the dynamics of SOC in response to management and climate. Agroecosystem models combine carbon and nitrogen turnover routines with algorithms that represent soil water dynamics, greenhouse gas emissions and crop growth and allow the simultaneous prediction of multiple variables from different domains of the soil–plant–atmosphere continuum. Due to the explicit simulation of crop growth in response to weather, site conditions and crop management, agroecosystem models exhibit a clear advantage over simple soil carbon turnover models, as they can translate any changes in conditions that affect biomass input to soil, e.g. climate or crop management, directly to responses in SOC. Potentially, they could also simulate any feedback from those changing SOC dynamics to crop growth, for example, via nitrogen release or immobilisation. The disadvantage of this is that complex models require an accurate reproduction of the carbon input variability, whereas simpler carbon models estimate the carbon inputs from measured yields. Processes that are not described by the model, such as the effects of phosphorus and other nutrients and the effects of plant diseases on plant growth, could further constrain the ability of complex models to simulate adequate carbon inputs. Additionally, complex agroecosystem models are inherently black boxes that make it difficult to understand the various factors and feedback loops. They have as of yet rarely been tested for their ability to do exactly this, and they have just as rarely been compared to simple soil carbon turnover models to identify whether their advantage in applicability comes with a reduced predictive power (Smith et al., 1997; Farina et al., 2021). In any case, to ensure the adequate predictive capability of such models, complex or simple, a comprehensive validation process using measured data with different soil types, management practices and weather conditions is required.

MONICA (Model of Nitrogen and Carbon dynamics on Agro-ecosystems; Nendel et al., 2011, 2014) is a process-based agroecosystem model that simulates the effects of spatiotemporal environmental variability and agricultural practices on crop yields and soil condition of mineral soils (Nendel et al., 2014). It is a successor of the HERMES model, which was initially developed to estimate N dynamics in the soil–crop system (Kersebaum, 1995, 2006). A comprehensive SOC turnover submodule based on the Daisy model was added, replacing the simpler approach in the HERMES model (Hansen et al., 1991; Abrahamsen and Hansen, 2000). The holistic principle of the MONICA model supports research questions on the interlinkage between site conditions, plant development and C and N dynamics in a changing climate. The performance of the MONICA model showed the capability to adequately reproduce crop yields for different crops (Asseng et al., 2013; Bassu et al., 2014; Kothari et al., 2022; Rötter et al., 2012; Salo et al., 2016) and more complex crop rotations (Kollas et al., 2015; Kostková et al., 2021). For SOC dynamics under bare fallow treatment, MONICA performed similarly to simple C turnover models (Farina et al., 2021) but furthermore demonstrated also good performance when simulating short-term high-resolution CO₂ exchange in a soil–plant system (Specka et al., 2016). The Daisy model, another agroecosystem model, had already been tested against long-term soil carbon experiments (Smith et al., 1997) and performed comparably well to soil carbon turnover models, such as RothC (Coleman and Jenkinson, 1996) or CANDY (Franko et al., 1995).

The question now is whether the MONICA model with its carbon turnover submodule is able to adequately simulate changes in SOC stocks when considering different agronomic practices. It is also important to establish how the performance of a complex model like MONICA compares with simple C turnover models such as RothC, ICBM, C-TOOL, CENTURY and CCB. To answer these questions, we have the following four objectives: (i) to validate the SOC submodule of the MONICA model with data from 46 agricultural long-term monitoring sites in Germany; (ii) if necessary, to enhance its applicability in view of the diversification of the dataset used for validation; (iii) to identify the main drivers of C sequestration and decomposition in the model; and (iv), if necessary, to improve model performance.

2 Materials and methods

2.1 Data composition and site characteristics

For calibration and validation of the SOC turnover subroutine of the MONICA model, long-term monitoring data collected by the State Office for Mining, Energy and Geology (Landesamt für Bergbau, Energie und Geologie, LBEG) were used. In 1990, the state government of Lower Saxony (Germany) established a programme to permanently observe agricultural and forest soils in the federal territory. Its goal was to monitor soil changes and identify adverse ecologic and economic implications and thus predict the risk of imminent damage for the sustainable utilisation of soils.

The examined area is part of the North German Plain, with highlands up to 971 m in its southern region. It is located at latitudes between 51 and 54^∘ N and longitudes between 7 and 11^∘ E (Fig. 1a). It is characterised by a maritime climate, with mean annual temperatures ranging from 8.9 ^∘C to 10.4 ^∘C. Average annual precipitation ranges from 574 to 887 mm, with a positive climatic water balance in spring, autumn and winter and a negative balance in summer.

The database consists of 48 agricultural sites, of which 42 are conventional agricultural systems and 6 are organic. Two sites were removed due to their observation periods being shorter than 10 years, leaving 46 entries for analysis. The observation data consist of 33 sites that have been used for agriculture for more than a century. The remaining 13 sites were converted from forests, meadows or peatlands fewer than 5 decades before the start of the observation period. Soil sampling was generally conducted at regular intervals of between 1 and 5 years (Höper and Meesenburg, 2021). They were sampled at a depth of up to ∼30 cm and consisted of four replicates, which were combined as mixed samples.

In addition, soil horizons were sampled every 10 years at a minimum depth of 1 m. During the soil monitoring programme, a total of 333 soil carbon measurements (replicates not included) were taken at the selected sites from 1992 until the end of 2015. Between the sites, the mineral composition was diverse, with soil textures predominantly sand (20 sites) or silt (21 sites) and less frequently clay soils (five sites; Fig. 1b). The key properties of the examined soils were highly variable between the sites, with SOC ranging from 0.84 % to 4.91 %, $\mathrm{C}/\mathrm{N}$ from 5.25 to 23.66, bulk density from 1.14 to 1.76 g cm³ and pH from 3.9 to 7.9. A special distinctiveness of the investigated area is the occurrence of so-called black sands, which are characterised by comparably high, stable C contents in sand-dominated locations (Vos et al., 2018). Site-specific details are summarised in the Appendix (Table A1).

Figure 1 Location (a) and soil classification (b) of the 46 agricultural long-term monitoring sites in Lower Saxony, Germany.

The soil and crop management at the sites was undertaken by the landowners or their associates. Therefore, there is a high degree of variability in the agronomic practices performed at the sites, such as tillage depth and type, fertiliser application, irrigation, crop rotation, harvest and crop residue management. The reported management reflects common agricultural practices for the site-specific regions (more details are given by Höper and Meesenburg, 2021). Control treatments were not available for the sites. Simulated soil temperatures and moisture contents were evaluated using data from 11 German long-term field experiments (Table A2; Flessa et al., 1995; Frolking et al., 1998; Schmädeke, 1998; Leidel et al., 2000; Wolf et al., 2014; Flessa et al., 2017; Mallast et al., 2021; Winkhart et al., 2022).

2.2 Model description

The MONICA model is comprised of subroutines that characterise plant growth, C and N dynamics, soil water and soil temperature (for a more detailed description, see Nendel, 2014). Execution of the programme requires comprehensive input data, including weather, soil and management information, as well as parameter values describing model settings. Every model iteration simulates soil and plant conditions in a daily time step for soil depths up to 2 m. Mathematically, the model is one-dimensional (point model), and the spatial representativeness of the model output depends primarily on the calibration data. Generally, the spatial representation ranges between 1 m² and 1 km². The soil profile is partitioned into 20 layers, each with 10 cm thickness of homogenous soil conditions. Soil input data are uniformly extrapolated to the relevant depth if the entered soil horizons are shorter than 2 m.

According to the Daisy model, the SOC content is partitioned into three groups composed of dead native soil organic matter (SOM), living soil microbial biomass (SMB) and added organic matter (AOM), each of which contains one rapidly (SOMf, SMBf and AOMf) and one slowly decomposing subpool (Fig. 2; SOMs, SMBs and AOMs; Hansen et al., 1991; Bruun et al., 2003). These compartments are calculated separately for each soil layer.

MONICA is usually initialised using a standard distribution of the SOC pools, which should correspond to an equilibrium at the beginning of the simulation (Bruun and Jensen, 2002). During the initialisation, inert SOC (IOM) is subtracted from the total SOC, based on equations adopted from Falloon et al. (1998). Then the SOM and SMB pools are determined by multiplying the manually entered SOC contents by the corresponding distribution parameters. The AOM pools are initialised and altered by added C contents from crop residues and organic fertilisers. Organic fertiliser and crop types are parameterised with specific N contents (NH${}_{\mathrm{4}}^{+}$ and NO${}_{\mathrm{3}}^{-}$), decomposition rate coefficients and pool distributions. Furthermore, plant biomass accumulation is divided into five groups, namely root, leaf, shoot, fruit and sugar. For perennial plants, an additional permanent structure is considered. Each group has its own parameters controlling the growth rate and the biomass that can potentially be harvested or added to the AOM pools. Following a harvest command, accrued crop residues left on an agricultural field are immediately distributed to the AOM pools. Root remains are fragmented into the soil layers, depending on the root mass distribution in the soil profile. During tillage, the SOM content is distributed homogenously to the soil layers affected by tillage depth.

The initialised pools are managed by first-order degradation due to decomposition, with each controlled by related parameters. This process follows a decay rate that is proportional to the pool size (Hansen et al., 1991). In addition, the decomposition rate of organic matter is affected by soil temperature, soil water and clay contents (Bruun et al., 2003; Farina et al., 2021). These environmental conditions influence each pool as rate-modifying factors, with the exception of clay, which solely affects SMBs. Decomposed AOMf and AOMs are subsequently transferred to the SMB pools, from where fractions are respired to the atmosphere as CO₂ or alternatively converted to SOMs and SOMf, following a pool-specific carbon use efficiency coefficient. Biologically degraded SOMs and SOMf goes through a similar process and are either respired or returned to the SMB pools. A detailed description of the SOM turnover subroutine is depicted in Hansen et al. (1991) and Nendel (2014).

Figure 2 Conceptual SOC pool distribution in MONICA, subdivided into groups of added organic matter (AOM), soil microbial biomass (SMB), soil organic matter (SOM) and inert organic matter (IOM). Except for IOM, each group consists of a slow (s) and a fast (f) pool.

2.3 Model modification

In total, 16 crop plants (amaranth, asparagus, buckwheat, carrot, flaxseed, field bean, hemp, lupine, serradella, spelt, strawberry, summer onion, sunflower, turnip and vetch) and 8 types of organic fertilisers (biogas slurry, calf slurry, dried poultry manure, hoof and horn meal, lime-stabilised sewage sludge, sheep manure and vinasse) were added to the model. This enabled the reproduction of all occurring land management characteristics in the dataset. Parameter values for crop growth (Tesar, 1984; Jeuffroy and Ney, 1997; Allen et al., 1998), fertiliser nutrient contents (Möller and Schultheiß, 2015; Wiesler et al., 2016), crop residues and organic fertiliser decomposition and partitioning coefficients (Gilmour et al., 1998; Thuriès et al., 2001; Lashermes et al., 2009; Peltre et al., 2012; Semenov et al., 2019) were taken from the literature. Furthermore, a specific C content parameter (CorgContent) was defined for each crop and organic fertiliser, providing a more adequate replication of C inputs in comparison with the original constant (Möller and Schultheiß, 2015; Wiesler et al., 2016; Ma et al., 2018). In addition, an overhaul of existing crop and fertiliser types was undertaken, ensuring distinct and more reasonable values for essential parameters. Crop growth was calibrated according to the measured yield data, while the residue-to-crop product ratio was adapted from the literature and an ensemble of allometric functions (Bolinder et al., 2007; Franko et al., 2011; Jacobs et al., 2020; Taghizadeh-Toosi, 2020). The differences between the allometric functions served as indicators of uncertainty in the estimation of C inputs from crop residues.

The functions that describe the influence of soil temperature, water and clay content on decomposition were revised and parameterised. This enabled all relevant carbon turnover factors to be tested for uncertainty and subsequently calibrated. The optimal values of the parameters were used as the evidential basis for the prior probability distributions, which is necessary for model calibration by means of Bayesian inference. For further explanations, see Sect. 2.5 and Table 1.

Extensive research provided evidence of the influence of temperature on physical soil processes, microbial metabolisms and enzymatic activities (Pietikäinen et al., 2005; Haddix et al., 2006; Meyer et al., 2018). The novel temperature function is described by the following function Eq. (1):

$$\begin{array}{}\text{(1)}& f_{\mathrm{TOD}}={\displaystyle \frac{\mathrm{1}}{(\mathrm{1}+e^{T_{\mathrm{d}}-\left(e+T_{\mathrm{opt}}\right)}{)}^{Q_{\mathrm{10}}/\left(T_{\mathrm{opt}}/\mathit{\pi }\right)}}}\left(-\mathrm{1}+Q_{\mathrm{10}}^{\frac{T_{\mathrm{d}}}{\mathrm{15.76}}}\right),\end{array}$$

where T is the daily soil temperature in ^∘C, Q₁₀ is the Q₁₀ temperature coefficient specified by the parameter QTenFactor, and T_opt is the optimal temperature for decomposition in ^∘C defined by the parameter TempDecOptimal (Fig. B1). From a literature analysis, a range of 38.2±7.5 ^∘C was identified for the optimal temperature (Pietikäinen et al., 2005; Richardson et al., 2012; Menichetti et al., 2015; Liu et al., 2018; Čapek et al., 2019) and 2.6±0.6 for the Q₁₀ temperature coefficient (Haddix et al., 2006; Vanhala et al., 2007; Conant et al., 2008; Gillabel et al., 2010; Meyer et al., 2018).

Microbial aerobic respiration is often driven by changing soil water contents up to a point where water becomes a restricting property for the availability of oxygen (Curiel Yuste et al., 2007; Schimel et al., 2007). Based on an analysis of incubation and field experiments, a Gaussian function was proposed, described by the following Eq. (2):

$$\begin{array}{}\text{(2)}& f_{\mathrm{MOD}}=e^{-\mathrm{18}({\text{WFPS}}_{\mathrm{d}}-M_{\mathrm{opt}}{)}^{\mathrm{2}}},\end{array}$$

where WFPS_d is the daily water-filled pore space in percent, and M_opt is the optimal water-filled pore space (WFPS) for decomposition in percent defined by the parameter MoistureDecOptimal (Fig. B2). Optimal WFPS for respiration was determined to be around 59.5±9.9 % (Linn and Doran, 1984; Doran et al., 1990; Liebig et al., 1995; Aon et al., 2001; Gabriel and Kellman 2014; Zhang et al., 2015).

Several studies have indicated that increasing soil clay contents have a negative effect on respiration due to the reinforced inaccessibility of SOM for microorganisms (Six and Paustian, 2014; Churchman et al., 2020; Fomina and Skorochod, 2020). In the MONICA model, the modulating effect of clay on the decomposition rates of SMBs is by default adjustable through the parameter LimitClayEffect. An inverse logistic function was proposed, described by the following Eq. (3):

$$\begin{array}{}\text{(3)}& f_{\mathrm{COD}}={\displaystyle \frac{(\mathrm{1}-\text{LCE})}{(\mathrm{1}+e^{-\mathit{\pi }+C\cdot \mathrm{16}})+\text{LCE}}},\end{array}$$

where C is the site-specific clay content (in kg kg⁻¹), and LCE is the limit at which the clay content affects the mineralisation defined by the parameter LimitClayEffect (Fig. B3). Laboratory evidence suggests that soils amended with clay minerals reduce mineralisation rates by $\sim \mathrm{1.5}\pm \mathrm{1.1}$ % for every percent of clay added to the soil, followed by a gradual reduction in the effectiveness after a certain clay content is exceeded (Nguye and Marschner, 2014; Schapel et al., 2018; Riaz and Marschner, 2020; Liddle et al., 2020).

2.4 Screening of variables relevant to soil carbon turnover

For the purpose of understanding the relationship between input variables, parameters and the predictive proficiency of the MONICA model, a comprehensive sensitivity analysis was conducted using an improved Morris elementary effects screening method (Fig. C1; Morris, 1991; Campolongo et al., 2007; Pujol, 2009). The algorithm demonstrates the capability to determine the importance of parameters in relatively short computing times (Campolongo et al., 2007; Confalonieri et al., 2010). The effect of 51 variables was examined, of which 15 were environmental, 6 were related to the crop, 4 to the management and 26 to SOC turnover. Of the turnover parameters, six were unique for each fertiliser and crop residue type.

Depending on the tangible properties of the variables (factual or conceptual), the upper and lower bounds of each factor were either derived from the literature or intuitively inferred from model characteristics. After building the parameter space, all the factors were scaled to a range from 0–1. During each iteration, one variable was changed to a random value in the defined space (levels = 5; grid jump = 3), while the remaining variables were unaltered. Subsequently, the effect of the modified factor on SOC stocks was assessed by comparing the present iteration with the previous one. Based on the magnitude of the changes and the degree of the interaction with other factors, a parameter importance ranking was identified using the two measures μ^∗ and σ (Campolongo et al., 2007).

Common recommendations for the number of model iterations are between 10 to 50 for each variable considered (Campolongo et al., 2007). However, several studies suggest a much higher number of iterations to achieve robust results (Sarrazin et al., 2016; Vanuytrecht et al., 2014; Pianosi et al., 2016). As a base value, 10³ iterations and three runs were used in total to accomplish adequate coverage of possible input spaces and to estimate the uncertainty in the elementary effects.

The programming language R (version 4.1.0) with the graphical user interface RStudio (version 1.4.1717) and the sensitivity package (Iooss et al., 2020; version 1.26.0) were used for screening of variables relevant to soil carbon turnover in the MONICA model.

2.5 Parameter calibration

Bayesian inference techniques are an efficient approach for the calibration of process-based models (Vrugt, 2016; Van de Schoot et al., 2021). Common methods for prior elicitation, likelihood specification and posterior distributions have been extensively researched (Van de Schoot et al., 2021). For instance, Bayesian calibration has been applied successfully to forest biomass growth (Van Oijen et al., 2005; Svensson et al., 2008) and agroecosystem models (Lehuger et al., 2009; Gurung et al., 2020).

As a first step, the experimental data were divided into one subset for calibration and one for validation. In the process, the data were clustered into three groups based on their initial SOC stocks in order to achieve a general representation of the variability between the sites. From each of the groups, two-thirds of the sites were randomly sampled for model calibration (totalling 31 sites) and one-third for model validation (totalling 15 sites). The parameter selection process for the calibration was adapted from the results of the sensitivity analysis and the recommendations from Minunno et al. (2013), excluding local and 25 % of the least important parameters.

The calibration of the 14 selected parameters was undertaken by means of Bayesian statistical inference. In Bayesian calibration, the determination of possible parameter values is based on evidential probabilities (Van de Schoot et al., 2021). Accordingly, it was important that the upper and lower limits of the parameters were based on observed, realistic values. Therefore, a comprehensive literature analysis combined with expert knowledge provided the basis for the prior probability distributions, as shown in Sect. 2.3 (Oakley and O'Hagan, 2007). A weakly informative prior probability distribution was constructed, with the mean of all identified measurements as the most likely value and the observed lower and upper bounds as the variance for each parameter (Table 1; Van de Schoot et al., 2021). At the beginning of the calibration, a random initial value (within the prior probability distribution) was generated for each parameter, creating a parameter vector for each iteration. Selection for subsequent parameter values was limited by the minimum and maximum value of the prior probability distribution and a randomly determined step size corresponding to 1 %–25 % of the parameter space (Van Oijen et al., 2005). Accordingly, the subsequent parameter value changes from the previous value by a minimum of 1 % and a maximum of 25 % of the entire parameter space. During each iteration, MONICA was executed using the created parameter vector. Afterwards, the generated output was compared with the measured data using a log–Laplace likelihood function, suggested by Vrugt (2016), which estimated the association between simulated and observed values. This association was a representation of the posterior probability of each parameter vector. The Laplace distribution showed robustness against outliers and random variations, resulting in less biased and more consistent parameter estimates (Schoups and Vrugt, 2010). The Gelman–Rubin diagnostic (Gelman and Rubin, 1992) was used as a measure for the convergence of the calibration and was conducted with the R package coda (Plummer et al., 2006).

Table 1 Description of the 14 parameters chosen for the calibration of the C submodel in MONICA. The prior distribution depicts the upper and lower bounds of possible parameter values inferred from the literature; the posterior distribution describes the median value of the posterior and the corresponding standard deviation.

2.6 Quantitative methods and model evaluation specifics

The statistical methods for evaluating the sampling distributions and model performances were conducted with the R packages hydroGOF (Zambrano-Bigiarini, 2020; version 0.4-0) and car (Fox et al., 2019; version 3.0-12). Shapiro–Wilk (Shapiro and Wilk, 1965) and Levene's tests (Levene, 1960) were used on the measured fluxes to test the population for normality and variance equality. The Shapiro–Wilk test resulted in rejection of the null hypothesis for the variables (p<0.001), inferring that the examined datasets were not normally distributed. Levene's test resulted in acceptance of the null hypothesis, implying that the variance in the samples was homoscedastic (p=0.88). Using the interquartile range criterion (Jargowsky and Yang, 2005), 19 outliers were identified in the data, which represents ∼6 % of the cases. Q−Q plots confirmed the assumption that the data points were nonlinear distributed.

To improve representativeness, the years 1992 to 2000 were excluded from the analysis, reducing the observed number of cases from 333 to 223. The reason for this approach was justified by the inherent risk of statistical biases from overfitting of measured and simulated values due to the initialisation of simulated SOC stocks based on measured values, as significant changes in SOC stocks usually take several years. Furthermore, a spin-up period of 20 simulated years was included to ensure a realistic reflection of available fresh matter in the soil, since the AOM pools were initialised equal to zero. The chosen period corresponds to the maximum duration of mineralisation of most AOM (Thuriès et al., 2001; Semenov et al., 2019).

Based on the reduced sample size and the abovementioned statistical assumptions, Kendall's rank correlation coefficient (T; Kendall, 1948) was selected to assess the statistical dependence between the observed and simulated values. Kendall's τ is the preferred method for measuring the ordinal association in data consisting of small sample sizes and outliers and is more robust and efficient than the Pearson correlation coefficient and Spearman's rank correlation coefficient (Bonett and Wright, 2000). In addition, the model results were evaluated with statistical methods like the mean absolute error (MAE; Willmott and Matsuura, 2005) for an analysis of the average error between observed and simulated pairs and with the Nash–Sutcliffe model efficiency coefficient (NSE; Nash and Sutcliffe, 1970) to assess the predictive skill of the simulation model (Fig. E1). Furthermore, the slope intercept was considered in order to better assess the behaviour of the models.

3 Results

3.1 Simulation of soil and plant growth variables

Soil temperatures were well predicted, with an average error of just under ∼1.8 ^∘C (Table 2). Low temperatures were slightly overestimated, and high temperatures were underestimated. In comparison, soil water contents were predicted somewhat less favourably, with an average error of ∼5.4 %. Nonetheless, the predictive power was better than the mean of the observations. The model tended to slightly overestimate moderate water contents and underestimate the measured variance. A comparison of the default with the calibrated plant growth parameters was not performed, since new crops were added to the model in order to simulate all possible systems.

The calibrated plant growth model achieved a good representation of measured yields, with an average error smaller than ∼1.6 Mg DM ha⁻¹ for all crops. Overall, there was a linear relationship between the calibrated and measured yields. However, there was a high variability in the estimation quality between different crop types. Some crops, such as triticale and potatoes, were modelled less adequately than silage maize and sugar beet, which was partly explained by the differences in the number of observations and the variance in yields. The C inputs were compared with the estimated mean of the allometric functions. The model produced an average error that was similar to the calculated difference between the allometric functions (∼1 Mg DM ha⁻¹).

Table 2 Statistics describing the performance of the MONICA model regarding soil and plant growth variables. The statistics depicted are Kendall's rank correlation coefficient (T), Nash–Sutcliffe model efficiency coefficient (NSE), mean absolute error (MAE), the slope and intercept for characterising the relationship between measured values (except values for C inputs, which were calculated with allometric functions) and simulated values.

${}^{***}$ p<0.001.

3.2 Quantification of uncertainty and parameter calibration

A strong relationship was identified between SOC stock change rates and C inputs through organic fertilisation and crop residues (Table F1; OrganicFertilization, AOM_DryMatterContent and CorgContent), given the upper and lower bounds of possible values. Additionally, these factors were less affected by other variables, noticeably by the lower σ on average in comparison with μ^*. The climate variables of temperature (Temperature) and precipitation (Daily_Precipitation) were also of greater importance to the SOC turnover rate but had a higher standard deviation of elementary effects. Except for solar irradiance (Global_Radiation), temperature and precipitation, all other climatic factors had an insignificant effect on SOC stocks. Of all the variables that were initialised at the beginning and remained static during the model run, the initial SOC stock (SoilOrganicMatter) was the most influential, followed by the site-specific field capacity (FieldCapacity), soil bulk density (SoilBulkDensity) and pore volume (PoreVolume). Other site-specific variables (Sand, Clay, PermanentWiltingPoint, pH, and Sceleton) were comparably less important. The management factors of tillage (Tillage), mineral fertilisation (MineralFertilization), irrigation (Irrigation) and N content in the added organic matter (CN_Ratio_AOM_Slow, CN_Ratio_AOM_Fast, CN, AOM_NH4Content and AOM_NO3Content) had a small effect on the SOC turnover rate. σ values were generally higher than the overall importance of the input factors, indicating that there is probably an interaction and non-linear effect between the environmental and management variables. In comparison, the sensitivity of the turnover variables was less pronounced and revealed three particularly important parameters, namely one C pool parameter (SOM_SlowDecCoeffStandard) and two weather-regulated decomposition factors (QTenFactor and MoistureDecOptimal). They were similar in significance to precipitation, initial SOC and field capacity. Changes in the parameters QTenFactor and MoistureDecOptimal were more important than in LimitClayEffect, which was in accordance with the influence of clay being limited to the SMBs pool. The effect of the factor describing the optimal temperature for decomposition (TempDecOptimal) was negligible but had a considerable increase in importance when a higher range of possible values for the variable temperature was included. Crop- and fertiliser-specific turnover parameters were also sensitive and coincided with the importance of C inputs for SOC stock change rates. Less pronounced sensitivities involved the parameters describing the microbial dynamic (besides the parameter SMB_UtilizationEfficiency), distinguishable by the SMB in the name. Non-linear effects and interactions were more strongly pronounced in the SOC turnover parameters.

Figure 3 shows the results of the calibration process, with the posterior distribution curve representing the frequency of accepted values and the vertical red line indicating the optimal value as the median of the posterior. Bayesian inference reduced the uncertainty for all parameters, converging the posterior distribution to a normal or lognormal distribution. Throughout the random sampling, the entire plausible parameter space was covered and curve progressions showed that the optimal values were probably within the specified limits. During the calibration process, the posterior distribution of most parameters became narrower in comparison with the prior, confirming the convergence to the optimal values. Similar to Lehuger et al. (2009), the acceptance rate for the proposed parameter values was ∼25 %. Based on the results of the Gelman–Rubin diagnostics, adequate convergence was achieved for all parameters (Fig. G1), signifying the unbiased outcome of the Markov chain Monte Carlo (MCMC). The correlations across parameters remained below a coefficient of 0.2, with the exception of QTenFactor and SOM_SlowDecCoeffStandard (Fig. H1).

Figure 3 Posterior probability distribution (vertical bars) of 14 C submodel parameters calibrated over the range of the prior probability distribution (density curve). The posterior is depicted as the distribution of accepted values determined by random sampling from 10⁵ proposed values. The vertical red lines represent the optimal values as the median of the posterior and the blue lines the initial values of each chain.

3.3 Model performances

The statistical methods showed that the default MONICA model (MONICA_def) was inferior to the RothC, ICBM and C-TOOL models and produced results similar to CENTURY and CCB (Tables 3 and I1; Fig. 4). The mean absolute error was the highest of all the models but was within the margin of error compared with CENTURY and CCB for all sites and with CCB for the validation sites. Kendall's rank correlation coefficient and Nash–Sutcliffe model efficiency coefficient values were within the margin of error for all models and sites and therefore showed insignificant differences in the predictive skill of the models.

Table 3 Statistics describing the performance of the default (MONICA_def), the modified uncalibrated (MONICA_mod) and the modified calibrated (MONICA_cal) soil carbon submodel in MONICA in comparison with five C turnover models for all sites. The first 10 years were excluded from the analysis (n=223). The following statistics are depicted are Kendall's rank correlation coefficient (T), Nash–Sutcliffe model efficiency coefficient (NSE), mean absolute error (MAE) and the slope and intercept for characterising the relationship between measured and simulated C stocks.

${}^{***}$ p<0.001.

MONICA_def estimations deviated on average by more than ∼ 10 % from the measurements. Inclusion of a C content parameter for different crop residues and organic fertilisers in the MONICA model (MONICA_mod) improved the accuracy of predictions for C stock changes by ∼16 % for all sites and by ∼8 % for the validation sites, while the additional global parameter calibration (MONICA_cal) improved MONICA's average estimation error by ∼29 % for all sites and by ∼30 % for the validation sites. Compared to MONICA_cal, only RothC achieved a lower MAE for all sites. The differences between RothC and MONICA_cal were within the margin of error for the validation but not for all sites. MONICA_cal, MONICA_mod, RothC, ICBM and C-TOOL produced errors within the average standard deviation between the replicates of the C measurements for the validation (7.76 Mg C ha⁻¹) and for all sites (5.01 Mg C ha⁻¹). In contrast, MONICA_def, CENTURY and CCB had errors that were not within the average standard deviation for all sites. Each model was capable of producing results corresponding to the 95 % confidence interval for the validation sites (6.4 Mg C ha⁻¹). No model achieved the same for all measurements (3.3 Mg C ha⁻¹).

Based on the outcomes of the NSE method, all models performed better than the mean of the measurements, with MONICA_def achieving the lowest and RothC the highest scores. The differences between MONICA_cal and MONICA_def were within the margin of error. The Kendall rank coefficient displayed a significant statistical dependence between the simulated and measured data, with MONICA_cal, RothC and ICBM achieving the best results and CCB and CENTURY the weakest results. Overall the T coefficients were high, with MONICA_def achieving average results and MONICA_cal above-average results, which is in contrast with the other models. Of the six models, MONICA_cal generally produced better results than CENTURY, C-TOOL and CCB. The slope and intercept of the statistical analysis indicated that MONICA_def generally overestimated decomposition, while MONICA_cal produced results closer to the measured data. Overall, MONICA_def underestimated the SOC stocks at 24 sites and overestimated them at two sites, leaving 20 sites that were simulated within the standard deviation of the measurements. In comparison, MONICA_cal underestimated the SOC stocks at six sites and overestimated them at nine sites, leaving 31 sites that were simulated within the standard deviation of the measurements. The calibration resulted in a model behaviour that was comparable with RothC.

Figure 4 Linear regression between simulated and observed SOC stocks for 46 agricultural long-term monitoring sites, all years and 333 measurements for the MONICA (a), RothC (b), C-TOOL (c), ICBM (d), CENTURY (e) and CCB (f) models.

4 Discussion

4.1 Differences in models and performances

The MONICA model featured a plant growth, nitrogen and soil water submodel, while RothC, CCB, ICBM, CENTURY and C-TOOL did not. Smith et al. (1997) have argued that coupling C turnover with agroecosystem submodels will be necessary to simulate the impact of global change on agricultural production and SOC changes. However, they also assumed that simpler models have an advantage in simulating various land uses due to the lack of adequate data and the necessary parameterisation of complex models. Other literature shows that C turnover models coupled with sophisticated plant growth models do not result in a better performance, highlighting that complex models are more error prone (Mueller et al., 1996; Parton, 1996; Li et al., 1997). One explanation could be that some processes, such as the effect of phosphorus and other nutrients or the impact of plant diseases on crop growth, are not included in the models. In addition, information on crops and plant residues is often not available, so that measured values from other publications must be used to determine the required parameter values for crop development and ultimately the C inputs. This also applied to the data requirements of the MONICA model, which led to the necessary parameterisation of new crop and fertiliser types in order to simulate all sites. Since the data available for this publication did not contain additional information on the organic fertilisers and crops, allometric functions and values of differing genotypes and growing conditions had to be used to compare ingoing C inputs. Although we were able to adequately replicate yields as a proxy for modelling C inputs from crop residues, some discrepancies between model and observed yields remained. This could have led to an inaccurate reproduction of C inputs, since the specific properties of organic fertilisers and crop residues are characterised by high standard deviations (Möller and Schultheiß, 2015; Wiesler et al., 2016; Ma et al., 2018). In contrast, the C input estimation for the simpler models was based on observed crop yields for each site and allometric functions evaluated for Germany in previous studies (Dechow et al., 2019; Riggers et al., 2019). We assumed that a C input estimation based on measured yields should better describe the variability in the C inputs between sites. However, the results of the present study did not confirm the assumption that model complexity is an indicator of predictive ability in process-based models. Notable differences were driven by model characteristics, like SOM pool initialisation or model-specific response functions, some of which were controlled by parameterisation. Each model showed weaknesses depending on the specific conditions at each site (Figs. J1 to J7). For instance, MONICA_def overestimated decomposition in ∼52 % of all sites in the dataset but managed to predict the C losses at the Fuhrberg and Vinnhorst sites better than MONICA_cal. These sites had undergone recent land use change and had on average higher initial SOC stocks and losses over the measured period. In temperate agroecosystems, conversion from grassland or forest to arable land is considered to be one of the main reasons for C release (Vos et al., 2019; Poeplau et al., 2020). Correspondingly, sites with recent land use change showed the strongest decline in C stocks and were therefore a good indicator of the downward trend in SOC stocks. Removing sites with a cultivation period <50 years (n=10) resulted in a balanced C budget (0.02 t C ha⁻¹ a⁻¹), which also produced an even better performance by the calibrated MONICA model. The magnitude of changes in SOC stocks for all sites was comparable with other studies conducted in temperate climates, representing a general trend for German arable land (Steinmann et al., 2016; Sanderman et al., 2018; Vos et al., 2019; Seitz et al., 2022).

All models underestimated the magnitude of the variation in the natural SOC dynamics. The best simulation results were generally achieved at sites with slowly changing SOC stocks. Some sites distinguished by dramatic changes in SOC were especially problematic for the models to simulate, in particular the Dinklage site (BDF033_L), which was characterised by high C losses, and in contrast the Reinhausen site (BDF051_L), which was characterised by increasing C stocks. CENTURY and CCB, which under-performed in comparison with the other models, were the best at predicting these highly variable sites. The other models overestimated the SOC stock for BDF033_L and underestimated it for BDF051_L. In MONICA_cal, these two sites caused ∼10 % of the mean absolute error.

At site BDF033_L, the farmer was cultivating winter barley and corn as energy crops but experienced declining corn yields over the measured period. This probably led to a negative C balance, even though the field was extensively fertilised with pig and bull slurry. In MONICA, organic fertilisation was the best indicator of SOC changes, while there was no significant relationship between the amount of organic amendments and SOC stock changes in the measured data. Accounting for the quantitative variations in manure application supports the assumption that organic fertilisation is a weak predictor of SOC turnover rates and was correspondingly overestimated in the model. However, numerous studies clearly show a positive effect of organic amendments on SOC stocks (Koishi et al., 2020; Gross and Glaser, 2021; Roß et al., 2022). Since no control treatments were included in the experimental design of the long-term monitoring sites, it was not possible to carry out a comparative analysis of the individual sites with regard to the impact of organic inputs. Another consideration concerned site-specific uncertainties in the necessary model input. The properties of organic fertilisers and residues were unknown and thus had to be supplemented with standard values from the literature. Each fertiliser and crop type had fixed C ratios in the models, while in reality these inputs have a highly variable composition and can vary greatly depending on farm management, animal diet, climate and soil conditions (Schnug et al., 1996; McCartney et al., 2006; Cajamarca et al., 2019). There are three possible reasons for the disparity between the models and observations at site BDF033_L: (i) the fertiliser and crop-specific nutrient contents were inaccurate, leading to an overestimation of C inputs; (ii) the C losses from harvesting were underestimated; or (iii) the site and environmental characteristics and their effect on SOC decomposition were incorrectly simulated.

In comparison, at site BDF051_L, the farmer had transitioned to a management with reduced tillage intensities and no cultivation of root crops and managed to increase the C content in the upper soil layer by ∼30 Mg C ha⁻¹ during the measured time period. Although the effect of conservational tillage is still being debated, numerous studies support the effect measured at site BFD051_L (Freibauer et al., 2004; Luo et al., 2010). In the MONICA model, tillage uniformly distributed SOC in the affected soil layers. This was supported by the sensitivity analysis, which showed that tillage had an effect when there was heterogeneity in the soil profile. However, based on this study's data, it is likely that the model underestimated the effect of tillage on SOC decomposition. Possible reasons for the disparity were either an imbalance of C inputs into the topsoil or missing functions that describe the processes involved, such as the effect of tillage on soil aggregation, hydraulic properties and root growth (Mondal et al., 2020). The effect of tillage on SOC stocks is currently absent from other C turnover models, such as RothC. Jordon and Smith (2022) suggest modifying the decomposition rate to achieve a realistic representation of the effect of tillage on SOC, although they assume a small effect of tillage on SOC stocks.

Another probable reason for the differences between the model and measurements was the climate variables. In the data, only precipitation correlated significantly with SOC stock change rates, while in MONICA this was important but less so than temperature. Sierra et al. (2015) argue that functions in process-based models tend to overestimate the temperature effect on decomposition rates, and this imbalance cannot be offset by the moisture effect. This was also likely in the MONICA model, since the effect of temperature on decomposition was generally twice as strong as moisture. The importance of climate variables is still being debated, with Luo et al. (2017) and Carvalhais et al. (2014) finding major effects on SOC changes, while only minor effects have been established by Fujisaki et al. (2018) and Vos et al. (2019).

4.2 Implications for improving MONICA's soil carbon turnover model

Most biogeochemical models were developed many decades ago. This is also true of the C turnover subroutine in MONICA, which is based on the Daisy model (Hansen et al., 1991; Abrahamsen and Hansen, 2000). Therefore, it is based on the conceptual design of multiple SOC pools with pool-specific turnover rates that are not associated with clear, measurable soil properties (Bruun and Jensen, 2002). This circumstance makes it difficult to clearly define how the model could be optimised using response functions but also the initialisation of the model. Bruun and Jensen (2002) argue that it can be important to use plausible assumptions to initialise the SOM subroutine using radiocarbon measurements or the land use history. Nonetheless, functions in C turnover models need to be improved in order to estimate SOC dynamics across a wide range of spatiotemporal and agronomic conditions (Smith et al., 1997). Correspondingly, general insights into soil processes are needed to justify the introduction of additional or improved functions. There needs to be further examination of existing mechanisms and the interactive effects of site-specific characteristics but also missing factors that affect microbial abundance, diversity and activity (Louis et al., 2016). Microorganisms in particular play a significant role in the soil carbon cycle by controlling the rate of mineralisation and promoting plant growth by improving mineral nutrition or stimulating immune and stress reactions in crops (Tardy et al., 2015; Berg et al., 2017). This could be advantageous for representing the processes like priming or SOC responses on warming. However, microbial models, as reviewed in Chandel et al. (2023), vary widely in terms of model structure, processes considered and parameters required, which makes it clear that there is currently no consensus on which approach is best suited to represent the biosphere. Some of the models are evaluated and compared in Sulman et al. (2018), using experimental data from litter input studies and warming experiments. In evaluating these models, they find high variability in model results and conclude that first-order models already produce divergent projections due to parameter uncertainties and that structural diversity among models would exacerbate these uncertainties. We believe that increasing differences between models due to structural diversity does not represent a degradation in predictive ability but rather a more accurate estimate of predictive uncertainty (Bradford et al., 2016; Lovenduski and Bonan, 2017). It is unclear if the adoption of one of the various microbial approaches could increase the model accuracy of MONICA. Because of the increased model complexity and number of parameters, the evaluation and calibration of an incorporated microbial model for regional conditions would also require a more constraining dataset, preferably including those events where microbial models might outperform models based on first-order kinetics.

This study attempted to introduce two changes to the default MONICA model: (i) a revision, parameterisation and calibration of the functions that describe the influence of soil temperature, water and clay content on decomposition and (ii) an integration of a crop- and fertiliser-specific C content parameter. Alteration of the temperature function involved mathematically reformulating the original exponential function into a surge function. The reason for this change was that the default equations did not consider the decreasing activity of soil microorganisms at temperatures above 30 ^∘C (Liu et al., 2018; Fang and Moncrieff, 2001). This could lead to an overestimation of mineralisation, especially under globally increasing temperatures caused by climate change (DWD, 2021). Further improvements could be achieved by incorporating changes in the optimal temperature for decomposition based on site-specific varieties of microbial communities and C qualities (Fierer et al., 2006).

The clay function was changed to an inverse logistic function on the basis of the original equation. Many studies suggest that increasing soil clay contents have an adverse effect on SOC decomposition due to the reinforced inaccessibility of organic substances for decomposers (Six and Paustian, 2014; Churchman et al., 2020; Liddle et al., 2020). Adversely, several studies have found no effect or even a positive effect (Wei et al., 2014; Fissore et al., 2016; Zhang et al., 2021). There are multiple explanations for this disparity, e.g. different clay types, residue qualities and microbial diversities (Nguye and Marschner, 2014; Fissore et al., 2016). Without knowing the exact underlying processes to explain the variability and limitation of the available data, only minor changes were made when modifying the clay function. Removing the clay factor might reduce the complexity of the model but could also negatively impact the predictive capability of the SOC turnover routine. Further evaluation is required to assess the importance of the clay factor for decomposition and to credibly change the function.

Furthermore, the soil moisture factor, which was originally based on the pressure potential derived from the van Genuchten model (Van Genuchten, 1980; Abrahamsen and Hansen, 2000), was simplified. The default moisture function was directly controlled by the soil water content, soil texture, bulk density and SOC stocks, while the revised Gaussian function was only regulated by the soil water content. It is argued that the effect of texture on the decomposing community could be overestimated, since clay as a factor was already included in a function limiting mineralisation. In addition, the influence of soil texture on decomposition proved to be inconsistent, therefore requiring a better understanding of underlying processes (Doran et al., 1990; Scott et al., 1996), while the influence of bulk density was limited by model characteristics as it remained static after initialisation, hampering the reproducibility of management-related soil disturbances on the moisture function. Linn and Doran (1984) showed that the effect of ploughing on bulk density increased the microbial activity and decreased the optimum soil moisture content, while other studies have also determined an impact of soil structure changes on the mineralisation rate (Franzluebbers, 1999; Vilkiene et al., 2016). Another particularity of the default moisture function was that the optimal water-filled pore space for microbial activity was estimated to be between 10 %–25 % for sandy soils and 30 %–35 % for clay soils, while in the literature a general value of ∼60 % is considered the most probable (Linn and Doran, 1984; Doran et al., 1990; Liebig et al., 1995; Aon et al., 2001; Gabriel and Kellman, 2014; Zhang et al., 2015). Nevertheless, the simplification is probably less error-prone and easier to understand but at the same time not optimal for depicting the variance in microbial diversity and their environmental preferences. The inclusion of a more complex moisture function should be based on definite results.

As shown with the crop- and fertiliser-specific C content parameter, increasing the complexity of a model can be important for improving performance. According to the sensitivity analysis and the results, the addition of a C content parameter was necessary for the representation of realistic C inputs. However, the limited availability of data hinders the determination of precise C contents in plant residues and organic fertilisers. As mentioned in Sect. 4.1, values adopted from the literature could have a detrimental effect on the model's predictive capability. Therefore, measurements of C contents in crop residues and organic fertilisers should be part of the experimental design when investigating changes in SOC stocks. To further optimise the model, it is recommended that the C inputs from crop residues be divided into aboveground and belowground biomass. Currently there is no distinction between roots, shoot and leaves in the AOM pool, even though there is good evidence that roots in particular have a much longer mean residence time in the soil (Poeplau et al., 2021; Sokol and Bradford, 2019).

The above suggestions for improvement relate solely to the soil carbon model, but in complex models such as MONICA, the different subroutines can be interdependent, and changes in one could have a relevant impact on another function. Especially with respect to the feedbacks between the C and N cycle, plant growth and soil water contents should be considered in future studies.

5 Conclusions

The evaluation of the original MONICA model revealed weaknesses in estimating SOC trends. Using the default parameterisation, decomposition was overestimated, and the statistical performance was inferior compared with the C turnover models RothC, ICBM and C-TOOL. MONICA's capabilities of simulating SOC trends increased when more realistic functions were introduced that described the effects of soil temperature, moisture and clay on decomposition, coupled with a better specification of C inputs from organic fertilisation and crop residues. A sensitivity analysis displayed the importance of C inputs and the environmental factors of temperature and precipitation, as well as the corresponding parameters on SOC variation in the model. A successful calibration by means of Bayesian inference improved the parameterisation of the C turnover subroutine. The modified and calibrated MONICA model outperformed CENTURY, C-TOOL and CCB and produced similar results in comparison to RothC and ICBM. With the exception of certain sites, adequate reproduction of SOC stock change rates was achieved. However, none of the investigated models was capable of simulating each site satisfactorily. All the models underestimated the variation in measured SOC dynamics that occurred at sites with particular properties and management practices. This paper reveals that even complex, biogeochemical models are capable of performing similarly or even better than the more simplistic C turnover models. It is likely that with better data availability and continuing optimisation of functions representing natural processes, MONICA's performance can be further enhanced.

Appendix A: General properties of each monitoring and experimental site

Table A1 Description of the long-term monitoring sites (n=46) used for calibration and validation of the carbon turnover routine in the MONICA model. SOC is soil organic carbon. MAT is mean annual temperature. MAP is mean annual precipitation.

^a Sites used for model validation (n=15). ^b Predominantly applied fertilisers; N values for organic amendments were estimated. UAN is urea ammonium nitrate. U is the urea. CAN is calcium ammonium nitrate. AS is ammonium sulfate. AN is ammonium nitrate. POM is poultry manure. PIS is pig slurry. PIM is pig manure. LFR is liquid fermentation residues. SS is sewage sludge. CAS is cattle slurry. CAM is cattle manure. SHM is sheep manure. CP is compost. HOM is horse manure.

Table A2 Description of the long-term experimental sites (n=11) used for the validation of simulated soil temperatures and water contents. SOC is soil organic carbon. MAT is mean annual temperature. MAP is mean annual precipitation.

Appendix B: Revised environmental decomposition factors

Figure B1 The effect of soil temperature on the soil organic carbon decomposition rate. The solid red line describes the standard function that is integrated in the MONICA model. The solid blue line describes the revised function as the best fit curve to the observed data, while the dashed blue lines describe the possible properties of the equation according to the prior probability distribution of the parameters QTenFactor and TempDecOptimal.

Figure B2 The effect of soil moisture on the soil organic carbon decomposition rate. The solid red line describes the standard function for a clay soil (clay = 60 %, sand = 10 %, bulk density = 1.4, and soil organic carbon = 2 %), while the solid yellow line describes the standard function for a sandy soil (clay = 4 %, sand = 80 %, bulk density = 1.4, and soil organic carbon = 1 %). The solid blue line describes the revised function as the best fit curve to the observed data, while the dashed blue lines describe the possible properties of the equation according to the prior probability distribution of the parameter MoistureDecOptimal.

Figure B3 The effect of clay on the soil organic carbon decomposition rate. The solid red line describes the standard function that is integrated in the MONICA model. The solid blue line describes the revised function as the best fit curve to the observed data, while the dashed blue lines describe the possible properties of the equation according to the prior probability distribution of the parameter LimitClayEffect.

Appendix C: Morris elementary effects screening method

Table C1 Mathematical formulation of the sensitivity measures of the Morris elementary effects screening method, where r is the number of input factors, and EE is the elementary effect of each input.

Appendix D: Metropolis–Hastings algorithm

Table D1 Mathematical formulation of the Metropolis–Hastings algorithm and the log–Laplace likelihood function, where θ is the parameter vector, Y the dataset used for calibration, t the number of iterations, σ the standard deviation inferred jointly from the error between observation and simulation, ω the model input, f(ω_i;θ_i) the model run with the parameter vector θ, α the acceptance ratio and δ a random vector from a multivariate distribution.

Appendix E: Model evaluation methods

Table E1 Quantitative methods used for data and model evaluation, where S is the predicted value, O is the observed value, n is the number of cases, and k is the number of divided subgroups from n.

Appendix F: Results of the Morris elementary effects screening method

Table F1 Results of the Morris method depicting the effects of 51 variables on soil carbon changes. The variables are subdivided into site, management, crop- and fertiliser-specific and carbon turnover factors. The parameters CorgContent (C content of crop residues and organic fertilisers), MoistureDecOptimal (optimal water-filled pore space for decomposition), TempDecOptimal (optimal temperature for decomposition) and QTenFactor(Q₁₀ temperature coefficient) were included as modifications to the model. The overall importance of an input factor is specified by two sensitivity measures, namely the mean (μ^*) and the standard deviation of the elementary effect (σ) of each parameter. Higher values in the two sensitivity measures imply a greater importance of the variable for changes in soil carbon stocks. The uncertainty is presented as the standard deviation between the results of three separate Morris elementary effects screening runs.

Appendix G: Gelman and Rubin shrink factor

Figure G1 Evolution of the Gelman and Rubin shrink factor. The 50 % (solid line) and 97.5 % (dashed red line) quantiles show the convergence of the calibration.

Appendix H: Correlations between C turnover parameters

Figure H1 Correlation matrix quantifying the statistical relationship between the prior probability distributions of the SOM submodel parameters.

Appendix I: Model performance

Table I1 Statistics describing the performance of the default (MONICA_def), the modified uncalibrated (MONICA_mod) and the modified calibrated (MONICA_cal) soil carbon submodel in MONICA in comparison with five C turnover models for 15 validation sites. The first 10 years were excluded from the analysis (n=70). The statistics depicted are Kendall's rank correlation coefficient (T), Nash–Sutcliffe model efficiency coefficient (NSE), mean absolute error (MAE) and slope and intercept for characterising the relationship between measured and simulated C stocks.

${}^{***}$ p<0.001.

Appendix J: Model performances

Figure J1 Performance of the default MONICA model for each individual site. The first 10 measured years were excluded.

Figure J2 Performance of the modified and calibrated MONICA model for each individual site. The first 10 measured years were excluded.

Figure J3 Performance of the CCB model for each individual site. The first 10 measured years were excluded.

Figure J4 Performance of the CENTURY model for each individual site. The first 10 measured years were excluded.

Figure J5 Performance of the C-TOOL model for each individual site. The first 10 measured years were excluded.

Figure J6 Performance of the ICBM model for each individual site. The first 10 measured years were excluded.

Figure J7 Performance of the RothC model for each individual site. The first 10 measured years were excluded.