Carbon sequestration in agroforestry systems can be assessed by accounting for the different C stocks and flows. When trees are planted on a field, C is allocated in their above- (trunk and branches) and belowground biomass (roots) during their growth. The rate of tree growth and biomass accumulation depends strongly on the species growth strategy, and may thus differ strongly among tree species. Meanwhile, trees shed their leaves every year, and leaf C is gradually transferred to the soil and converted into SOC upon the leaves’ decomposition. The pattern of leaf fall around the trees is also highly species dependant as well as the decomposability of the leaves and thus the rate at which the leaf C is converted into SOC. On top of that, C contained within crop residues that are left on the field are also converted into SOC upon decomposition. CARAT is designed to simulate these C turnover rates as well as the spatiotemporal variation in these above- and belowground C stocks for a customized agroforestry field. CARAT does not yet account for C inputs to the soil from root decomposition or root exudates.

## 2.1. Field design and species selection

In practice, the field is defined by the user with given dimensions, and subsequently divided into cells of 1 m². For each of these cells, the amount of leaf litter and associated C inputs to the soil are modelled. These calculations are repeated through time with an interval of one year. The evolution of soil C stocks can thus be estimated for a given time period (currently maximum 30 years) at any point in the field. Secondly, the C stocks in the woody biomass of trees on the field can be simulated for the same time period. The tool is now developed for nineteen different tree species that typically occur in temperate agroforestry systems (Table 1). Alternatively, a standard tree species can be selected for which the properties are based on the average of all tree species. The user can determine the exact position of individual trees in the field.

Table 1

Overview of the nineteen different tree species included in CARAT with their respective growth rate, diameter at breast height (DBH)-age relationship and biomass equation as well as the source from which the different equations were compiled. The DBH of the trees is expressed in cm, the age (t) in years and the biomass (B) in kg.

Species | Growth rate | Allometric equation | Source |

*Acer pseudoplatanus* | Fast | DBH = 61/(1 + exp(2-0.059*t)) | Monteiro et al. (2017) |

| | B = 520*(0.0019421*DBH^1.785) | Shafer et al. (2019) |

*Alnus glutinosa* | Fast | DBH = 728/(1 + exp(4.62 − 0.035*t)) | Van De Berge et al. (2021) |

| | B = 0.00079*DBH^2.28546 | McPherson et al. (2016) |

*Aesculus hippocastanum* | Medium | DBH = 52.8/(1 + exp(2.02–0.069*t)) | Lukaszkiewicz & Kosmala (2008) |

| | B = 580*(0.0002431*DBH^2.415) | Shafer et al. (2019) |

*Corylus avellana* | Medium | DBH = 89.9/(1 + exp(2.28 − 0.04*t)) | Lukaszkiewicz & Kosmala (2008) |

| | B = exp(-2.2118 + 2.4133*ln(DBH)) | McPherson et al. (2016) |

*Fraxinus excelsior* | Fast | DBH = 68.7/(1 + exp(2.7-0.072*t)) | Lukaszkiewicz & Kosmala (2008) |

| | B = 530*(0.0005885*DBH^2.206) | Shafer et al. (2019) |

*Juglans regia* | Medium | DBH = 89.9/(1 + exp(2.28 − 0.04*t)) | Lukaszkiewicz & Kosmala (2008) |

| | B = exp(-2.2118 + 2.4133*ln(DBH)) | McPherson et al. (2016) |

*Malus domestica* | Medium | DBH = 89.9/(1 + exp(2.28 − 0.04*t)) | Lukaszkiewicz & Kosmala (2008) |

| | B = exp(-2.2118 + 2.4133*ln(DBH)) | McPherson et al. (2016) |

*Populus* × *canadensis* | Very fast | DBH = 6.47185 + 16.32306*log(t) | Carton (2022) |

| | B = 374*(exp(-2.14 + 2.26*ln(DBH))) | Hjelm & Johansson (2012) |

*Prunus avium* | Fast | DBH = 89.9/(1 + exp(2.28 − 0.04*t)) | Lukaszkiewicz & Kosmala (2008) |

| | B = exp(-2.2118 + 2.4133*ln(DBH)) | McPherson et al. (2016) |

*Pyrus communis* | Medium | DBH = 89.9/(1 + exp(2.28 − 0.04*t)) | Lukaszkiewicz & Kosmala (2008) |

| | B = exp(-2.2118 + 2.4133*ln(DBH))) | McPherson et al. (2016) |

*Quercus petraea* | Slow | DBH = 60.2/(1 + exp(2-0.061*t)) | Van De Berge et al. (2021) |

| | B = 580*(0.0002431*DBH^2.415) | Shafer et al. (2019) |

*Quercus robur* | Slow | DBH = 60.2/(1 + exp(2-0.061*t)) | Van De Berge et al. (2021) |

| | B = 580*(0.0002431*DBH^2.415) | Shafer et al. (2019) |

*Robinia pseudoacacia* | Fast | DBH = 1452.5/(1 + exp(5.66 − 0.041*t)) | Van De Berge et al. (2021) |

| | B = 4.13741*DBH^2.17752 | McPherson et al. (2016) |

*Salix* sp. | Fast | DBH = 1452.5/(1 + exp(5.66 − 0.041*t)) | Van De Berge et al. (2021) |

| | B = 0.454*10^(0.8856 + 2.0552*log(2.54*DBH)) | McPherson et al. (2016) |

*Sorbus aucuparia* | Medium | DBH = 89.9/(1 + exp(2.28 − 0.04*t)) | Lukaszkiewicz & Kosmala (2008) |

| | B = exp(-2.2118 + 2.4133*ln(DBH)) | McPherson et al. (2016) |

*Sorbus torminalis* | Medium | DBH = 89.9/(1 + exp(2.28 − 0.04*t)) | Lukaszkiewicz & Kosmala (2008) |

| | B = exp(-2.2118 + 2.4133*ln(DBH)) | McPherson et al. (2016) |

*Tilia cordata* | Medium | DBH = 89.9/(1 + exp(2.28 − 0.04*t)) | Lukaszkiewicz & Kosmala (2008) |

| | B = exp( -2.6788 + 2.4542*ln(DBH)) | McPherson et al. (2016) |

*Tilia platyphyllos* | Medium | DBH = 89.9/(1 + exp(2.28 − 0.04*t)) | Lukaszkiewicz & Kosmala (2008) |

| | B = exp(-2.6788 + 2.4542*ln(DBH)) | McPherson et al. (2016) |

*Ulmus* sp. | Fast | DBH = 61/(1 + exp(2-0.059*t)) | Monteiro et al. (2017) |

| | B = 460*(0.0018*DBH^1.869) | Shafer et al. (2019) |

## 2.2. C stocks in woody biomass

To estimate the aboveground woody biomass and associated C stocks for the different tree species on the field, state-of-the-art allometric relations were used based on the diameter at breast height (DBH; in cm). First, the DBH of each tree was modelled for every time step using a logistic regression (Eq. 1).

\(DBH=\frac{a}{{e}^{b-c t}}\) [Eq. 1]

Where \(a\), \(b\) and \(c\) are species-specific model parameters, and \(t\) the age of the tree (in years). The model parameters were estimated per tree species based on literature data: Van Den Berge et al. (2021) for trees growing in hedgerows in Flanders, Lukaszkiewicz & Kosmala (2008) for urban trees in Poland, Monteiro et al. (2017) for urban trees in Great Brittan, and Carton (2022) and Pardon (2018) for line plantings of poplars in Flanders. From these studies, the raw data was used and the parameters in Eq. 1 were estimated by minimizing the Root Mean Square Error (RMSE). For some species, no data was found in literature, and the DBH-age relationship of a species with similar growth rate was chosen (see Table 1). For some other species, two or more relationships could be constructed from literature but then the one representing its growth in a solitary context (e.g. urban trees) was chosen (**Figure S1**).

The estimated DBH values are then used as input for species-specific allometric equations that model the trees’ aboveground biomass for each time step. These equations were again extracted from literature: Shafer et al. (2019) with biomass equations for trees in food forests, McPherson et al. (2016) with volume equations for urban trees in the US, and Hjelm & Johansson (2012) with volume equations for poplars in Sweden. For some species, no allometric equations could be found in literature, and hence the equation of a species with similar growth rate was used (see Table 1). If only an allometric equation for aboveground volume was available, the biomass was calculated by multiplying the volume with the wood density of the specific species (values extracted from the Wood Density Database; Zanne et al. 2009).

To obtain the total woody biomass of each individual tree, a root-to-shoot ratio of 0.26 was assumed (average value for temperate deciduous forests; IPCC 2006). It should be noted, however, that root systems vary with species, age, climatic zone and other regional environmental conditions, and may also differ for trees growing in a solitary or linear context (Ye et al. 2019). Given the current lack of data, this variation could not yet be incorporated. Finally, the woody biomass of all individual trees on the field was summed and multiplied with a conversion factor of 47% (typically used for trees in the temperate zone; IPCC 2006) to estimate the total C stocks (expressed in kg C ha− 1).

## 2.3. C stocks in leaf litter

To estimate the C input via leaf litter to the soil, a spatial leaf distribution model was used. In this tool, we opted for a simple leaf model by Ferrari and Sugita (1996) that does not account for wind direction (Eq. 2). This model estimates the amount of leaf litter (\({LF}_{ij}\); g dry matter m− 2 year− 1) for each tree species j at a location i in the field.

\({LF}_{ij}=\frac{{\alpha }_{j}{\gamma }_{j}^{2}}{2\pi }\sum _{k}{DBH}_{jk}^{{\beta }_{j}}\text{e}\text{x}\text{p}(-{\gamma }_{j}{d}_{ijk})\) [Eq. 2]

Where \({\alpha }_{j}\), \({\beta }_{j}\) and \({\gamma }_{j}\)are specific model parameters for tree species j, \({DBH}_{jk}\) the DBH of the kth individual of tree species j (in cm), and \({d}_{ijk}\) the distance of location i to the kth individual of tree species j (in m). This model assumes that the amount of leaf litter reduces exponentially with the distance to the tree. The three specific model parameters (α, β and γ) were determined for each tree species in this tool based on Ishihara & Hiura (2011). This study reports the estimated model parameters for 17 different tree species in the temperate deciduous forests of Japan. By comparing the shape and size of the tree species included in this tool with the tree species from the Japanese database, we completed the corresponding parameters (**Table S1**). The DBH of each tree increases over time, and was estimated using the logistic model described above. This way, the amount of leaf litter produced by every individual tree could be estimated at any point in the field and for any point in time. The total amount of leaf litter per location (grid cell) and time step was then calculated by summing the leaf litter of all individual trees on the field.

Subsequently, the total organic C content in the leaf litter was again calculated using the conversion factor of 47% (IPCC 2006). Once the leaves have fallen, the amount of C contained within them will over time be partially converted to the soil. To model this process, the organic C content of the leaves is divided into three fractions: labile (\({f}_{DPM}\)), recalcitrant (\({f}_{RPM}\)) and humified (\({f}_{HUM}\)) organic matter. These fractions are determined using the approach described by Peltre et al. (2012). First, the \({I}_{ROC}\) indicator value is calculated, which mirrors the amount of organic matter that is contained in the soil for multiple decades (Eq. 3).

\({I}_{ROC}=44.5+0.5 SOL-0.2 CEL+0.7 LIC-2.3 {C}_{3d}\) [Eq. 3]

Where \(SOL\), \(CEL\) and \(LIC\) are the soluble, cellulose and lignin fractions [%] in the leaf material, respectively (Van Soest 1991). \({C}_{3d}\) represents the fraction of organic matter that has been mineralised after three days in the lab. This variable is strongly correlated with the lignin fraction in the leaf material. Next, the data from Peltre et al. (2012) were used to model a linear regression between the \({C}_{3d}\) and \(LIC\) values (Eq. 4). This equation could then be adopted to calculate the \({C}_{3d}\) fractions of other types of leaf material given that its Van Soest fractions are known.

\({C}_{3d}= -0.2716 LIC+10.769\) [Eq. 4]

Using the regression from Peltre et al. (2012), the fractions of labile (\({f}_{DPM}\)) and recalcitrant (\({f}_{RPM}\)) organic matter in the leaves could also be calculated (Eq. 5 and Eq. 6, respectively).

\({f}_{DPM}=-1.254 {I}_{ROC}+115.922\) [Eq. 5]

\({f}_{RPM}=0.979 {I}_{ROC}-8.928\) [Eq. 6]

Eq. 3, 4, 5 and 6 were then used to model the organic matter fractions in the leaves of different tree species included in De Schrijver et al. (2012) (*Alnus glutinosa*, *Fagus sylvatica*, *Tilia cordata*, *Quercus rubra*, *Prunus avium*, *Populus* × *euramericana* and *Fraxinus americana*), and calculate their respective \(DPM:RPM\) ratio (Eq. 7 and **Table S2**).

\(DPM:RPM= \frac{{f}_{DPM}}{{f}_{RPM}}\) [Eq. 7]

Using these values, a distinction was made between trees with slow (DPM:RPM = 0.62) and fast (DPM:RPM = 0.80) decomposing leaf material. The nineteen tree species included in CARAT were then assigned to one of these two categories based on the decomposition rate of their leaf material (from literature and expert knowledge) (**Table S3**).

## 2.4. C stocks in soil

The evolution of SOC was simulated using the RothC-model, which uses a monthly time step to calculate turnover of organic carbon in non-waterlogged topsoils (Coleman and Jenkinson 2014). The model considers five organic C pools: decomposable plant material (DPM), resistant plant material (RPM), microbial biomass (BIO), humified organic matter (HUM), and inert organic matter (IOM). Organic inputs, both from plants and manure, are split between the DPM, RPM, and HUM pools based on the origin of the material. Each time step, the model calculates the turnover between the different pools, or the mineralization into CO2, using specific decay rate constants. This turnover is also a function of the soil clay content, and of rate modifying factors that account for soil water content, air temperature, and soil cover. For each active C pool containing Y ton C ha− 1, the amount of material that decomposes each timestep is calculated according to Eq. 8.

\({Y}_{Decomposed}=Y(1-{e}^{-abckt})\) [Eq. 8]

Where a, b, and c are the rate modifying factors for temperature, moisture, and soil cover, k is the decomposition rate constant for each pool, and t is 1/12 to account for the monthly timestep.

Model inputs are both the amounts and the DPM:RPM-ratio of the incoming organic material. For CARAT, three land-use scenarios are pre-defined which are common in temperate Europe: (1) a rotation of potato, winter wheat, sugar beet, and maize, with a manure dose of 25 ton ha− 1 year− 1 pig slurry, (2) monoculture maize with grass as a cover crop and a manure dose of 35 ton ha− 1 year− 1 cattle slurry, and (3) perennial grassland with a manure dose of 25 ton ha− 1 year− 1 cattle slurry (**Table S4**). To account for the effect of leaf fall in an agroforestry system, the monthly leaf fall calculated with the spatial leaf distribution model was added as an additional input. To differentiate between agroforestry and conventional agriculture, simulations are carried out both with and without the organic input from tree leaf fall.

Further input for RothC includes data concerning the soil characteristics (initial SOC content, soil clay content, and depth of the modelled layer), and monthly climactic data (air temperature, cumulated precipitation, and evapotranspiration). Evapotranspiration is calculated from daily climatic data using the Penman-Monteith equation (Allen et al., 1998).

Initial SOC content was divided between the DPM-, RPM-, HUM-, and BIO-fractions based on a 100-year warming-up simulation of the RothC-model for common crop rotations (potato-sugar beet-winter wheat and monoculture silage maize) and soil types (sand, sandy loam and silt) in Flanders, Belgium. The crop rotations considered were potato-sugar beet-winter wheat and monoculture silage maize (**Table S5**). The final average fractions of DPM, RPM, HUM and BIO for Flanders, Belgium were respectively 1.4%, 13.3%, 83.2%, and 2.1%.

## 2.5. Online interface

The CARAT model was implemented in R (R Core Team 2023), and an online interface was developed using the Shiny-package (Chang et al. 2024). The interface is divided into four parts, which are separated into tabs (Fig. 1). First the ‘Field’-tab allows the user to create a rectangular field with custom dimensions, or to select an existing agricultural parcel on a map of Flanders. Subsequently, different tree species can be chosen from a database of nineteen different species and placed on the field (illustration for a fictitious agroforestry system on an existing field in Flanders in **Figure S2**). Trees can be positioned individually or according to a grid with predefined inter- and intrarow distances. In case of a mistake, trees can also be removed from the field, either individually or per species. In the ‘Time series’-tab, the soil parameters needed for the RothC-model are specified by the user (soil clay content, bulk density, depth of the modelled layer, initial SOC, and crop rotation) and the simulation can be carried out for a period of 30 years. Upon completion, the model output is visualized in this tab as a time series which shows the build-up of SOC under agroforestry and conventional management as a function of distance to the trees. This output is further visualized as a raster plot for a chosen timestep in the next tab, ‘Raster output’. This tab allows the user to choose a year within the 30-year simulation period and display the expected leaffall (kg C ha− 1 year− 1), expected SOC (%) and expected SOC increase due to leaffall (kg C ha− 1 year− 1) for this timestep (**Figure S3** illustrates the output after 30 years for the fictitious field agroforestry system specified above). Finally, in the last tab ‘Summarizing table’ an overview table is generated showing the increase in SOC due to agroforestry, and the total C in the trees (both in kg C ha− 1 and in kg C field− 1).

The CARAT application V1.0.0 is freely available online at https://shiny.ilvo.be/PLANT/CARAT_v1.0.0, and a detailed practical manual including examples of simulations can be accessed via https://ilvo.vlaanderen.be/uploads/documents/CARAT-manual.pdf.