Implementation and initial calibration of carbon-13 soil organic matter decomposition in Yasso model

. Soil carbon sequestration has gained traction as a mean to mitigate rising atmospheric carbon dioxide concentrations. Veriﬁcation of different methods’ efﬁciency to increase soil carbon sink requires, in addition to good quality measurements, reliable models capable of simulating the effect of the sequestration practices

It is essential to understand which soil processes are the most important for soil C sequestration. The soil C pool can be divided into different fractions based on their chemical composition, physical characteristics or assumed turnover or residence times (Poeplau et al., 2018). Soil processes in general are complex as biological, chemical and physical drivers act simultaneously. For modelling purposes, the fate of carbon-13 isotope ( 13 C) gives valuable additional information of the parameter values of used model as the 13 C signatures are sensitive indicators of changes in processes. Soil organic matter (SOM) consists 25 of molecules with different carbon isotopes. In theory, molecules with lighter 12 C atoms have lower activation (kinetic) energy requirements than those with 13 C. This leads to easier decomposition of 12 C-bearing compounds and enrichment of 13 C in residual organic molecules (Fry, 2006). By estimating 13 C in different fractions of SOM or varying residence times and adding 12 C/ 13 C reaction kinetics into the models would allow verification of the model functioning, and improve model predictions.
13 CO 2 measurements associated with gas flux measurements provide a promising way to link soil models to ecosystem mod-30 els and allow further integration to earth system models where 13 C isotopes are used to detect large scale C cycling patterns (Flanagan et al., 2005).
There are a multitude of ways to improve MRV (Smith et al., 2020), but in our experience one method has not been as readily examined -that of carbon isotope composition in the soil and in heterotrophic respiration. The reason behind the lack of such examinations is simple, such approaches require a model that can reliably represent the soil organic carbon (SOC) dynamics 35 for different carbon isotopes while still retaining relatively straightforward structure. The latter is especially important when we take into account the lack of good-quality calibration and validation data.
In this paper we introduce a simple 13 C isotopic circulation into the recently re-calibrated SOC model Yasso (Viskari et al., 2021(Viskari et al., , 2020Tuomi et al., 2011). In our approach, the decomposition of 13 C-specific soil organic matter ( 13 C-SOM) is linearly dependent on the default Yasso model parameters, the carbon isotope fraction 13 C/ 12 C and a new scaling factor θ, 40 that represents change to the decomposition rate between the carbon isotopes. The underlying hypothesis behind this design is that since 13 C has a larger atomic weight it is therefore not as reactive as 12 C, but environmental factors should still affect the decomposition of SOM, containing either isotope, similarly. We calibrate the new 13 C-related decomposition parameters (θ) and assess the model functionality both on short and long term (50-year simulation) basis.
Our aim is to improve Yasso20 model parameterisation (Viskari et al., 2021) to include 13 C/ 12 C reaction kinetics in the 45 model by using empirically measured SOM and 13 C data. We hypothesize that measuring 13 C in soil organic matter fractions 1) detects differences in the pool 13 C content supporting the 13 C-fractionation and enrichment theory, and 2) allows model development for significant improvements in SOM decomposition predictions.

50
The SOC measurements were derived from experiments described in Straková et al. (2012Straková et al. ( , 2011Straková et al. ( , 2010, where different types of plant litter was left to decompose inside litterbags in natural environment at Lakkasuo, a raised bog complex in Central Finland (61.8 • N, 24.3 • E, 150 m.a.s.l.). We utilised data detailing the conditions for pine branch and pine needle specific litterbags. In addition to determining the initial states for both litter types, 14 litterbags describe the soil conditions for pine branches and seven for pine needles at later stages of decomposition during the four-year-long experiment.

55
The litter was characterized by dividing it into carbon fractions by sequential extractions and hydrolysis according to Hilasvuori et al. 2013 (and references therein), also called AWEN extraction (acid, water, ethanol, non-soluble). In short, this included analysing the amounts of nonpolar extracts (corresponds to E), polar extracts (W), acid hydrolysable substances (A) and non-soluble Klason type substances (N). Air dried litter material was ground in a mill (Fritsch) to pass the 0.5 mm sieve and weighted into a centrifuge tube (35 ml). The amount of extractables was determined through the remaining mass after 60 shaking (2h or 18h; 250 rpm) with the different solvents followed by filtering through glass crucibles (Robu, Borosilicat 3.3 por. 4). At the start of the extraction procedure 0.5 g litter mass was used. Dichloromethane (CH 2 Cl 2 ; 15 ml; repeated twice) was first used to remove the nonpolar extractives. 0.35 g of the remaining dried (105 • C) solid sample was weighted again into a centrifuge tube and hot water (80 • C; 15 ml) was added and kept in a water bath (80 • C; 18 h). After centrifugation (1500 × g) the pellet was washed with 30 ml hot water to finish the extraction for polar extractives. In all cases the respective extractives were combined and dried. Evaporation was used for the nonpolar fraction and warming (50 • C) followed by freeze drying was polar fraction. 0.1 g oven dried (105 • C) material from the residue left after the hot water extraction was weighted into a centrifugation tube and 1.25 ml 72% sulphuric acid (H 2 SO 4 ) was added and shaken in room temperature (1 h; 250 rpm). Thereafter 35 ml water was added and incubated in a water bath (95°C; 18h) followed by filtration. The remaining mass (Klason lignin) was washed once with hot water (95 • C; 30-40 ml) and the mass was dried o/n in 105 • C. Each fraction ie. the 70 original litter, the solid remains after dichloromethane, water and acid extraction and from the evaporated nonpolar and polar extractants, subsamples were analysed for their relative 13 C/ 12 C ratios as δ 13 C values. The definition of δ 13 C is given below, where ( 13 C 12 C ) standard = 0.01123720 is the Vienna Pee Dee Belemnite (Craig, 1957, VPDB).
The isotopic composition of carbon was measured on a NC2500 elemental analyzer coupled to a Thermo Scientific Delta V The meteorological variables required to run the Yasso model were extracted from a nearby weather station measurements (Kolari et al., 2009), located at Hyytiälä (61.85 • N, 24.29 • E, 180 m.a.s.l.). We gathered monthly temperature and annual precipitation from the beginning of year 2005 to the end of 2008. Additionally we calculated averaged monthly temperature and averaged annual precipitation from years 2000-2014 to be used in simulating long-term carbon decomposition.

Yasso model
We generate the soil carbon pools utilising the Yasso SOC model (Liski et al., 2005). The underlying model version is the Yasso20 (Viskari et al., 2021) with parameter values given in appendix The state vector (x t ), representing the C content in AWENH pools at time t, is calculated by operating the state transition 100 matrix (M ) on the state vector of the previous time step (x t−1 ) and adding litter input (b t ), which in our simulations is set to zero (as we assume no litter is added into the litterbags after the beginning). The model initial state (in our simulations) is set to match the first measurements. The matrix M determines the decomposition of SOM and the flow of carbon between the different pools and it is dependent on various parameters as well as temperature and precipitation. We introduce 13 C-SOM decomposition into the Yasso model by adding separate 13 C-specific storages for each AWENH pool and including an 105 additional 13 C-specific SOM decomposition step. The input data (essentially carbon content) is first separated into 13 C-specific content and the rest, which we call 12 C although it also contains 14 C. The Yasso model is first run normally with 12 C as in Eq. 2, which is followed by 13 C decomposition using a modified version of the state transition matrix M . We modified the diagonal elements of M (these determine the SOM decomposition within each pool) by replacing the parameter (α) affecting the diagonal element with: Essentially, we include a dependency for the mass ratio of the carbon isotopes ( 13 C/ 12 C) as well as a free parameter θ for each AWEN pool separately. We didn't include a parameter for the humus pool (H) as we did not have measurements to calibrate the related parameter. We also note here that the diagonal elements are further dependent on temperature and precipitation, but these model aspects were not modified. Further details on Yasso model can be found in Viskari et al. (2021).

Model calibration
We calibrated the four θ parameters related to the decomposition of each AWEN pool 13 C-SOM. The objective function (f ) of the calibration is the cumulative squared error of the observed and modelled δ 13 C values: Here the summation is taken over all AWEN pools and available litterbag measurements (with measurements indicating zero 120 concentration for total carbon content removed from the calculations). The unnormalised (pointwise) parameter likelihood is Since we had only four parameters to calibrate, we set similar initial limits for the parameters and a suitable increment that determined how densely the parameter values were distributed. We ran the model with every member of the parameter "grid" to get an estimate of the overall shape of the parameter likelihood and to test for reasonable limiting values for the parameters.

125
This process was repeated several times with refocused grid and readjusted increment. as likelihood values started to plateau). The elongated shape of the likelihood dependencies between two parameters indicate some correlation between the parameters (Fig. 1). When we examine the parameter combinations with highest likelihoods (top 250 values), the strongest correlations (≈ 0.77) are present between θ A and θ W , θ A and θ N as well as θ W and θ N . The default and optimised parameter values were used to generate SOM decomposition and related C, 13 C and δ 13 C timeseries from the given initial states (Fig. 2). The differences between the simulated 13 C concentrations are too small to be 135 evident (C concentrations are identical), but we get a clear signal from the δ 13 C values. Overall, both model versions tend to underestimate the speed of SOM decomposition (the C and 13 C concentrations) at Lakkasuo for A and N pools and overestimate for W pool. The default Yasso model is reducing the relative 13 C content (reducing the δ 13 C values) for A and W pools and deviating from the observations whereas the optimised model version seems to be increasing the relative 13 C content and The Lakkasuo initial states and generated average year (averaged monthly temperature and annual precipitation) from years 2000-2014 were used to simulate a 50-year long carbon decomposition (Fig. 3). This simulation can be compared to Lakkasuo peat column δ 13 C values at different depths (Hilasvuori et al., 2013,  to the E pool and polar extracts to the W pool. A noteworthy detail is that on short term (Fig. 2) the default model increased the relative 13 C content (δ 13 C values) of E and N more than the optimised version, but on longer timescale this situation is reversed for E and projected to reverse for N (Fig. 3).  Table 2) at different depths.
We have introduced simple modifications to the Yasso model in order to account for 13 C-SOM decomposition. Incorporation of δ 13 C on SOM decomposition models is a necessary step towards integration of Earth system and dynamic land ecosystem models. The δ 13 C values of different organic compounds or chemical fractions of mixed organic material can be used as natural tracers which provide a unique tool to investigate and uncover complex decomposition processes in soil.
In the current study, we introduced new θ parameters to account for 13 C-SOM decomposition in the Yasso model. The calibration of these parameters only depend on the δ 13 C values, i.e. the relative carbon isotope concentrations, and revealed unimodal distributions for all four AWEN pool related parameters. Considering the acquired optima and taking into account that generally the ratio 13 C/ 12 C ≈ 0.01, then the new 13 C-SOM decomposition utilises values that differ at maximum 3‰ (for θ A ) from the default decomposition parameter values. Therefore, it is not surprising that both default and optimised model versions generate nearly identical SOM decomposition both on a short (Fig. 2) and long term (we did not provide this image 160 as it provides no additional value).
The acquired optima for θ A ,θ E and θ W are all negative, which in the model translates to reduced 13 C-SOM decomposition rate (Fig. 2). Likewise, the positive value for θ N implies increased 13 C-SOM decomposition when compared to the default model. However, the reduction in δ 13 C values, when compared to the default model version, is only true on shorter timescales ( Fig. 2) as each pool has a clear trend to increase relative 13 C content during the 50-year long simulation (Fig. 3). This is 165 due to the reduced 13 C-SOM decomposition in other pools -as there is more 13 C present in these pools, there is more 13 C available to be transported into the N pool, which compensates for the increased decomposition. This situation is not ideal, but understandable as Yasso also operates in the terms of "bulk" soil -there are no soil layers or differences in the effects of temperature, precipitation or Q10 at different soil depths (see e.g. Fig. 1 in Hilasvuori et al. (2013)).
The straightforward changes to the Yasso model have improved the model capabilities in reproducing observed δ 13 C values 170 in short (Fig. 2) and longer timescales (Fig. 3). Results from the 50-year long simulation seem to corroborate the initial hypothesis for A,W and N pools that the relative 13 C content in soil (larger δ 13 C values) increases with time. The optimised model even yields a positive trend for E pool δ 13 C whereas the default model tends to converge the δ 13 C values of all pools to roughly -30. The peat decomposition at different depths (Hilasvuori et al., 2013, Table 2) can be naively approximated to be of different ages at 10-year intervals. The optimised model behaviour follows the trend of these measurements and the results are 175 highly encouraging, even though the model is driven with a single averaged year representing the meteorological conditions from the beginning of the 21st century.
As estimation and modelling of soil organic matter decomposition, but also C sequestration, are current scientific chal- straightforward to produce, or meta-analysis using literature-based values could be also used for further evaluation across varying scales (local, regional, global).

Conclusions
We have demonstrated how to incorporate 13 C-SOM decomposition into the Yasso model and calibrate it. The model modifications were simple and straightforward and resulted in greatly improved simulated δ 13 C values. The capability of a model to simulate soil 13 C content and to better simulate SOM decomposition improves the applicability of Yasso-C13 model to scale process from ecosystem level to regional and global using δ 13 C as a tracer. Conceptually the presented work is on solid ground, 190 but the lack of suitable calibration and validation data urges further studies that produce new, precise experimental δ 13 C data suitable for Yasso-C13 model calibration and validation. Likewise further efforts should be made to include soil layers into the Yasso model.
Code and data availability. The data required to reproduce the results is available at Zenodo portal (Mäkelä, 2021a). The Yasso model source code and R interface are also available at Zenodo (Mäkelä, 2021b) or as "C13" branch at https://github.com/YASSOmodel/Ryassofortran. Appendix A: Yasso model parameters