A Bayesian modeling framework for estimating equilibrium soil organic C sequestration in agroforestry systems

Agroforestry is a form of productive land management that combines trees or bushes with annual crops or pasture, and it can bring benefits in terms of food security and increased carbon (C) sequestration compared with conventional agriculture. But agroforestry as a structured form of agronomic management is relatively new compared with well-established and widespread agronomic systems. Consequently, there is a lack of data and few models of soil organic carbon (SOC) have been developed specifically for agroforestry systems. Also, agroforestry SOC sequestration data measured in field experiments are often reported only as average linear sequestration rates over the study period. This approach, equivalent to zero-order kinetics, makes it difficult to compare results since, in reality, SOC sequestration rates are variable over time and change depending on the duration of measurements. Sequestration rates are also strongly dependent on former C stocks in the soil, further hampering comparisons between agroforestry systems established on different former land uses. To describe the SOC stocks variation over time, researchers often employ models considering at least firstorder kinetics. This approach can take care of the two above mentioned issues, considering both the variation of the sequestration over time and the effect of previous land use. However, the variability of agroforestry systems makes applying these models more challenging compared to simpler agricultural systems. To deal with this problem we propose to use detailed uncertainty estimation methods, based on stochastic calibrations that can deal with broad probability distributions. To do so, we adapted a first-order compartmental SOC model to agroforestry systems. It was calibrated within a Bayesian framework on global agroforestry data. Compared to linear coefficients, the model (ICBMAgroforestry) estimates equilibrium SOC stocks of different agroforestry systems probabilistically and is providing uncertainty bounds. These values are independent of initial land use and time duration of the experiments. ICBMAgroforestry can be used for rapid assessment and comparison of the maximum potential SOC stocks for different agroforestry systems and climatic zones. In this study, we could use our approach to estimate the global maximum C that can be sequestered by agroforestry systems at equilibrium, which ranged between 156 and 263 Mg C ha−1 on average, above but comparable with similar estimates for simpler agricultural systems.


Introduction
With an impending climate crisis, the world is actively seeking new methods to limit climate change and reduce the vulnerabilities of ecosystem services on which humans depend (Sterner et al., 2019). One promising method is agroforestry, which combines perennial woody plants with more conventional annual crops or pastures in a productive system. Agroforestry is an evolving concept used in many different environments to improve resource use efficiency and the resilience of traditional agricultural systems (Mbow et al., 2014;Weiwei et al., 2015). By combining plants that occupy different ecological niches, agroforestry production systems achieve the dual benefits of a) higher efficiency, due to better niche exploitation and synergies between species, and b) diversification of niches, making the system more resilient to extreme events (Lasco et al., 2014). Thus agroforestry can be an important practice in climate change adaptation (Verchot et al., 2007), particularly in regions where agricultural inputs are low and exposure to climate extremes is high.
On top of that, the niche stratification in agroforestry systems often increases carbon (C) inputs compared with conventional agricultural systems. Agroforestry systems can represent an improvement compared with agriculture in terms of soil organic carbon (SOC) stocks (De Stefano and Jacobson, 2018). The carbon that is absorbed from the atmosphere through photosynthesis enters the soil, and when the C inputs are bigger than the C output due to microbial mineralization we have a net negative C balance and a climate change mitigation effect (Corbeels et al., 2018;Lorenz and Lal, 2014;Ramachandran Nair et al., 2009). Estimating the magnitude of the C storage effect is essential for promoting C sequestration in agroforestry systems. Agroforestry systems can be quite diverse, as can their C inputs, which makes it difficult to estimate their impact on SOC, particularly on a large scale. Nevertheless, recent studies have attempted to assess the potential of agroforestry systems for C sequestration by aggregating data available from the literature Feliciano et al., 2018;De Stefano and Jacobson, 2018;Zomer et al., 2016). Such studies identified the great potential of these systems, but accurate estimation and reporting of C sequestration in agroforestry is still a barrier preventing its adoption for climate change mitigation efforts (Rosenstock et al., 2018).
There are a few studies over relevant time scales and most studies report results as the difference in C stocks before and after implementation of an agroforestry system (either with a synchronic or, more commonly, a less robust diachronic setup). This difference is usually divided by the number of years since the start of the study, giving an average "sequestration rate" usually expressed in Mg C ha −1 y −1 . This approach is used regardless of the initial state of the system and does not differentiate between initial land-use. However, since sequestration rates are, in reality, highly variable over time and space, the applicability of these average rates is limited to time scales close to the number of years over which they were measured. Field trials on agroforestry are relatively new, and few cover more than a couple of decades, with most lasting only 5-20 years Corbeels et al., 2018;Feliciano et al., 2018), whereas the scale for effective climate change mitigation policies is at least 20-50 years (European Commission, 2013).
There are more detailed approaches to estimate SOC sequestration. The conventional first-order SOC decay theory can easily describe most of the variation in soil C sequestration rates over time (Coleman et al., 1997;Andren and Kätterer, 1997;Powlson et al., 1998). This theory has been applied and tested for decades and can capture most of the observed variation in long-term agricultural field experiments (Cagnarini et al., 2019;Dechow et al., 2019;Franko and Merbach, 2017;Kröbel et al., 2016;Menichetti et al., 2019). Some other sources of variation exist, leading to higher-order kinetics, and the scientific community is actively experimenting with these (Moyano et al., 2018). However, soil C models are particularly sensitive to C inputs (Begum et al., 2017). These are already relatively uncertain in conventional agricultural systems, but a lot more in the case of agroforestry systems due to the multitude of functionally different plant species and their ecological interactions as well as differences in litter quality. Therefore, more accurate higher-order kinetics would not be an advantage in the case considered here, since the small improvements introduced by the second or higher order would likely not be detectable with the data available. The suggested new steady-state approach for Tier 2 and the Tier 3 SOC stock calculations currently used within IPCC are also based on first-order kinetics (Jim et al., 2006;Ogle et al., 2019). The main difference with the model proposed here is in the statistical approach considering the uncertainty of data in the model predictions.
Assuming an average sequestration rate over time, as suggested for Tier 1 calculations Jim et al., 2006), is instead equivalent to assuming zero-order decay kinetics, and this may severely overestimate C stocks when projected over long time scales since it would overestimate future sequestration rates. Furthermore, the zeroorder assumption is not accounting for the fact that SOC stocks converge towards an equilibrium depending on the balance between C input and output rates. Following a land-use change to agroforestry, the system will move towards a new, either higher or lower SOC equilibrium depending on the SOC stocks under the former land use. This new equilibrium is determined by the mass and quality of the new C input rates, and the rate of change is determined by difference between the initial state and the new steady state C stock. Therefore, SOC sequestration rates depend on previous land use (De Stefano and Jacobson, 2018).
By estimating the new SOC stocks at equilibrium with first-order kinetics, it is possible to obtain a value independent from past land use that also accounts for the variation of the sequestration rate over time. However, the main difficulty encountered in any attempt to estimate this equilibrium state for agroforestry systems is that such systems are remarkably diverse. This diversity causes a considerable variation in both input rates of different litter quality and recorded SOC sequestration rates (Corbeels et al., 2018;Feliciano et al., 2018). This variation is difficult to handle and represents an additional obstacle to the application of more refined approaches to estimate changes in C stocks in agroforestry systems.
In this study, we applied a Bayesian approach that integrates parameters and variables uncertainties stochastically. Bayesian calibrations deal with realistic probability distributions as input parameters, and combine and maintain such information in the model output without having to approximate it with given probability functions. Another advantage is that this approach makes it possible to fully utilize as input all the information in highly variable datasets. In order to keep model complexity low and focus on the statistical approach, we utilized the Introductory Carbon Balance Model (ICBM) (Andren and Kätterer, 1997), a twocompartment first-order SOC model often selected for its simplicity and ease of handling (Wutzler and Reichstein, 2013). The model is also particularly flexible in its climate scaling modules.
The aims of this study were: 1) to calibrate a version of ICBM specifically modified for agroforestry in a Bayesian framework on a global dataset of agroforestry experiments and 2) using the model for estimating the maximum SOC stocks (at equilibrium) of several agroforestry systems.
The calibrated version of the model (called ICBM Agroforestry ), including the result of calibrations as priors for updating the present calibration with local data, is available as supporting material.

Model structure description
The model relies on two separate decaying SOC pools, young (fast) and old (slow) (Andren and Kätterer, 1997). In our conceptualization, we considered SOC in all soil particles < 2 mm. The original model version includes two C input classes: plant material (above-and belowground combined) and organic amendments. We expanded them to four classes considering separately: woody roots from trees, fine roots from annual crops and trees, and above-ground litter from trees and annual crops. The model considers annual time steps. To describe agroforestry systems, we extended the fast pool to four different pools humifying in parallel: woody roots, fine roots, amendments, and litter from aboveground biomass (Fig. 1). Each of these pools humifies with different kinetics, but they all enter the slower pool. Then we updated the model parameters using a global database of agroforestry SOC data (Feliciano et al., 2018).
The mathematics of the model is relatively simple and based on linear derivatives. However, the approach still manages to capture most of the variance observed in SOC decay over several long-term experiments (Andren and Kätterer, 1997;Kätterer and Andrén, 2001;Kröbel et al., 2016). In the original formulation of ICBM, C inputs (I) are added to a decaying 'young' pool (Y). This pool decays with constant k y , according to the following equation defining its dynamics: The local external kinetic modifiers are all summarized in the climate scaling factor r e , which can vary over time as a function of measured forcing variables (such as air temperature and precipitation). In some cases, r e can also include edaphic modifiers (but not in this case, due to lack of data). The term r e is calculated based on time series and, ideally, the longer the time series the more it will be representative of the site-specific climate. The decay constant (k o ) of the 'old' pool (represented by 'O') is estimated from a soil under bare fallow in a longterm (> 50 years) field experiment that has not received any inputs since its establishment (Andren and Kätterer, 1997), and therefore represents only the decay of SOC. Another parameter in the model that is crucial to determining the C stocks is the proportion of material that is 'humified', or moved from the young to the old pool. This fraction is denoted h and called humification coefficient and it follows approximately the early definition made by Hénin and Dupuis (1945). Estimated values are available for several organic materials, based on longterm data series (e.g., Kätterer et al., 2011). The resulting function, describing the state of the old pool, is: The model can be solved analytically, thus defining the steady state (SS) of the two pools (Kätterer and Andrén, 2001) mathematically: where Y SS and O SS represent the steady-state C content of the young and old pool, respectively. The term I represents annual C inputs, which are assumed to be constant (an assumption that is necessary to postulate steady-state) and can be approximated by the average annual inputs. Due to the requirements of the software used for calibration, we wrote the equations in discrete time step form, starting from the model described above: where the C in the young pool at time t + 1 (Y t+1 ) corresponds to the C in the young pool at time t (Y t ) plus the inputs at the same time step (I t ), which is subject to exponential decay. The term r e,t can vary at every time step and is driven by climate time series (in this case, soil temperature and moisture). The decay is defined by the decay constant k y and by the climatic variable (temperature and moisture) as a function of time, summarized by r e,t . In the present application, since we did not have specific data on the variations in weather over time, the term r e,t was assumed to be constant over time (r e ). The C in the old pool at time t + 1 (O t+1 ) is defined in a similar way (O t ), but with inputs from the humification flux from each young pool : and: To consider several input classes, we introduced four different young pools (Y i ), where i represents the inputs origin classes fine roots including rhizodeposition, woody coarse roots, shoots, and organic amendments (i = r, w, s and a, respectively): n r w s a n t n t k r 1 , , , where n defines the four different input classes (i). All young pools are then summed together. The same is done with the fluxes according to Eq. (6). For simplicity, the same kinetic parameter (k y ) is assumed for all non-woody materials (r, s and a), but not the decay of woody roots (k y,w ). This is supposed to be slower due to higher chemical recalcitrance and time needed for decomposer organisms to colonize the whole root after plant death. Rates k y and k y,w have the same prior probability distribution but are calibrated independently. For more details about the mathematical foundations of the model, we refer the reader to former publications (Andren and Kätterer, 1997;Andrén et al., 2004;Kätterer and Andrén, 2001;Kröbel et al., 2016).
The steady-state (SS) equations that express the total SOC stocks at equilibrium then become:

Fig. 1.
A schema of the model (shoots denotes the aboveground biomass production, roots the belowground biomass production).
L. Menichetti, et al. Agriculture, Ecosystems and Environment 303 (2020) The total SOC stocks at equilibrium are the sum of these three states (Y SS , r+a+s , Y SS,w , and O SS ).
Because of the nature of the steady-state equations, steady states absolute values are determined mostly by the inputs and humification coefficients (at the numerator in Eqs. (9)-(11)) and the decay kinetics. However, since the decay rate of the young pool is always much greater than that of the old pool, the latter is far larger. The steady states are, therefore, mainly influenced by the humification factors and old pool decay kinetic (modifiable by management only indirectly through influencing moisture or temperature) and by the inputs (which can be directly modified by management).

Agroforestry classes by climate region
We used a data from Feliciano et al. (2018) for calibrating the model. The classes to categorize the different agroforestry systems were also from this database, as follows: 1) agrosilvicultural systems, 2) cocoa-based systems, 3) coffee-based systems, 4) systems with fast growing trees, 5) home gardens, 6) systems with intercropping, 7) parklands, 8) silvopastoral systems, and 9) woodlands. These classes are further differentiated by geographical area into 1) arid, 2) semiarid, 3) tropical, and 4) temperate. Based on these classes, we combined the data from different experiments on belowground and aboveground production. After filtering this database for complete time series, we found sufficient data from 64 experimental sites covering the following system × climate categories: agrosilvicultural-arid, agrosilviculturaltropical, home garden-tropical, shadow system-tropical, silvopastoraltemperate, silvopastoral-tropical, woodlots-semiarid and woodlots-tropical (Fig. 2). Feliciano et al. (2018) report SOC to 30 cm depth, which is, therefore, the depth considered in this study.

Climate scaling factor
The climate scaling factor r e r e has been estimated over the year with various approaches, all revolving around an exponential temperature and moisture response function that are combined by multiplication (Andrén et al., 2008). In our application we used another moisture function able to represent the optimum reached before full water saturation, and the subsequent decrease in microbial activity due to anoxic conditions until full saturation (see Moyano et al., 2011 for details). In the present study, both functions were calculated based on global raster data, where data on temperature, solar radiation, and wind speed were obtained from WorldClim v2.0 (Hijmans et al., 2005), while data on moisture was obtained from the ESA CCI dataset (https:// www.esa-soilmoisture-cci.org/). Soil moisture values were averaged over one year of sampling. Although local anomalies might arise when using the high-resolution map, this is not a problem here considering the level of aggregation used in this study. The moisture data are based on a combination of active scatterometer and passive radiometer measurements, and might present some uncertainties due to disturbances added by vegetation (see Dorigo et al., 2017 for a more detailed discussion). As in the original ICBM (Andren and Kätterer, 1997), the moisture and temperature functions are normalized to 1 for Ultuna (Sweden), the long-term experimental site used for calibrating the kinetics of the two SOC pools. To include random spatial error, we rescaled the data on the average of the whole of Upplands County, where which Ultuna lies. The raster containing the rescaled climate factor r e (Fig. 3) was utilized together with a global Köppen climate map (Peel et al., 2007), and then reclassified (Fig. 4) to calculate average r e for each climate class considered in the calibration. Since data on exact site coordinates were lacking, we performed model calibration using discrete climatic values averaged for each class. A georeferenced raster is available for model users having known site coordinates.

Calibration and development of a predictive tool
We applied linear partitioning of the belowground flux in woody roots and rhizodeposits for each agroforestry class in three steps based on aboveground biomass.
Step 1) a tree:crop biomass ratio to determine the amount of shoot (aboveground) biomass in trees and the amount in other plants.
Step 2) a root:shoot ratio to determine how much C is allocated to the soil as roots.
Step 3) a root:rhizodeposition factor to calculate the amount of rhizodeposition produced by plants and trees. All these factors are calibrated (and resulting posterior distributions are included as supplementary electronic material). Roots of trees are also characterized by a mortality rate, which defines when that C first enters SOC (the young pool, i.e., particles < 2 mm), after a specific time representing aggregated tree mortality. This value is an approximation and includes implicitly the effect of coarse roots entering the soil and that of fine root mortality and rhizodeposition. Since tree fine root mortality is extremely uncertain according to the literature, we decided to calibrate it implicitly in this work, letting the calibration take care of the conceptual error by compensating for it. Crop fine roots were instead assumed to enter the soil every year. Aboveground inputs were calculated based on an empirical formula (Kumar, 2010). To keep the model structure as simple as possible, the change in total aboveground biomass over time was simulated by assuming constant biomass increase until a maximum, which thereafter remains constant. This produces a piecewise-defined function, which is a reasonable approximation of more complex saturation functions (such as a Michaelis-Menten).
The prior distributions for the decay constants of the young pools (k y and k y,m ) and for the old pool (k o ), and the four humification coefficients (h r , h w , h s , and h o ) were considered to be normally distributed since they have already been calibrated and validated extensively in the literature (Andrén et al., 2004;Kätterer and Andrén, 2001). In order to maintain compatibility with former calibrations for L. Menichetti, et al. Agriculture, Ecosystems and Environment 303 (2020) 107118 all these parameters, we selected normal distributions centered in the previous calibrated value, assuming an error of ± 0.05 and truncated in a range of ± 0.1 with respect to the former calibrated value. We selected as starting average value for h 0.35 for roots, woody roots and amendments, 0.15 for aboveground material, 0.8 for k y and k y,w and 0.00605 for k O . In this study, the large database makes the calibration relatively robust. The available information was assigned a weight proportional to the number of data points for the posterior uncertainty estimates. The error defining the SOC priors was assumed to be equal to half the value, while the error for the climate factor was assumed to be equal to its variance by climate class. The initial apportioning of SOC between the old and young pools, akin to the initialization of any compartmental model usually performed with spin-up, was assumed to be an uncertain parameter, with a normal distribution centered around 0.93 and with a variance of 0.05 based on a former model calibration on the Ultuna long-term experiment (cut between 0.8 and 0.98). Aboveground estimation parameters priors were modeled according to normal distributions with error either assumed or estimated from the data or literature sources, except for the tree:crop ratio, which was modeled as uniform distribution. Although there is no presence of organic amendments in the Feliciano et al. (2018) database, this parameter was introduced in order to maintain model flexibility for future applications. It was set to 0.0 Mg ha −1 y −1 with an error of ± 0.1 t ha −1 y −1 based on the assumption that amendments utilized in agroforestry systems are likely rather limited. Inputs in agrosilvopastoral systems are here not considered. This is not particularly problematic if  L. Menichetti, et al. Agriculture, Ecosystems and Environment 303 (2020) 107118 the C grazed and then returned to the soil is photosynthesized in situ and so also implicitly distributed into the other pools. The annual increment in aboveground biomass was taken from Cardinael et al. (2018), while the maximum value was taken as the maximum value in the dataset for each class with an error corresponding to the variance of the aboveground biomass in each class. The ratio of tree:crop biomass was assumed to be: 0.35 for agrosilvicultural systems, 0.5 for home gardens and intercropping systems, 0.6 for cocoa, coffee, parklands, and silvopastoral systems, 0.9 for fast-growing tree systems, and 1.0 for woodlots. These are our expert-opinion assumptions based on what we considered most reasonable for each system and we associated those with a high uncertainty having an error term of ± 0.3 (and cut between 0 and 1). We considered woody root mortality roughly proportional to the life cycle of trees (Hammar et al., 2014) and assumed 20 ± 10 years for woodlots and 35 ± 5 years for all other classes, irrespective of climate. The prior distributions for root:shoot ratios of trees and crops were based on Cardinael et al. (2018), by calculating the ratio between aboveground and belowground increments, while the root:rhizodeposition ratio was set to 0.5 (according to Pausch and Kuzyakov, 2018), with the error assumed to be ± 0.1. The model was calibrated according to Bayesian principles with a Markov chain Monte Carlo (MCMC) sampler (specifically a standard Metropolis-Hastings sampler) within the programming framework JAGS (Plummer, 2003). The calibration was performed on one single chain based on 10000 runs. We used the posterior probability distributions of the kinetic parameters (k y , k y,w, k o , h s , h r , h w , h a ) after calibration to develop a predictive tool based on Bayesian principles. The probability distributions from the calibration are a discrete nonparametric approximation utilizing directly the population of parameter values extracted from the whole Markov chain Monte Carlo.

Climate normalization
The climate normalization in ICBM Agroforestry is deterministic and outside the Bayesian framework, and is supplied as model input. In this study, our aim was to develop a rapid but relatively low-detail method for predictions, and for this application we were using the map in Fig. 3 (also available as georeferenced raster as electronic supplement). In order to minimize eventual error, we suggest averaging the values on a regional scale, smoothing out eventual local variation. Scripts and documentation for more localized calculation of the climatic scaling are available on request from the corresponding author.

The sensitivity of model predictions to variations in model kinetic and input parameters
After calibration, we used the matrix defined by all the parameter sets of the MCMC (where one row consists of a parameter set randomly chosen based on the parameter prior probability distributions) to determine how variations in parameter values influence predictions. Here we consider both kinetic parameters and input parameters. All the parameters related to the decay kinetics are general for all the classes, but the parameters related to the C inputs to soil are specific for each class so they have more degrees of freedom. Therefore, we compared these two groups of parameters separately and normalized between 0 and 1 among each group to assess a ranking in sensitivity. We calculated for each parameter set of this matrix the resulting average root mean square error (RMSE) for all sites, and used these data in a sensitivity analysis according to the method by Hornberg, Spear and Young (HSY) proposed by Beven (2008). The method is based on distinguishing the whole posterior MCMC in two bins of parameter sets, "behavioral" and "no-behavioral", based on an arbitrary criterium, and then measuring the distance in the probability distribution of each parameter value between these two bins. The more the value of a certain parameter or variable is distant between the two bins, the more the model is sensitive to that parameter or variable. The distance was expressed according to Beven (2008), as Kolmogorov-Smirnov distance. The sensitivity of the parameters related to the C inputs and all the associated errors associated were considered by analyzing the model sensitivity to the various input classes variation. Agrosilvicultural systems were particularly sensitive to C inputs from roots, woodlots to variation in aboveground biomass C inputs, while silvopastoral showed high sensitivity to C inputs from both aboveground biomass and roots. Homegarden and shadow systems did not show particular sensitivity to any parameter, reflecting the high variance of these systems (Appendix 1).
Concerning kinetic parameters, the model was particularly sensitive to the humification rate of wood, followed by the decay rate of the old pool and the decay rate of woody material. The model predictions were less sensitive to the humification rate of shoots and roots. The humification rate of amendment and decay constant of the young pools were almost not influencing the model predictions (Appendix 2).

The model posterior parameters
The posterior model kinetic parameters (Table 1) did not differ much from the priors (Appendix 3 and Andren and Kätterer, 1997) because of the choice of strong (or informative) priors and a relatively conservative range of exploration. For these parameters, we decided to use strong priors because we wanted to maintain compatibility with previous versions of the model. Given the many degrees of freedom in the calibration (in particular the C input estimation parameters), we believe that model calibration already has enough flexibility to represent the results. Moreover, the determination of the old pool kinetics in the original ICBM version is tied to real data from a long-term bare fallow and is considered robust enough to be generalizable (Andrén et al., 2008(Andrén et al., , 2012Menichetti et al., 2019). The posterior distributions of the kinetics of the young pools were similar to the prior distributions. The reason for this small difference in the information contained in the calibration is related to the interactions between these terms and the C inputs, together with the small impact of the young pools to the total SOC stocks. The information of the posteriors for the kinetic parameter values is therefore coming mainly from the literature, and although it comes mostly from agricultural systems, the decay processes of organic matter are for the most part universal . We must also remember here that all compartmental SOC models are approximations of reality, and in every case, the different "pools" are arbitrary discretizations of a continuum. Parameter values depend on the assumptions coming with this discretization and are therefore not representing with perfect correspondence a specific physical reality. However, they are useful to produce meaningful predictions.
The uncertainty of the C input estimation parameters had to be set very high, with resulting non-informative priors, and this is reflected in the quite uncertain posteriors. Uncertainty of both priors and posteriors was particularly high for tree:crop ratio (prior set as uniform between 0.3 and 0.7). This uncertainty is representing the high variance in each class since the classification adopted here was rather broad, and the Table 1 The posterior values of the reparameterized model. The two k n parameters are expressed in unit of C per year so in this case Mg ha −1 y −1 , while the h n parameters are dimensionless.  Menichetti, et al. Agriculture, Ecosystems and Environment 303 (2020) 107118 systems can still differ considerably from one class to other in terms of C inputs since they might present very different species associations. Reducing this uncertainty would require a more granular dataset or more precise information in order to constrain the prior probability distributions.

Global maximum SOC stocks of agroforestry systems
Our model, given all our assumptions, estimated total SOC stocks at equilibrium of agroforestry systems for all classes as roughly comparable to reported global mean SOC stocks at equilibrium for agricultural systems, which are around 82 Mg C ha −1 (Zomer et al., 2017). However, there was great variation in equilibrium SOC stocks of agroforestry both between classes and within classes. The mean SOC stocks at equilibrium were between 156 and 263 Mg C ha -1 and the mode value (i.e., the value appearing most often in the chain, Dutta and Goswami, 2010) between 44.2 and 124 Mg C ha −1 , depending on the system (Table 2). In general, with the data at our disposal, the equilibrium SOC stocks of agroforestry systems seems comparable with that of agricultural systems. In many agroforestry system the mode of the equilibrium SOC stocks were higher than that of agricultural systems, in particular for systems in tropical environments (shadow systems and woodlots). Considering the means all systems seemed instead to store more C than agricultural systems, with particularly large values for homegardens and agrosilvopastoral systems.
All classes presented a very broad distribution of posterior steady state SOC stocks, with a very long tail to the right ( Fig. 5 and 6, Appendix 4). Since the variation in equilibrium values depends only on the C inputs to each system and not on the initial SOC stocks, it is independent of the land-use history of each site. Our results suggest, therefore, that the choice of systems play a crucial role in determining the net SOC sequestration of agroforestry systems while site history does not matter much.
From the data at our disposal, the most considerable potential for SOC sequestration lies within homegarden and silvopastoral systems, which showed very high mean values but lower modes, indicating a particularly right-skewed distribution. The impact of animals on decay processes in silvopastoral systems was not considered here due to the scarcity of available data to use as a proxy variable of animal density. Animals might modify the cycle of the young shoot pool through grazing, but the part of this pool that is then humified in the soil should not change much since animals would metabolize only easily available C otherwise consumed by soil organisms. The presence or absence of animals would probably therefore, have a negligible effect on the steady state C in the soil, which is mainly determined by net primary production, humification coefficients, and decay of the old pool. The probability distribution of woodlot systems in tropical environments also displayed a rather thick tail to the right, indicating the many exceptions gaining much SOC.
In general, the very long right-hand tails of the probability distributions for all systems (Appendix 4) suggest that agroforestry might in potential have even higher equilibrium SOC stocks than what suggested by the average estimates. However, it is not possible to develop specific management recommendations on this database alone. Advances in obtaining more SOC data over long time scales in agroforestry systems are key for the optimization of these systems. Gathering more data and improving therefore our knowledge on agroforestry as a C sequestration tool could be essential in the coming decades to combat climate change and to increase the resilience of food production systems. Against this background, we consider our calibrated first-order SOC model as an introductory tool for analyzing the equilibrium SOC stocks of agroforestry systems.

Model validation
We tested the ICBM Agroforestry predictive model on published data from some agroforestry studies, allowing us to make a validaton for nine case-specific simulations (Abaker et al., 2016;Beer et al., 1990;Fernández-Núñez et al., 2010 andNorgrove andHauser, 2013). In general, the results (Fig. 7) confirm the validity of the model predictions when considering its uncertainty intervals. ICBM Agroforestry predictions had an overall average error of 12.3 %. The model systematically under-predicted the observed trends of SOC stocks in the four cases from Spain (Galicia) in a temperate oceanic environment (Fernández-Núñez et al., 2010). This is not surprising given that the data on aboveground productivity used to develop our model did not include similar climates (the closest of the three cases from temperate environments in it comes from Russia, Dixon et al., 1994). The second site in Abaker et al. (2016) and the site in Norgrove and Hauser (2013) each presented some local specificities not considered in our model. For example, at the sites described in Abaker et al. (2016), the authors reported different degrees of disturbances (cutting, wind throw, grazing and natural mortality). Concerning the site described by Norgrove and Hauser (2013), the authors reported particularly low productivity for Cacao in southern Cameroon, which might explain the overestimation in SOC stocks with our model. These results clearly show the obvious difficulty developing a general model able to account for all the possible site-specific factors and climatic conditions. However, ICBM Agroforestry provides reasonable generic predictions once we include in the statistical treatment also the associated model uncertainty.
A Bayesian tool for reporting and comparing the potential equilibrium SOC stocks of different agroforestry systems Being a first-order compartmental model, the mathematical definition of our model is very similar to the Century model utilized in IPCC Tier 2 and Tier 3 estimates (Ogle et al., 2019), although simpler (two pools instead of five, but due to the mathematical similarities the two models will produce a very similar dynamic representation). The Bayesian approach for model calibration applied in this study allows for a flexible estimation of the model error, which can include additional knowledge available at each site. This approach, applicable to any other compartmental first-order model, represents the central value of this study.

Actual model limitations and future possibilities
The term "agroforestry" is rather loose and encompasses a vast diversity of agroecosystems. A major challenge for any agroforestry SOC model is the estimation of C inputs for the diverse plant associations occurring in these systems. To improve this source of uncertainty there is a need for much more ecologically specific predictors than the ones associated with the data at our disposal (e.g., some functional classification of plants and specific root:shoot ratios for different system). The model proposed in this study deals with variability and consequent uncertainty by considering it as a statistical error. Such uncertainty is a fundamental part of model predictions and must always be included.
The model must also deal with a relatively short-term perspective of available data, making long-term predictions challenging and increasing uncertainty in the distant future. As more data from mid-and L. Menichetti, et al. Agriculture, Ecosystems and Environment 303 (2020) 107118 long-term agroforestry experiments around the world become available, the model can easily be recalibrated by updating the posterior probability distributions according to the new data. A larger global dataset that includes spatial coordinates for experimental sites could also be developed, together with global soil maps, to introduce in the model an edaphic function (e.g., decay kinetics a function of clay content; see Bosatta and Ågren, 1997). Edaphic properties were not considered in the present study since we could not find enough studies reporting consistently the soil properties needed to calibrate such a response function. A more refined possibility for the end user would be to combine the information from this study with the new information present at the local site the user wants to model into a new Bayesian calibration run over a number of years selected by the user. This approach would produce predictions making also use of the information available to the user, improving precision, to estimate the future SOC stocks and, most importantly, the SOC stock equilibrium for a specific agroforestry system. This procedure is similar to run a new deterministic calibration, with the difference that Bayesian statistics preserves and updates the available model information combining it with new information. The distinction between calibration and prediction is then no longer exact, and the whole approach lies conceptually between calibration and prediction. This latter approach requires at this point some programming skills, but the needed information about the probability distributions from the calibration in this study are included as supplementary electronic material. Another possible use of the model would be to reprogram it in a tool, accessible to users, and where they could enter new experimental data, and where the present calibration could continuously be updated. For instance, having a Bayesian calibration run on a centralized server through a web interface, the server would store the new information and updating the model, producing better predictions for each end-user and improving the model accuracy with more data.
As one can see from our demonstrative application (accessible at https://ilmenichetti.shinyapps.io/ICBMAgroforestry/), it is also possible to make the model more accessible for users without specific programming skills, but this requires, first of all, to define specific application cases for it. Each specific model application will present a different equilibrium in the trade-off between generalization power and accuracy, and requires to develop a specific program specific for a certain class or level of detail. Users might not always be able to determine all the parameters required to run a certain model version, or might instead have more data at disposal than what the current model can utilize wasting some precision. Different model applications might also require a revision of model assumptions. Future steps required for the introduction of our approach to a broader community will require first of all stakeholders and requirement analysis to define the use cases and initiate a new project aimed at developing the tool.

Conclusions
Our results shows that some classes of agroforestry systems have a L. Menichetti, et al. Agriculture, Ecosystems and Environment 303 (2020) 107118 potential for SOC sequestration superior to agricultural landscapes and could represent an essential land management strategy for climate change mitigation. However, as noticed in previous studies, such potential is highly uncertain due to the large variability between different types of agroforestry systems. The model developed and calibrated in this study, ICBM Agroforestry , accounts for this variability by using a Bayesian framework. With the data at our disposal, the model predicts that the most promising agroforestry systems in terms of SOC sequestration seem to be woodlot, homegarden and shadow systems in semiarid and tropical lands, but that potentially all agroforestry systems could sequester a similar amount of SOC compared to agricultural systems.
Our sensitivity analysis revealed that the C inputs were by far the most crucial factor determining model fit to the data, so we recommend that experimental efforts focus on accurate estimation of inputs and their variation over time (possibly over decades).
The model can be used to predict SOC stocks at equilibrium, which are independent of time and previous land-use history of sites, allowing our approach to make robust comparison of C balance between ecosystems and land-use changes. Whenever a specific time perspective is set, the model can also be used for robust C balance comparisons over a certain time (present and future). At this point, we propose the ICBM Agroforestry model mainly for large-scale applications. The statistical error is a fundamental part of the model predictions within this framework. The statistical approach utilized in this study could also easily be applied to similar models.

Declaration of Competing Interest
The authors report no declarations of interest.