Bayesian calibration of the ICBM/3 soil organic carbon model constrained by data from long-term experiments and uncertainties of C inputs

Abstract

Models with various complexity can asses soil C sequestration in agriculture. In this study, we updated the Introductory Carbon Balance Model (ICBM) with 28 years of additional data and included multiple long-term bare fallow experiments for old pool kinetics. We validated the model with data from a new long-term sister experiment. The new calibration included uncertainty in the estimation of below-ground C inputs to soil. The model now considers above- and belowground and external C inputs separately (ICBM/3). The underlying mathematical approach is the same, with two state variables and a climatic parameter, and such simple structure remained robust enough to describe soil organic C dynamics over six decades. We tested including an inert soil C pool in the model structure, but it did not decrease the observed variance. Similarly, an intercept in the functions for estimating belowground C input from crop yield was not useful. Results suggest that root C contributes more to the old organic C pool than aboveground C inputs. We also evaluated parameters interactions, in particular between C inputs and their transformation into more stable soil C interacted and decomposition kinetics. We also describe new functions for estimating the ICBM climatic parameter in a more user-friendly way.

Keywords:

• SOC decomposition
• SOC models
• Bayesian uncertainty estimation
• ICBM model

Introduction

Soil organic carbon (SOC) is the primary terrestrial C pool, and its management is crucial to the global C budget in the biosphere [1]. To manage it, several countries use models with different degrees of complexity to simulate changes in soil C stocks to determine their national contributions to global greenhouse gas balances [2]. The robustness of such predictions depends on the accuracy of models and is unavoidably associated with uncertainties. Uncertainty in SOC models is related to imperfect knowledge of processes, model structure, errors in the data from long-term field experiments (LTEs) used for their calibration, and uncertainty in input data. All these components reduce the robustness of predictions. While structural uncertainty is inherent to specific model assumptions, other uncertainty sources can be estimated with stochastic calibrations [3].

The annual C input to soil is the most important driver for changes in SOC stocks, and consequently, it is, in all models, one of the largest uncertainty sources [4,5]. This C input comes from exogenous organic amendments and aboveground (AG) and belowground (BG) crop residues; the latter comprises the root biomass left in the field at harvest plus root-derived C delivered to the soil during the growing season, usually referred to as rhizodeposition (RD) [6]. Many SOC models estimate such inputs with empirical approaches (i.e. allometric functions) based on shoot-to-root (SR) ratios to estimate the C input from root biomass left in the field at harvest. These ratios consider shoots as the total AG biomass produced, estimated from grain yields using data on harvest index (i.e. the ratio of grain to total AG biomass; HI) for annual crops [4]. Most models assume that a fraction of such C inputs is converted into more stabilized soil organic material (i.e. a humification coefficient; h) entering one or several SOC pools depending on model structure [7].

Of the large uncertainties associated with C inputs, the most important are relative to estimating BG C inputs. For annual crops, pioneering work by Barber et al. [8] and by Sauerbeck et al. [9] indicated an RD coefficient of 1.0 (i.e. RD contributes as much as the root biomass to BG C inputs), an assumption commonly used thereafter [10,11]. Because these estimates are affected by traits of plants, pedo-climatic conditions and fluctuate over the growing season, RD coefficients reported in the literature vary from much lower than 0.5 to as high as 2.0 [12,13], evidentiating the high uncertainty around root C inputs estimation. We also know that field-based estimates of SR ratios for crops are highly variable within different crops, such as small-grain cereals and maize, with coefficients of variation typically of at least 50% [6,14].

Furthermore, SR ratios and HI data vary with genotypes and can change over time due to plant breeding; this can cause problems when modeling LTEs or assessing historical changes in C inputs [15]. Finally, the most common allometric functions assume that SR ratios and HI data are linear with respect to production levels [4]. However, some recent studies suggest that estimating BG C inputs using fixed values for root biomass by crop type may be better, independent from yields [16,17] or using HI values that change with crop yields [18].

All these uncertainties affects the uncertainty of the humification (h) coefficients used in most SOC models, which are not measured but estimated with model calibration. The amount of C inputs to the soil and their decomposition (controlled in large part in compartmental SOC models by the h coefficients) interact directly in the models to determine the simulated SOC stocks over time. The change in SOC stocks over time in LTEs, the primary source of information to calibrate models, can result from many combinations of C inputs and humification. These direct interactions between the kinetics of decomposition and the C inputs, neither of which is ever precisely known, imply that many parameter combinations can explain the observed results (which is defined as equifinality) in the models [19] and make it very difficult to determine both SR and h values.

Historically, commonly used values of h for AG crop residues (h_s) vary between 0.05 and 0.15 [20–22]. However, the h-values for roots (h_r) were early recognized as higher than those for h_s. For example, Hénin and Dupuis [20] assumed the h_r coefficient for maize and small-grain cereals were 1.50 and 1.80 times higher than the corresponding values for h_s. Since then, several studies on LTEs and reviews have been quantifying this further, showing that total root-derived C (i.e. roots and RD) contributes from 1.3 up to almost four times more to SOC than AG [11,14,23]. Following these studies, some models have recently assumed a higher root contribution to SOC than contributions from AG vegetation. For example, Dechow et al. [24] modified the default h-value used for both shoots and roots (0.47) in the RothC model and found that using a lower range of values for h_s (0.27–0.35) compared to h_r (0.43–0.47) improved the model accuracy. On the other hand, in the AMG model, Clivot et al. [25] considered h_r values to be 26% to 77% higher than that for h_s.

Sweden uses a relatively simple first-order two-component model, Introductory Carbon Balance Model (ICBM) [26], to estimate SOC stock changes in the national reporting system [27]. This approach relies on linear allometric functions for estimating plant C inputs from yields. It includes a parameter modifying SOC decomposition rates as a function of climate, edaphic properties, and crop conditions. The ICBM model is also used within farm-scale advisory tools such as Holos [28] and other applications (Supplementary Material, Table 1). Recent studies of the Ultuna LTE, used to calibrate the original version of ICBM, and its sister LTE at Lanna [29,30] showed that roots contribute more to SOC than AG crop residues. Based on these findings, some attempts have already been made to modify the ICBM model [31,32] by dividing the crop residue C input pool into AG and BG components, assuming h_r 2.3 times the value of h (i.e. 0.125) in the original calibration. Kröbel et al. [28] suggested the introduction of another third component accounting for organic amendments (i.e. ICBM/3 version), recognizing the need to calibrate and validate this new model structure.

Table 1. A summary of the treatments (•) and topsoil texture for the long-term field experiments used in the calibration and validation of ICBM/3B.

Treatments	Ultuna^a	Lanna^b	Versailles^c	Grignon^c	Askov B3^c	Kursk^c
	Calibration	Validation
Bare fallow	• (A)	• (A)	•	•	•	•
Unfertilized control	• (B)	• (I)	–	–	–	–
Calcium nitrate	• (C)	• (B)	–	–	–	–
Ammonium sulfate	• (D)	• (C)	–	–	–	–
Calcium cyanide	• (E)	–	–	–	–	–
Straw	• (F)	–	–	–	–	–
Straw + N	• (G)	–	–	–	–	–
Green manure (grass)	• (H)	• (D)	–	–	–	–
Peat	• (I)	–	–	–	–	–
Farmyard manure	• (J)	• (E)	–	–	–	–
Farmyard manure + P	• (K)	–	–	–	–	–
Sawdust	• (L)	–	–	–	–	–
Peat + N	• (M)	–	–	–	–	–
Sawdust + N	• (N)	–	–	–	–	–
Sewage sludge	• (O)	• (F)	–	–	–	–
Sewage sludge + metals	–	• (G)	–	–	–	–
Compost	–	• (H)	–	–	–	–
Topoil texture (%)
Clay	36.5	45.5	17	30	7	30
Silt	41	46.5	57	54	11	65
Sand	22.5	8	26	16	82	5

^a For details, see Kätterer and Andrén [34], in parenthesis, the treatment ID used in this publication.

^b For details, see Kätterer et al. [30], in parenthesis, the treatment ID used in this publication.

^c Sites from the long-term bare-fallow network are described in detail in Barré et al. [33].

Furthermore, since its development, the original ICBM calibration still needs to be updated. More than two decades of new data are now available from the Ultuna LTE, and the sister experiment at Lanna, starting just after the original calibration, now has more than two decades of data. Moreover, throughout the years, we have further developed the functions for estimating the climatic parameter of ICBM, and such incremental improvements also need to be assessed and validated. The present study provides an update of the model, together with additional tools, such as the climatic module, that allow the model users to adapt this updated ICBM version to their workflow quickly and maintain the model being easily accessible. Given the expected equifinality of the model, we relied on a stochastic calibration approach, which should be more robust. Our objectives in the present study were as follows:

• Calibrate the ICBM/3 version on the Ultuna LTE using more than two decades of new data within a stochastic Bayesian framework, including the new functions for estimating the climatic parameter. This version is referred to as ICBM/3B (B for Bayesian).

• Test whether the longer time scale requires modifying ICBM/3B by including an additional inert SOC pool.

• Test the robustness of a linear allometric approach to BG input estimations by adding an intercept in the allometric functions representing a minimal total BG C input independent of crop yield.

• Analyze model uncertainties and equifinality, focusing on BG C inputs, including the available information for constraining SR ratios, hs, and hr from the literature.

• Validate ICBM/3B on the Lanna sister LTE, the decomposition of the old SOC pool in ICBM/3B, and the new functions for estimating the climatic parameter on European bare-fallow LTEs.

Materials and methods

LTEs used for calibration and validation

The two Swedish LTEs are located along an east-southwest gradient in central Sweden. Experiment treatments involve the effect of three different types of mineral N fertilizer and organic amendments (farmyard manure FYM, composts, sewage sludge SLU, sawdust SAW, peat, small-grain cereal straw, and green manure GM) or a combination of both. The Swedish classical frame trial (FRAME-56) was initiated in 1956 in small plots (2 × 2 m, separated by metal frames) at Ultuna [29], while its sister experiment (R3-0130) at Lanna started in 1996 in large (92 to 112 m²) plots [30]. The latter includes several treatments in the former site; both LTEs also include a bare fallow treatment (Table 1). Both sites replicate all treatments four times in randomized block designs, and the organic amendments are applied every second year, corresponding to approximately 2 Mg C ha⁻¹yr⁻¹ in the frame trial and 4 Mg ash-free dry mass ha⁻¹yr⁻¹ at the Lanna site.

We used the Ultuna LTE (59.82^oN, 17.65^oE) time series for all model calibrations (Figure 1). The field had been cultivated for centuries before the start of the experiment and is located on a clay loam, classified as a Typic Eutrochrept (USDA soil taxonomy) or Eutric Cambisol (FAO). The annual mean total precipitation is 542 mm, and the mean annual temperature is 5.8 °C. The organic amendments are added in the fall by carefully incorporating them into the topsoil down to a depth of 20 cm using a spade. The frequency of crops (number of years) between 1956 and 1999 was as follows: small-grain cereals (31), fodder rape (5), oilseed (3), fodder beet (2), and crop failure (3). Small-grain cereals were mostly spring barley (17) and oats (11), spring wheat (2), and winter wheat (1) were occasionally grown. Silage maize has been growing continuously since the year 2000. Every year at harvest, all AG plant biomass is cut at the soil surface, removed from the experimental plots, weighed, and sorted into crop products and AG crop residues; see Kätterer et al. [29] for more details. The climatic data come from the Ultuna meteorological station, located at about 1 km distance from the experimental field.

Figure 1. Soil organic carbon (SOC) stock time series for the treatments of the Ultuna LTE, on which we based all model calibrations.

The Lanna LTE (58.34^oN, 13.10^oE) was used for validating the model and is located on an Aquic Haplocryept developed on a Quaternary silty clay deposit [30]. The annual mean total precipitation is 636 mm, and the mean annual temperature 7.3 °C. The organic amendments are added in the fall by incorporating them through moldboard plowing into the topsoil to a maximum depth of 25 cm. The frequency of small-grain cereal crops (number of years) between 1996 and 2018 was as follows: oats (12), spring barley (8), winter wheat (2), and spring wheat (1). At this site, the stubble was left in the experimental plots after harvest, corresponding to 33% of the total amount of straw produced, measured for most of the years; see Kätterer et al. [30] for more details. The climatic data for Lanna come from the Swedish Meteorological Institute (https://www.smhi.se/weather/Sweden-weather/obsevartions).

We also used data from an LTE bare-fallow network in Europe to specifically validate the decomposition rates of the old SOC pool and the climatic module of the ICBM model. From this network, consisting of six LTEs, we selected five, including Ultuna (year of initiation): Askov in Denmark (1956), Grignon and Versaille in France (1959 and 1928), and Kursk in Russia (1965). All sites have been kept free from vegetation, and we retrieved weather data records from https://www.ncdc.noaa.gov/cdo-web/; for details on these LTEs, see Barré et al. [33].

SOC measurements and equivalent soil mass calculation

Measurements of SOC in the Ultuna and Lanna LTEs were made intermittently, more or less every second year, in the topsoil down to a fixed depth of 20 cm. For the time series of SOC data from the bare-fallow LTEs, measurements were taken ten times at Versaille and Grignon and six times at Kursk to a fixed depth of 25 cm, while measurements were available for almost every year at Askov to a fixed depth of 20 cm [33]. The mass of soil sampled to a fixed depth changes over time proportional to dry soil bulk density changes. Consequently, we corrected the reference depth of 20 cm used in this study using an equivalent soil mass approach described in detail by Kätterer et al. [29]. In short, treatment-specific bulk density values over time were estimated using linear functions based on the initial value measured at the start of the trial (1.44 Mg m⁻³) and those measured in the different treatments in 2009. The depth to which an equivalent amount of soil mineral mass is distributed was normalized for 20 cm depth in the control treatment. We then calculated soil carbon stocks by multiplying the bulk density and SOC concentration product with equivalent soil depth per treatment and year. The equivalent soil depth in the treatments at the Ultuna LTE varied between 19.5 and 27 cm [29].

The ICBM/3B model

The model we calibrated is derived from the original ICBM presented by Andrén and Kätterer [26]. The original ICBM model is a two-component first-order compartmental model with two state variables referred to as Young (Y) and Old (O) SOC pools, having specific decomposition rates (k_Y and k_O, respectively), connected by a humification coefficient (h). The annual C input to soil enters through the Y pool, representing organic matter not yet humified with a fast decay rate, while the O pool represents the organic matter stabilized in the soil that decays much slower. In the original calibration, k_Y (here referred to as k₁) was 0.8 year⁻¹, and k_O (here referred to as k₂) was 0.0061 year⁻¹ under standard climatic conditions at the calibration site in central Sweden [26].

In order to consider C inputs from different organic materials having specific h-values, which determine the fraction that enters the O pool, we introduced three separate Y pools, all decaying in parallel (Figure 2). These pools were AG crop residue inputs (Y_s), total BG inputs including root tissues and rhizodeposits (Y_r), and inputs from amendments (Y_a), which are associated each with different humification coefficients (h_s, h_r, and h_{a, i}, respectively, where i refers to each specific amendment in the experiment). In this application, Y_a represented farmyard manure, sewage sludge, peat, and sawdust (Table 1). The basic mathematics of the ICBM model has been extensively introduced in the literature (Appendix A). In particular, we refer the reader to Appendix 1 in Kätterer and Andrén [34], which thoroughly describes the equations. This model, from now on called the standard ICBM/3B model to distinguish it from the original version, adds two additional Y pools as described above, which are associated with different humification factors to account for different qualities of C inputs. This modification does not alter the mathematics foundations of the model, and, under constant input conditions, it is equivalent to considering an h value that is the weighted average of the individual h values of the various material added to the soil. However, it allows to decompose the different h values in calibration under non-steady input conditions, resulting in a more useful model when applied.

Figure 2. A schematic representation of the ICBM/3 model, adapted from Menichetti et al. [ 3].

The climatic functions

The climatic interactions with decomposition are implemented through a climatic parameter multiplying the two kinetic coefficients, k_1, and k₂. This parameter derives from the product of relative water content (r_wat) and soil temperature (r_temp) calculated daily and averaged to an annual value since the model works in annual time steps. Its calculation involves daily climatic data and soil and crop properties. In this application, the climatic parameter was for convenience called r_clim (please note that this corresponds to r_{e_crop} as defined in Bolinder et al. [35]).

In this study, we updated the calculations of r_wat and r_temp with the most recent functions we are using, as described in detail below. Specifically, we introduced a function with a peak in the middle at optimal moisture and with activity decreasing on both sides to represent the effect of soil moisture on biological decomposition activity at soil saturation. We also present in more detail the set of functions and information flow involved in the calculations of r_clim (Appendix B), and we compiled the relative functions into an R package available on the GitHub repository of the corresponding author (address at the end of the Discussion section).

The effect of soil temperature is based on a quadratic response function [27] with maximum ($T_{\mathrm{\text{max}}}$) set to 30 °C and minimum ($T_{\mathrm{\text{min}}}$) set to −3.8 °C. $$r_{\mathit{\text{temp}},t}=\frac{{\left(T_{\mathit{\text{soil}},t}-T_{\mathrm{\text{min}}}\right)}^2}{{\left(T_{\mathrm{\text{max}}}-T_{\mathrm{\text{min}}}\right)}^2}$$

This function presents a curve quite similar to the more conventional Arrhenius function [36]. Its main limit of applicability is the lack of a protein denaturation mechanism, which varies depending on different soils but happens at relatively high temperatures (>35 °C [37]), which are in general not reached in the soil in temperate climates. The daily mean soil temperature is calculated from the daily mean temperature based on an empirical relationship: $$T_{\mathit{\text{soil}},t}=\mathrm{\text{max}}\left(-2,0.92\bullet T_{\mathit{\text{air}},t}\right)$$

The effect of soil water content is more complicated to simulate since it is seldom available, and in general, it must be calculated from meteorological data. The response function for soil moisture is in this present model version adapted from Moyano et al. [38]: $$r_{\mathit{\text{wat}},t}=f\left(x\right)=\{\begin{array}{c}0,\wedge {\theta }_{\mathit{\text{th}}}<{\theta }_m\\ \frac{{\theta }_t-{\theta }_{\mathit{\text{th}}}}{\mathrm{\text{min}}\left({\theta }_{\mathit{\text{opt}}},{\theta }_{\mathit{\text{fc}}}\right)-{\theta }_{\mathit{\text{th}}}},\wedge {\theta }_{\mathit{\text{th}}}\le {\theta }_t<\mathrm{\text{min}}\left({\theta }_{\mathit{\text{opt}}},{\theta }_{\mathit{\text{fc}}}\right)\\ 1,\wedge \mathrm{\text{min}}\left({\theta }_{\mathit{\text{opt}}},{\theta }_{\mathit{\text{fc}}}\right)\le {\theta }_t\le \mathrm{\text{max}}\left({\theta }_{\mathit{\text{opt}}},{\theta }_{\mathit{\text{fc}}}\right)\\ 1-\left(1-\rho \right)\frac{{\theta }_t‐\mathrm{\text{max}}\left({\theta }_{\mathit{\text{opt}}},{\theta }_{\mathit{\text{fc}}}\right)}{{\theta }_s‐\mathrm{\text{max}}\left({\theta }_{\mathit{\text{opt}}},{\theta }_{\mathit{\text{fc}}}\right)},\wedge ‐\mathrm{\text{max}}\left({\theta }_{\mathit{\text{opt}}},{\theta }_{\mathit{\text{fc}}}\right)<{\theta }_t\le {\theta }_s\end{array}$$ where ${\theta }_{\mathit{\text{opt}}},$ the optimal water content for decomposition, is calculated according to ${\theta }_{\mathit{\text{opt}}}=0.2+1-26\bullet {\left({\theta }_{\mathit{\text{fc}}}\right)}^{2.03}$ and ${\theta }_{\mathit{\text{tg}}},$ the lower boundary water content for decomposition, is calculated according to ${\theta }_{\mathit{\text{th}}}=0.0965\bullet \mathit{\text{ln}}\left({\theta }_{\mathrm{\wp }}\right)+0.3,$ ${\theta }_{\mathit{\text{fc}}}$ is the water content at field capacity, ${\theta }_{\mathrm{\wp }}$ is the water content at wilting point, and $\rho$= 0.5. However, this function needs, as input, besides soil hydrological constants, the soil water content. Soil water content is seldom available but can be calculated based on a simulated daily soil water balance once precipitation and evapotranspiration are known. For a specific soil layer of a given thickness, the water balance at a given time (${\Omega }_t$) can be calculated according to: $${\Omega }_t={\Omega }_t-1+{\mathit{\text{PPT}}}_t-{\mathit{\text{ET}}}_t-P_{\mathit{\text{intercepted}},t}-{\Omega }_{\mathit{\text{bypass}},t}$$ where ${\mathit{\text{PPT}}}_t$ is the total precipitation in a given day $t,$ ${\mathit{\text{ET}}}_t$ is the actual evapotranspiration on the same day, P_{intercepted,t} is the water intercepted by the crop and ${\Omega }_{\mathit{\text{bypass}},t}$ is the excess water (lost by percolation). The ${\mathit{\text{ET}}}_t$ is calculated from reference evapotranspiration (${\mathit{\text{ET}}}_0$) through two coefficients according to ${\mathit{\text{ET}}}_t={\mathit{\text{ET}}}_o{K}_c{K}_r,$ where $K_c$ is the coefficient to calculate the potential evapotranspiration on the specific plot from the reference (${\mathit{\text{ET}}}_c={\mathit{\text{ET}}}_o{K}_c$). The factor $K_c$ is calculated from the green area index (GAI) according to $K_c=1.5-0.5\bullet e^{-0.17\bullet \mathit{\text{GAI}}}.$ The intercepted water $P_{\mathit{\text{intercepted}},t}$ is calculated using the GAI according to $P_{\mathit{\text{intercepted}},t}=\mathrm{\text{min}}\left({\mathit{\text{PPT}}}_t,{\mathit{\text{ET}}}_{c,},\mathit{\text{GAI}}\bullet 0.2\right).$ The factor $K_r$ is a coefficient to consider the effect of water stress on the crop, and it is calculated according to $K_r=\mathrm{\text{max}}\left(0,{\left(1-\frac{0.9\bullet {\theta }_{\mathit{\text{fc}}}-{\theta }_t}{0.9\bullet T_{\mathit{\text{field}}}-0.7\bullet {\theta }_{\mathrm{\wp }}}\right)}^2\right).$ Wilting point (${\theta }_{\mathrm{\wp }}$) and field capacity (${\theta }_{\mathit{\text{fc}}}$) are calculated according to pedotransfer functions developed for Sweden [39].

The climatic parameter is calculated for each site and treatment separately and is normalized to the unit for the straw + N treatment (G) at the Ultuna site since this represents the most normal agricultural practice [26,40]. For this treatment, the updated calculations described above yielded a r_wat × r_temp factor of 0.1057 for the time series 1956 to 2019. Thus, we were scaling all r_clim values by dividing them by this factor, which implies that r_clim became 1 for treatment G at Ultuna, and r_clim for all other treatments at Ultuna, Lanna, and the long-term bare fallow sites have values relative to this normalized calibration site and treatment. Since r_clim is multiplying the kinetic coefficients of the ICBM SOC pools, this implies that a higher value increases decomposition and that for a given amount of annual C inputs to soil, the steady-state SOC mass will be lower.

Testing possible modifications of the ICBM/3B model

Since the original ICBM calibration [26], the trends in SOC stocks have been changing for most treatments. Consequently, it is expected that the decomposition rates have changed by adding two more decades of data, particularly in the bare-fallow treatment used initially for estimating the decomposition rate of the old SOC pool (k_O), assuming a single exponential decaying pool may be less appropriate. Indeed, on very long time scales, the quality changes of SOC may not be well described by a single exponential law because the assumed exponential decay rate is slowing down over time [33,41]. Increasing the number of pools, adding, for example, a very slow or inert pool, can mitigate this problem. We also calibrated the model by introducing an inert pool to test if we needed to update the original model structure with only two SOC pools (Y and O) for time scales longer than four decades. This is a simple approach for increasing model complexity to represent a broader range of SOC quality. The size of the inert pool was treated as a parameter in the Bayesian framework.

The long-term experiment at Ultuna is located more or less in an urban area, and in order to avoid attracting people visiting the site to poach corn on the cobs, the maize was seeded late (between May 30 and June 23) and harvested relatively early (between September 2 and 29). Therefore, the maize at this site did not fully exploit the potential growing season. Some of the treatment plots are also extremely acidic. Consequently, maize yields are relatively low, varying between 1.9 and 8.4 and averaging 5.5 Mg dry matter ha⁻¹ across all treatments from 2000 to 2019. Although maize primarily grows in southern Swedish regions, national average maize yields (2014–2018) usually are slightly higher than 10 Mg dry matter ha⁻¹ (www.scb.se/JO0601). Therefore, we tested if an intercept in the allometric functions representing a minimal total BG C input could be more appropriate for these conditions, particularly for maize. $$I=\left(1+\varrho \right)\bullet \gamma \bullet \left(\mathrm{\alpha }+y_w\right)\frac{1}{{S:R}_w}$$

Where I is the mass of C inputs from roots each year, $\varrho$ is the amount of roots exudates as a fraction of root mass, and we set it to 0.65 [6], $\gamma$ is the average ratio of C in the total organic mass and was assumed between 0.40 and 0.51 with uniform prior, $y_w$ and ${S:R}_w$ represents each crop’s yields and shoot-to-root ratio, respectively, and $\mathrm{\alpha }$ is the intercept. Units are in Mg C ha⁻¹. The intercept was treated as a parameter in the Bayesian framework (i.e. it is zero in the standard ICBM/3B version).

Calibration and initialization of the ICBM/3B model

The model was calibrated within a Bayesian framework in JAGS language [42] and ran in R [43]. The technique is nowadays broadly adopted, and we refer to specific publications for details, such as Gelman et al. [44] and Kruschke et al. [45]. Shortly, the approach is a stochastic implementation of Bayesian statistics, achieved by running a large number of simulations with parameter sets sampled from their prior probability distributions and then updating such distribution based on the information from the data. The properties of the sampler used (Metropolis-Hastings) are Markovian (meaning that each new state is correlated only with the former), and this makes it so that the mean of each resulting chain (Markov chain Monte Carlo, shortly MCMC) converges to the optimal calibration. At the same time, the distribution of the MCMC represents the updated posterior probability distribution. Assuming that some distributions might be skewed, we generally relied on the mode rather than the mean (still reporting both) for estimating the best prediction. The climatic parameter r_clim was left outside the Bayesian framework and utilized as an external driving variable. We also considered the RD coefficient (0.65) outside the Bayesian framework. We chose this because any uncertainty related to this parameter, being a multiplier of root biomass, will be included in the root biomass uncertainty since the two terms interact directly.

The initial ICBM calibration [26] was based on a sequential approach, where different treatments at the Ultuna LTEs were used to infer specific model parameters. More specifically, the bare-fallow treatment, not receiving C inputs since the beginning of the experiment, was used to determine the decay rate of the old SOC pool (k₂). This assumption implies that the decomposition of this pool was considered independent of C inputs, which may be questionable, particularly for some organic amendment treatments. In the present study, we relaxed this assumption, and instead of a sequential calibration, we calibrated all parameters simultaneously using all the treatments. We recognize that eventual modification of the kinetic coefficients can occur through SOC interactions with C inputs (e.g. priming). However, we assume they are not important enough to make the process of SOC decomposition radically different between treatments and that the Bayesian calibration can take care of minor modifications of the decay rates with its error model. This versatility is one of the advantages of the calibration method utilized in this study, together with estimating the uncertainty.

Model predictions were based on a combination of the whole MCMC parameter sets, uniformly resampled 5000 times over 25 independent chains. The model was running for all these 5000 × 25 sets, and the uncertainty intervals were drawn where the maximum and minimum of all these runs are present.

We initialized the model by including the proportion between the young and old (Y, O) SOC pools, which prior was assumed to be between 80% and 100% O (and the remaining Y), in the Bayesian framework. We assumed a fixed proportion for redistributing the C mass assigned to each Y pool (50% to the shoots, 50% to the roots, and none to the amendments). This way, the initialization of the SOC pools could be treated as any other parameters and calibrated within the Bayesian framework.

Since one of the study’s aims was to calibrate the model with more than two decades of new data, we conducted the calibration on two different time series. The first ended in 1999 (close to the original calibration, which ended in 1991), and the second ended in 2019. The latter allows us to assess the influence of silage maize since the shift from C3 to C4 crops that occurred in the year 2000.

Estimating model uncertainties: model priors

Initially, ICBM was designed in the beginning of the 1990s and aimed to produce a simple yet versatile first-order kinetic model for describing and predicting SOC stock changes [26]. While recognizing that these simple approaches are not explicitly accounting for all the mechanisms of SOC formation and turnover [46,47], they have proven efficient in capturing a large amount of the variance in SOC decay [41]. A large part of the model uncertainty is inherent to the structure common to all models in the same class (first-order, compartmental). On top of that, the assumption of a particular discretization carries some uncertainty, which we could include in the Bayesian framework as a parameter with a more abstract generalization of first-order models such as Q [41]. However, model structural uncertainty was outside the scope of this study, and we focus here only on parameter uncertainty given the assumed model structure.

For the model version with a constant inert pool, we assumed its priors conservatively as a uniform distribution between 0 and 15 Mg SOC ha⁻¹, considering that former estimates were between 2 and 7% of total SOC stocks [33].

Priors for the decomposition rates and humification coefficients

In the present study, prior probability distributions for ICBM parameters were mainly taken from the literature [26,34,48] concerning the mean values, but we introduced the uncertainties with different approaches for the SOC pool kinetics and the humification coefficients. We considered that the decomposition of the two SOC pools was relatively well known since k₁ was estimated from empirical values from the literature and k₂ was estimated from the bare-fallow treatment at Ultuna (see Andrén and Kätterer [26]). We, therefore, utilized normal priors for the kinetic coefficients of the Y and O pools, with a standard deviation of 5% of the initial value.

For the priors regarding the BG humification coefficient (h_r), we were considering the ratio of average contribution of total root-derived C compared to that of AG crop residues (h_s) derived from studies in LTEs, mostly between 10 to 50 years duration [11,14,29,30]. Based on these studies, we considered that the upper limit of the priors for h_r was four times higher than that for h_s and assumed that both had uniform distributions. Similarly, since the previous humification coefficients for amendments (h_a) were derived from inverse modeling, we also assumed a uniform distribution with wider priors, varying around ±100% of the formerly estimated values [29,30]. Priors for green manure and straw were considered the same as for shoots and connected with the same model parameter. The C content of organic matter was assumed between 40% and 51% (Appendix E), with uniform prior, contributing to the posteriors’ uncertainty estimate.

Priors used for the allometric functions

For the Ultuna calibration site, all AG plant biomass is cut at the soil surface and removed from experimental plots, and there is no need to have any priors concerning harvest index data. Thus, the AG crop residues that are not measured consist only of very short stubble parts below the soil surface, including stem bases, and its prior was estimated conservatively as a uniform distribution between 1 and 8% of total AG biomass. This estimate was based on results from an adjacent field trial, where BG stem bases accounted for between 2.5 and 3.8% of shoot biomass at harvest for wheat under different experimental treatments [49].

In the present study, we used allometric functions assuming SR ratios relating linearly to crop yields, one of the most commonly used approaches [6,50], eventually adding a calibrated intercept (see “Testing possible modifications of the ICBM/3B model” section for an explanation of the rationale behind). However, BG estimations are possibly the most uncertain element of any SOC model. Therefore, we focused on estimating this approach’s uncertainty, conditional to the uncertainty in other parameters of the ICBM model. We were basing the priors for SR ratios on literature studies from Europe and North America evaluating SR ratios in field studies at or close to maturity and assuming uniform distributions. For small-grain cereals, it was based on Bolinder et al. [6,51], for maize on Amos and Walters [12] and Bolinder et al. [6,14], while priors for root crops were set accordingly to Bolinder et al. [52]. Priors for oilseed crops were based on the following studies: Ilola et al. [53], Barraclough [54], Pietola and Alakukku [55], Gan et al. [56], and Williams et al. [57].

We used a fixed proportion for rhizodeposition (RD) of 0.65, suggested as an approximation by Bolinder et al. [6], meaning that the plants are delivering RD C to the soil during the growing season in an amount equal to 65% of the total root biomass (i.e. 0.65 × root biomass measured at or close to maturity). This RD consists of root turnover (root hairs and fine roots), cell sloughing of epidermal root tissues, and root secretion of soluble low-molecular organic compounds (exudates) during the growing season. Our coefficient for RD is similar to the average value found for maize of 0.56 by Pausch et al. [58]. Our approximation of 0.65 more or less also equals the mean value (0.66) from the studies on maize when considering only the measurements made about four months after planting in the review by Amos and Walters [12], and it is similar to the value of 0.67 for wheat, maize, and soybean estimated at physiological maturity by Buyanovsky and Wagner [59].

For the model version with an intercept ($\alpha$) in the allometric functions, we assumed for this intercept a uniform prior between 0 and 5 Mg C ha⁻¹; this lies within the range of general estimates for total root biomass [60].

Results

Model calibration

The impact of 20 years of new data on the model

The present model calibration differs from the approach used in the original calibration in 1997 [26], and differences are expected even when considering only the initial time series. These differences could be due to the calibration approach, which relaxes many assumptions made with the sequential calibration approach, but also to differences in the estimations of BG C inputs or the updated climatic functions. In particular, the decomposition of both the old and the young pool are different, with the young pool ($k_1$) decaying slower, and the old pool ($k_2$) decaying faster.

The ICBM mean annual r_clim parameter in Ultuna showed high inter-annual variation, with values ranging from as low as 0.71 to as high as 1.32 (Figure 3). The inter-annual variation followed the same general pattern for all treatments. There were minor differences between the treatments in the average values of r_clim over the whole calibration period (1956 to 2019), from 1.00 for the straw + N reference treatment G (which is by definition scaled to unity) to the highest value of 1.05 for the bare-fallow treatment A. The r_clim parameter tended to be higher when maize was grown (2000–2019), with an average value for all treatments of 1.11, compared to 0.97 during the first period ending in 1999.

Figure 3. The ICBM climatic parameter (r_clim) in Ultuna.

The introduction of new information in the calibration with more than two decades of data slowed down $k_1$ compared to the calibration on the shorter time series and even more when compared with the original calibration from 1997, and in the same way, it sped up $k_2$ (Figure 4). The effect of the new data is particularly evident when comparing the posterior probability distributions of $k_1$ and $k_2.$ When calibrating only on the bare fallow treatment, the differences between prior and posterior probability distributions become not noticeable, also because the uncertainty of the probability distributions increased greatly.

Figure 4. The prior and posterior probability distributions of the model parameters representing the kinetics of the two pools for the standard ICBM/3B model, calibrated on the time series until 1999 and 2019 and on all the plots or only on the bare fallow (BF) plot. Densities are estimated with kernel density estimation, N = 3750.

The model did not seem to require the introduction of an inert pool (Figure 5), which was always calibrated close to zero.

Figure 5. The prior and posterior probability distributions of the intercept ( $\mathbf{\alpha })$ and inert pool parameters other than the kinetic terms for the ICBM/3B model (modified with intercept and inert, respectively) calibrated both on the time series until 1999 and 2019. Densities are estimated with kernel density estimation, N = 3750.

The model residuals in Ultuna for the model without an inert pool and intercept in the allometric functions (Figure 6) slightly increased when calibrating on the longer time series but with no difference in their general pattern.

Figure 6. Average residuals of the simulation in Ultuna considering the standard ICBM/3B version, calibrated both on the time series until 1999 and 2019.

Model residuals (measured minus predicted values) were generally well distributed around the mean (Appendix C), but the model presented some explicit biases toward different treatments, which could be related to processes specific to the treatments (Figure 6). In particular, the model seemed to be overpredicting the calcium nitrate treatment C, and most of the amended treatments not receiving N inputs (B – control, H – Green manure, I – Peat and J – Farmyard manure), except for L (Sawdust), which was still pretty close to neutral, while N (Sawdust + N fertilizer) were also slightly overpredicted. Among the treatments without N addition the bare fallow treatment (A) and straw (F) were instead undepredicted.

The humification coefficients

The posterior mode values of the humification coefficients (Figure 7) for the amendments on the time series to 1999 were highest for peat (0.75), followed by that for sewage sludge and manure, and lowest for sawdust (0.27). The more recent data considering the whole time series to 2019, slowed down the humification coefficients slightly. The ratio of the humification coefficient for roots compared to shoots (i.e. hr/hs) was 2.6 and 4.0 when calibrating on the time series to 1999 and 2019, respectively. This difference was due to a lower value for h_s in 2019 (0.09) than when the model was calibrated on the time series ending 1999 (0.14), whereas h_r remained at 0.36 for both calibration periods. The posterior mode and mean values for all humification coefficients for both time series were similar (Table 2).

Figure 7. The prior and posterior probability distributions of the model parameters other than the kinetic terms for the standard ICBM/3B model calibrated both on the time series until 1999 and 2019. Densities are estimated with kernel density estimation, N = 3750.

Table 2. Results from the ICBM/3B calibration at the Ultuna long-term field experiment without an inert soil organic carbon pool and an intercept for minimal total belowground C inputs for maize.

	Time series 1956–1999			Time series 1956–2019
Parameter	Mode	Mean	95% CI	Mode	Mean	95% CI
k1 (year^−1)	0.36	0.37	0.27;0.45	0.26	0.28	0.24;0.31
k2 (×10^−3year^−1)	10	9.9	9.2;11.2	9.7	9.5	9.3;11.0
hr	0.36	0.31	0.22;0.37	0.36	0.33	0.34;0.37
hs	0.14	0.13	0.09;0.15	0.09	0.09	0.07;0.1
hFYM	0.31	0.31	0.27;0.34	0.26	0.26	0.24;0.28
hPEAT	0.75	0.75	0.71;0.79	0.66	0.65	0.64;0.68
hSAW	0.27	0.27	0.24;0.31	0.23	0.23	0.21;0.25
hSLU	0.57	0.58	0.51;0.62	0.41	0.41	0.37;0.44
S:R cereals	4.81	4.65	3.7;5.49	4.62	4.70	4.18;5.5
S:R root crops	25.9	30.0	2.5;57.8	43.8	33.4	5.99;57.8
S:R oilseed	1.11	7.11	0.76;15.6	7.29	8.14	1.2;15.7
S:R maize	NA	NA	NA	12.1	14.1	9.64;23.1

Note: CI: lower and upper 95% confidence interval; NA: not applicable; FYM: farmyard manure; SAW: sawdust; SLU: sludge; S:R: shoot to root.

The simulation of the bare fallow treatment

All the models considered simulated the Ultuna bare fallow treatment well, but when calibrating only on the bare fallow treatment, the model had a more accentuated initial decrease. This pattern was true for the standard and the model with an inert pool. This different model behavior is mainly due to the combination of a different initialization and a slower decomposition described by the $k_2$ parameter (Figure 4), the most important in the bare fallow case since it describes the decomposition of the older organic matter. The $k_2$ parameter when calibrating only on the bare fallow treatment was slower when calibrating until 1999 and when calibrating over the whole time series and was closer to the original value from 1997 (i.e. 0.0061).

Also in Lanna, although the bare fallow was not simulated well due to processes unaccounted for (most likely inputs above zero since the measurements record an increase in SOC), the different models had almost no difference in the residuals (Figure 8).

Figure 8. The model simulation in Lanna, considering only the standard ICBM/3B model on the time series until 1999 and 2019 (the shaded area represents the uncertainty of the model calibrated until 1999) (panels A, C, E, G, I, J, L, and N) and the distributions of the RMSE in Lanna of all the calibrated model versions (panels B, D, F, H, J, K, M and O). The sewage sludge + metals treatment in Lanna was not considered for validation since it presents conditions that are far from the scope of the ICBM model.

The uncertainties of belowground C inputs estimation

Unsurprisingly there were strong interactions between the terms controlling the C inputs, particularly BG, and the humification factors. There was a relatively thin but elongated space of combinations of S:R and $h_r$ values for cereals with relatively good fitness (Figure 9) and a broader range for maize (Figure 10). In both cases, there was a relatively vast region of high probability density with multiple local optima and little difference in the fitness of many different parameters’ combinations. This explains the broad posterior distributions of the S:R parameters (Figure 7).

Figure 9. The values of cereals shoot:root ratio and the humification coefficient for root for each element of the MCMC chains (after random resampling) for the standard ICBM/3B model calibrated on the time series until 2019 showing the equifinality between the two parameters.

Figure 10. The values of maize shoot:root ratio and the humification coefficient for root for each element of the MCMC chains (after random resampling) for the standard ICBM/3B model calibrated on the time series until 2019 showing the equifinality between the two parameters.

Model predictions and validation

The model represented, in general, the measured SOC in Ultuna relatively well, with an average RMSE of 3 Mg ha⁻¹ when calibrated on the whole time series and 2.6 Mg ha⁻¹ when calibrated for data until 1999. The less extreme treatments, as well as the treatments losing C, are represented particularly well by the model, while it struggles more with the treatments receiving amendments that lead to SOC accumulation (Figure 11). The residuals were evenly distributed when considering the whole calibration, but biases were specific to each treatment (Figure 6). The last decade presented a regular bias in many treatments (Appendix C).

Figure 11. The model simulation in Ultuna considering all three model versions (all calibrated on the time series until 2019).

Interestingly, the inert pool did not seem needed regardless of how long the calibration time series was, with its probability distribution resulting in both calibrations being very small and skewed toward zero. The intercept in the allometric functions seemed helpful to the model only to some extent but was still skewed toward zero (Figure 5).

When considering the error of model predictions in the independent validation on the Lanna site (Figure 8), the intercept improved the model’s performance only in the two mineral fertilization treatments. The minimal impact of the intercept on residuals was even more evident when the model was calibrated on the shorter time series.

The standard ICBM/3B model had mixed results in Lanna, representing relatively well most treatments but struggling in particular with the bare fallow and the sludge treatments (Figure 8).

The four European long-term bare fallow sites considered were instead represented all pretty well except Kursk (Figure 12), but the initialization had to be tuned for Versailles, with a larger initial proportion of young material at the start of the experiment to represent the land use change (from grassland to cropland). There was no appreciable difference between the model calibrations on the two time series.

Figure 12. The simulations from the standard ICBM/3B model, calibrated on the time series until 1999 and 2019, on four long-term bare fallow experiments in Europe [ 33]. Please notice for Versailles an additional simulation with a different initialization ratio (0.8) other than the one resulting from the calibration (0.952).

Discussion

Comparing the new model calibration with the former version of ICBM

The new model calibration presents some differences compared to the original model from 1997 that deserves some discussion. The original calibration used a stepwise parameterization process to derive r_clim [26]. Briefly, this was assuming h was constant (i.e. h_s and h_r were the same) and determined on only four of the treatments at Ultuna (B, C, F, and G), where a linear regression against yield differences provided the diverse r_clim values under the assumption that higher plant biomass gives higher transpiration and thereby lower soil water content and r_clim values. The principle used in our calibration estimating r_clim values based on climatic response functions (i.e. r_temp and r_wat) was introduced later [27]. Using the climatic response functions showed that the bare-fallow in Ultuna still had the highest r_clim value among treatments but did not provide as high a difference in r_clim between the treatments as the former approach (Table 1 in Andrén and Kätterer [26]). This is partly explained by the fact that all treatments constrained our Bayesian calibration and allowed the other parameters (k- and h-values and allometric functions) to explain the variance between treatments. In particular, a faster decay rate of the old SOC pool (k₂) compensated the lower r_clim value compared to the original calibration was compensated. The slightly lower r_clim value for maize compared to other crops is mainly due to the shading of the maize plants.

Estimating the climatic parameter with climatic response functions has the advantage that this approach is more easily generalizable. For example, applied to more contrasting climatic conditions, such as the eight agricultural regions within the Swedish national inventory reporting system [40], it provided r_clim values ranging from 0.67 in the north to 1.30 in southern Sweden. Similarly, the approach has also helped describe climate’s influence on SOC decomposition at much broader scales [35,61], and the approach can be easily tuned to these broader ranges.

Another issue to consider is the equifinality that this class of models (first-order SOC models, driven by C inputs) presents [19,62]. The optima for the related parameters are, therefore, not well constrained. Because of that, slight differences in methods or assumptions can shift the optimal calibration values. The calibration approach adopted in the present study differs from the former version since it does not consider optima for single treatments to constrain parameters in sequence. Instead, in the current calibration approach, since all parameters are calibrated simultaneously on all treatments, the parameter values will represent optima between all the data points. With the former sequential approach, the SOC dynamics in the bare fallow treatment was assumed to be the most reasonable approximation of the old SOC kinetic coefficient, which was then used to define $k_2$ in all the other treatments (thus making one more assumption and reducing the amount of possible optimal parameter sets). The overall optima for the SOC decomposition rates could differ from the bare fallow optimum for several reasons, including substrate interactions such as priming effects), differences in microbial community composition, or effects of soil structure on biogeochemical processes.

Testing modifications of the ICBM/3B model

When considering the residuals of the three model versions, the only appreciable differences when introducing an intercept in the allometric functions occurred for the mineral fertilizer-only treatments (Figure 11). The addition of an inert pool seemed to make very little difference in all treatments, confirmed by the small size of the calibrations assigned to the inert pool (Figure 5). Besides, guessing the size of the “inert” pool is difficult since it represents a coarse approximation as a concept [41] and is connected to the slower cycling part of SOC, which is likely dependent on local edaphic factors [63,64].

The standard ICBM/3B model, similar to the original ICBM, performed relatively well in Lanna concerning most treatments but failed to represent the high C input extremes (Figure 8). In particular, the model struggled to represent sludges and, to a lesser extent, compost. The model assumes that all amendments have the same kinetic coefficient (i.e. k₁), which remains a simplification. Furthermore, concerning sludge, its properties depend on the process used in the sewage treatment plants, which can be variable over time and from one plant to another, depending on the technical solutions adopted. In the lowest C inputs extreme (i.e. bare fallow), it is possible that the experiment at Lanna had C inputs from weeds since weeding was mechanical, contrary to the Ultuna bare-fallow plots, where weeding is done by hand and more frequently.

Testing the model representation for the decay of the old SOC pool

Despite the model’s $k_2$ being optimized on all treatments rather than just on the Ultuna bare-fallow as in the 1997 model version, the model performed well (Figure 12) in other long-term bare-fallow experiments [33]. The model’s excellent performance in these experiments is not extremely surprising since the kinetics of old organic matter depends mainly on time, temperature, and moisture and very little on other edaphic factors [41,65]. Nevertheless, it represents a good validation of the general robustness of the kinetic coefficient of the old pool. In the absence of C inputs, the young pool rapidly more or less disappears, and the model becomes quickly equivalent to a one-compartment model represented by the old pool. Consequently, given the broad climatic gradient of the bare-fallow experiments, it is also a good validation of the climatic module since r_clim is the only driver under such conditions. The model worked well for Grignon and Askov, while it struggled more for Versailles and Kursk. The Kursk site is a chernozem, a particularly rich soil where local processes have accumulated organic matter. These processes are probably not fully captured by our model. The site at Versaille is older than the other sites and has a different history. This site was a former grassland [33], and the higher biomass productivity or at least root C inputs of this ecosystem, compared to cropland, seems to have influenced the quality of the initial SOC [41]. When we tuned the model with an 80% initial size of the old pool (and therefore a 20% young pool), compared with the 98% old pool assumed by default in the other sites (i.e. set to the same initial proportion than Ultuna), the model represented the site in Versailles closely. This suggests that the only discrepancy of the Versailles site is not in the decomposition rate of the old pool but probably in the higher amount of young material at the beginning of the experiment, which can be represented by site-specific model starting conditions.

Model uncertainty

Model predictions residuals

Despite a relatively high uncertainty of model parameters (Figure 4, Figure 7), model prediction boundaries were well constrained (Figure 8, Figure 11). Such relatively narrow predictions’ uncertainty suggests that the model structure makes it robust and suitable for describing the SOC balance even if parameters are not precisely known since predictions are relatively insensitive to parameter uncertainty.

Most of the treatments not receiving N were overpredicted except for sawdust, calcium nitrate and bare fallow. All other underpredicted treatments receive some form of quickly available N. This pattern might be related to N-related processes but not to the C/N ratio (which does not correlate with the model residual). More likely, a more readily available N form might slightly slow down the SOC decomposition in a way not considered by the model by interacting with the soil microbiota. This effect is clear and well known for N-limited environments such as forests [66], and it is mentioned in the literature for agricultural systems, although less clear [67]. The introduction of maize seemed to increase the difficulties for the model, probably because of additional processes not captured by it, for example the already mentioned possible shading effect, or other similar effect on the micrometeorology of the plots, or effects related with inputs estimation (maize inputs are estimated lower than previous crops).

While in temperate agricultural systems [68] the C:N ratio seems to drive decomposition with no effect of N fertilizers, N addition plays a more important role in determining SOC decomposition in poorer systems [69]. In our case, long-term treatments without additional N added might have little bioavailable N compared with more conventional agricultural treatments and present some N limitations. Regarding the absolute error, the model struggled with treatments receiving amendments where the material is already partially decomposed and can present some very specific chemical quality, such as peat. The model considers the specificities of these amendments with different humification coefficients, but the amendments might also lead to differences in the chemistry of the remaining SOC (conceptualized by the O pool), which may explain the error.

Model equifinality and belowground C inputs

The strong interactions between the mass of C inputs and SOC decomposition are well recognized for this particular class of SOC model [3,19,70,71]. Decomposition rates are uncertain due to limitations in the measurement techniques on the relevant time scales. These results also confirm the high uncertainty in the estimation of C inputs, in particular due to the limitations when measuring BG biomass and fluxes, which constitutes a major limitation in SOC modelling [4]. These uncertainties affect all SOC models and are a reason to keep the model as simple as possible since a more detailed model (with more pools) would face the same problem but only distribute the uncertainty among more parameters.

The equifinality caused by the interactions between the BG C input estimates and the kinetics of both SOC pools make the humification coefficient of roots and the SR ratios uncertain. The broad spread of the various chains in the traceplots of the relative parameter (Appendix D) also seem to be a consequence of such equifinality, which becomes visible when plotting the values of the cereals’ SR ratio and the $h_r$ coefficient for each element of the MCMC chains (Figure 9 and Figure 10). This interaction is particularly important for small-grain cereals and maize, the dominant crops in the time series at Ultuna, representing 53 and 39%, respectively.

The mode of posteriors for the humification coefficients for crop residues at Ultuna (Figure 7 and Table 2) were within the range of those estimated in LTEs involving mostly small-grain cereals and maize, with values between 0.08 to 0.20 and 0.14 to 0.44 for AG and BG crop residues, respectively [10,29,30,72–75]. The ratio we obtained for the humification coefficient for roots compared to shoots (i.e. hr/hs) was also within the range of previously estimated values (i.e. 1.3 to 3.7) obtained in LTEs and reviews [11,14,23,29,30]. Compared to similar studies using the RothC [24] and AMG [25] models, our calibration highlights the usefulness of attributing a higher C contribution of root C to the old SOC pool also in the ICBM model. For the organic amendments, humification coefficients were similar to those calculated earlier by Kätterer et al. [29] using a first-order kinetics approach and linear approximation. For manure, h_FYM calculated on the time series ending 1999 is almost identical to that used in the Hénin and Dupuis [20] model (i.e. 0.30), whereas h_FYM derived using the whole time series is similar to the average value (i.e. 0.26) in a recent review on data from Canadian LTEs [76].

The mode of posterior SR ratios for the dominant crops, small-grain cereals and maize, were 4.8 and 4.6 for small-grain cereals when calibrated on the time series up to 1999 and 2019, respectively, and 12.1 for maize (Figure 7 and Table 2). This is close to the average SR ratio that was found for small-grain cereals (i.e. 5.6) in a review study on Scandinavian data [77]. The relatively high SR ratio for maize is probably because its mean growth period at this site (see “The climatic functions” section) was only 95 ± 10 days, whereas the potential growing season in the area is about 140 days. However, the SR ratio from our calibration is similar to that reviewed for maize by Amos and Walters [12], where SR ratios were between 11.1 and 14.5 for the measurements made 73 to 105 days after planting. Other recent studies on varieties used in Europe are also reporting SR ratios for maize well above 10.0 [13,78].

The impact of almost three decades of data on the calibration

The modifications in the two kinetic coefficients caused by two more decades of data align with the expectations since decomposition naturally slows down with time. Since the decomposition, initially well approximated by a single exponential decaying pool [33], departs from such simple decomposition kinetics the more time passes, and in theory, the longer the time scale, the more pools would be useful. Although there is no real “inert” material in the soil [41], the decay tends to slow down over time more than what is described by a single exponential law because the SOC quality changes gradually over time toward more recalcitrant material. This slowing down in the decomposition of the old SOC pool is why we introduced in this calibration the possibility of explaining part of the variance with a constant inert pool, a relatively simple modification that would allow the model to better represent this decrease in decomposition over time. However, this inert pool was not particularly crucial (Figure 8), suggesting that a two-pool model can still capture most of the variance of SOC over time, even approaching the scale of a century. Residuals also were not particularly high, although they did increase from an RMSE of 2.6 to 3.0 Mg ha⁻¹. Adding the last decades of data shifted the optimum toward reduced humification coefficients and slightly faster old pool decomposition following the slightly diminished C inputs from maize (Appendix F). Specifically, the average root C inputs, representing the largest organic contribution other than amendments, evolved from an average of roughly 0.6 to 0.32 Mg C ha⁻¹ y⁻¹.

The regular bias present in many of the treatments for the points of the last decade (Appendix C) might also be related to some other effect connected with maize. For example, since the experimental plots are adjacent and relatively small, in some treatments, plants might be shaded by the plants in another (maize grows much taller than all other crops previously grown in the experiment) in a more favorable position in the randomized block design.

When calibrating the model only on the bare-fallow treatments, the effect of the last two decades of data becomes clear (Figure 4). The term $k_1$ is mostly not constrained by the data from the bare-fallow treatment since, with no inputs, the young material disappears quickly from the system, but the term $k_2$ is definitely affected by the new data. It is not surprising that the decomposition of SOC in the bare-fallow seems to slow down. Still, given how well the standard ICBM/3B model calibrated on all treatments performed on all long-term bare-fallow experiments, this effect is small.

Description of the material attached to the study

We included as Supplementary Material:

1. The updated functions for calculating the climatic scaling are downloadable at https://github.com/ilmenichetti/reclim and possible to install in R through the package devtools.

2. A table describing all the applications of ICBM in the literature that we are aware of, with references

Conclusions

The ICBM model was calibrated with new data and considering a new structure (ICBM/3B), which allowed us to test if the model assumptions still hold on a longer time scale than the original model calibration. The assumptions still hold; in particular, the model does not need to consider more than two pools to represent the observed changes in SOC stocks. Introducing an intercept in the allometric functions was also not useful.

The model was calibrated with a multi-objective Bayesian approach, resulting in slightly different kinetic coefficients than the former sequential calibration approach, with a slightly faster decomposition for the old pool that nevertheless seems realistic. In particular, the calibration made evident the strong equifinality of the model due to the interactions between the decomposition of SOC and uncertain C inputs, which is common for compartmental SOC models. The lack of precise knowledge about BG C inputs is confirmed as the highest uncertainty for this class of SOC models. At the same time, the calibration pointed out the relative robustness of model predictions to model equifinality and uncertainty of BG C input estimates.

Keeping the model as simple as possible seems an effective way to deal with such equifinality since added complexity would only increase it without substantially increasing the model performance.

The ICBM model is one of the simplest compartmental SOC models available, and, nevertheless, it captures most of the SOC variance over time, making it useful for applied studies involving SOC decomposition. This simplicity is one of the main assets of the model. There is still a need to evaluate model performance with new experimental results and measurements of SOC dynamics. In particular, future work will focus on a refined regional multi-site calibration using data from other LTEs to test further and improve the generalizability of modeling results.

Acknowledgments

We are thankful to our Faculty for providing resources to maintain the Swedish long-term experiments and especially to Per Hillström, who carefully managed the Ultuna frame trial for the three recent decades, and to Karl-Gunnar Söderdahl, who took over the management in recent years. We are also grateful to the persons and institutions who provided data from the long-term bare fallow network sites (see Barré et al. [33]). In particular, we acknowledge the inputs from Vladimir Romanenkov (Kursk site, Lomonosov Moscow State University, All-Russian Institute for Agrochemistry), Folkert van Oort and Sabine Houot (the two French sites, French National Institute for Agricultural Research and Agroparistech), and Bent T. Christensen (Askov Experimental Station, Dept. Agroecology, Aarhus University).

Disclosure statement

During the preparation of this work, the author(s) used Grammarly in order to improve the readability of the manuscript. After using this tool/service, the author(s) reviewed and edited the content as needed and take(s) full responsibility for the content of the publication.

Data availability statement

The data that support the findings of this study are available on request from the corresponding author, LM. The data are not publicly available due to being property of the author’s institution, which requires following particular policies for disclosing them.

Additional information

Funding

This work was funded by Swedish Government Research Council for Sustainable Development (FORMAS, grant number 2022-00214) within the project “Carbon sequestration in Swedish cropland soils” and the Swedish Farmers’ Foundation for Agricultural Research (SLF, grant number O-18-23-141). The SMED (Swedish Environmental Emissions Data) consortium also provided funding for this research.

Appendix A

The model has been adapted from Kätterer and Andrén [34]. In order to have a climatic parameter $r_{\mathit{\text{clim}}}$ variable over time, we started from the step version of the ICBM models.

The initial model: $$Y_{t+1}=\left(Y_t+I_t\right)\bullet e^{-k_y\bullet r_t}$$ $$O_{t+i}=\left(O_t-h\frac{k_y\bullet \left(Y_t+I_t\right)}{k_o-k_y}\right)\bullet e^{-k_o\bullet r_t}+h\frac{k_y\bullet \left(Y_t+I_t\right)}{k_o-k_y}\bullet e^{-k_y\bullet r_t}$$

For simplicity, we substitute with one single term for the flux between the pools: $$\phi =h\frac{k_y\bullet \left(Y_t+I_t\right)}{k_o-k_y}$$

Obtaining a more compact writing: $$O_{t+i}=\left(O_t-\phi \right)\bullet e^{-k_o\bullet r_t}+\phi \bullet e^{-k_y\bullet r_t}$$

To consider the three input classes, we now introduce three Young pools with three humifications (but with the same kinetic $k_y$), R (Roots), S (Shoots), and M (Manure): $${Y^R}_{t+1}=\left({Y^R}_t+{I^R}_t\right)\bullet e^{-k_y\bullet r_t}$$ $${Y^S}_{t+1}=\left({Y^S}_t+{I^S}_t\right)\bullet e^{-k_y\bullet r_t}$$ $${Y^M}_{t+1}=\left({Y^M}_t+{I^M}_t\right)\bullet e^{-k_y\bullet r_t}$$

We therefore have also three fluxes: $${\phi }^R=h^R\bullet \frac{k_y\bullet \left({Y^R}_t+{I^R}_t\right)}{k_o-k_y}$$ $${\phi }^S=h^s\bullet \frac{k_y\bullet \left({Y^S}_t+{I^S}_t\right)}{k_o-k_y}$$ $${\phi }^M=h^M\bullet \frac{k_y\bullet \left({Y^M}_t+{I^M}_t\right)}{k_o-k_y}$$

Which all end up in the Old pool: $$O_{t+i}=\left(O_t-{\phi }^R-{\phi }^S-{\phi }^M\right)\bullet e^{-k_o\bullet R_t}+\left({\phi }^R+{\phi }^S+{\phi }^M\right)\bullet e^{-k_y\bullet R_t}$$

The total SOC is given by the sum of the Old and Young pools: $${{\mathit{\text{SOC}}}_t=O}_t+Y_t$$

To which we can add a constant Inert pool ($G$): $${{\mathit{\text{SOC}}}_t=O}_t+Y_t+G$$

Appendix B

A diagram of the information flow in the climatic module.

Appendix C

All the residuals of the two calibrations. Empty points represent the calibration until 1999, and solid fill points until 2019.

Appendix D

The traceplots of the MCMCs of selected parameters.

Appendix E

The prior and posterior probability distribution of the C content of organic matter.

Appendix F

The C inputs from roots estimated by the Bayesian calibration.