Modeling persistence of coarse woody debris residuals in boreal forests as an ecological property

,


INTRODUCTION
The role of forests in meeting carbon targets Because of the central roles of forests in the global C cycle, a precise quantitative assessment of the carbon (C) balance of forests is particularly important for the global greenhouse gas balance.Key targets set by the European Commission for 2030 are at least 40% cuts in greenhouse gas emissions compared with 1990 levels (European Commission 2013).Forests can provide bioenergy, seen as an important component of future energy policies in Sweden (B örjesson Hagberg et al. 2016) and globally (Souza et al. 2017).In the EU, at least 32% of the energy should come from renewable sources by 2030 (European Commission 2013).However, the assumption that bioenergy is carbon-neutral has been questioned in recent years (McKechnie et al. 2011, Schulze et al. 2012, Gustavsson et al. 2015 Chen 2015, Liu et al. 2018).Wood has also many other uses, from short-term ones such as paper, cardboard, and low-cost furniture (which ultimately might end up in energy production for recovery, Medeiros et al. 2017) to longer-term C storage such as timber for construction or fibers for insulation (Pe ñaloza et al. 2016).The C balance of all these possible applications of plant materials, including being left in the forest, depends on the residence time of the C involved.
Life cycle assessment studies suggest that the fluxes resulting from changes in forest C stocks, influenced by forest management, are one of the key determinants of the climate impact of forest products (Kilpeläinen et al. 2011, Jäppinen et al. 2014).Therefore, a proper assessment of such impact must include the emissions from the utilization of the raw materials, but also the emissions that would be generated if raw materials would be left decomposing in forests (Repo et al. 2012).Depending on how the forest C balance is considered, the environmental impact of various forest components can change crucially (McKechnie et al. 2011, Stendahl et al. 2017).There is a large reported variation in decomposition rates of organic material in forests, which can affect the whole forest C balance and is not considered in the current decomposition models.

Ecosystem services from coarse woody debris
Among the forest components that can be extracted or left on-site, coarse woody debris (CWD) refers to dead trees and large branches and stumps on the ground.Here, we refer to CWD as >4.5 cm diameter, including dead standing trees (snags).If removed, this residual biomass is a potential energy source.A considerable part of biomass that would become CWD can be removed at the final harvest, as logging residuals, but it is usually left in the forest during thinnings.However, harvesting systems collecting even small-diameter stems for bioenergy have been designed (Sängstuvall 2018).
Recent forest management approaches have the two opposite objectives of maximizing the amount of CWD left in forests, and at the same time extracting as much raw material as possible (Kabrick et al. 2019).Generally, CWD harvest negatively affects forest biodiversity (Persson et al. 2013, Kline et al. 2015), especially wood-dependent organisms (such as beetles, fungi, and lichens), under-storey vegetation (Ranius et al. 2018), and birds (Sauerbrei et al. 2017).Pest saproxilic species are favored more compared to other saproxilic species by some forms of deadwood (Ranius et al. 2018).Snags in particular have a positive impact on birds, while the amount of CWD does not seem to affect mammals, reptiles, and amphibians (Riffell et al. 2011).
Coarse woody debris also plays a crucial role in nutrient cycles (Wiebe et al. 2014), and acts as a nutrient buffer maintaining long-term soil fertility (Herrmann and Bauhus 2018).In particular, CWD can act as a nitrogen (N) sink during the decomposition process (Palviainen et al. 2010), promoting initial N immobilization followed by a slow-release increasing the N efficiency of the stand (Wiebe et al. 2012), and increasing total N stocks (Wiebe et al. 2014).Additionally, CWD improves soil characteristics such as pH and physical properties (Moghimian et al. 2020).

The challenges in modeling the variability of CWD decomposition kinetics
Decomposition models can be used for quantitatively considering the trade-offs for CWD dynamics between different management options.The available models, e.g., Q or Yasso (Hyv önen et al. 1998, Liski et al. 2005), are based on similar kinetics, influenced by the time since tree death and initial diameter, but do not consider the position (standing or lying) of the wood ( Ågren andHyv önen 2003a, Liski et al. 2005).The position of the wood has a direct influence on the local microhabitat for the decomposers and affects their species' abundance and distribution (Parisi et al. 2018).The diameter of the decomposing material (Hyv önen et al. 1998), wood properties of different species (Djukic et al. 2018, Wang et al. 2018), and CWD position (Mäkinen et al. 2006) all affect decomposition kinetics.Neglecting these interactions can reduce the accuracy of predictions for deadwood pools, as Didion et al. (2014) found for the Yasso07 model.Differentiating the decomposition of CWD would help to develop harvesting strategies based on the specific decomposition rates of the various CWD classes.
On top of being useful for drawing conclusions, the main point of a parametric model is the possibility of using it for extrapolating quantitative predictions.For example, answering the question "how much faster," rather than just a yes/no (as with a statistical approach, also a simple case of the model), and hence opening up for the possibility of comparing numerically different management alternatives.

Aims of this study
For merely testing a hypothesis (falsifying it with a yes or no), a statistical approach would still be definitely appropriate and robust and it has been applied in the past to understand the impact of tree species and wood position on the decay of CWD in forests.The aim of this study was to assess such impact quantitatively developing a model that can be then used for further extrapolation.We utilized the Q model (Hyvönen et al. 1998) and calibrated it on two datasets (Tarasov andBirdsey 2001, Mäkinen et al. 2006).The data by Mäkinen et al. (2006) come from a set of longterm thinning experiments in Finland and the data by Tarasov and Birdsey (2001) from a Russian chronosequence study.The Finnish dataset contains information on different decomposition stages of CWD, tree species, namely Scots pine (Pinus sylvestris L.), Norway spruce (Picea abies (L.) Karst.), and silver birch (Betula pendula Roth.), and the position of the material with respect to the ground.The Russian dataset provides information about different decomposition stages for spruce.We then compared two calibrations of the model, one considering all this information for separate classes and one considering the same kinetic for all classes, within a Bayesian framework (Carpenter et al. 2017).The overall intention was to contribute to the development of new forest management strategies enabling timber harvest while optimizing the amount of ecological ecosystems services from the CWD left on site.We developed a model calibration considering CWD classes that can be used as a tool for developing such management strategies.

The calibration dataset
The data used in the calibration consisted of deadwood data obtained from Finnish thinning experiments for birch, pine, and spruce (Mäkinen et al. 2006), integrated with spruce data from Russia (Tarasov and Birdsey 2001;Fig. 1).
The Finnish dataset originated from 58 longterm thinning experiments in southern and central Finland was established by the Finnish Forest Research Institute in the 1960s and 1970s and a few older experiments.The stands were even-aged and were mostly pure Scots pine, Norway spruce, or silver birch.At the time of the sampling, the stand age ranged from 29 to 131 yr.The dataset contained information about: wood type (standing:snag, lying:log), tree species, time since death, initial and final wood density, initial diameter, and geographical coordinates of each stand (of which we utilized latitude).The Russian dataset contains similar information, but it is based on chronosequence studies, so wood parts that were fully decomposed (and therefore disappeared at the time of sampling) were not identified and recorded.We, therefore, considered it as a separate class.
In total, the datasets contained 1386 data points (Fig. 1).We reclassified them into coarser classes compared with the more detailed classification in Mäkinen et al. (2006) (which distinguishes various types of logs) and kept only a distinction between logs, snags where <1/3 of the whole length was broken and snags where ≥1/3 of the whole length was broken.These classes were combined with the tree species to produce 10 CWD decomposition classes (Fig. 1): 1. Logs, snags, and snags with <1/3 of the length broken from the pine.2. Logs, snags, and snags with <1/3 of the length broken from the spruce.3. Logs, snags, and snags with <1/3 of the length broken from the birch.4. The spruce data from Tarasov and Birdsey (2001), which were considered separately.

The model
The model ( Ågren andBosatta 1998, Hyv önen et al. 1998) considers decomposition of organic matter (OM) based on the assumption that microorganisms are C-limited and their growth and death rates depend on the availability of C sources and the assumption that C exists for decomposers in different availability states over a continuous quality gradient.Quality is proportional to the distribution of the mean residence time of all organic compounds in the decomposing material.The model describes how the average distribution of quality of the organic C in a system changes as decomposition proceeds (while part of the C is lost as respiration).This is a generalization of the theory on which all compartmental C models are based (Coleman 1996, Liski et al. 2005) but does not require specific C "pools" to be defined.The dynamic change in the average C quality is instead used to represent this.
We integrated a Bayesian error model into our predictive model.In the subsequent text, the term "priors" denotes the prior probability distributions while the term "posteriors," the posterior probability distributions (Carpenter et al. 2017).
The model describes the decomposition of OM controlled by a set of parameters.First, decomposers are described by a certain production to assimilation efficiency, e q , that determines the fraction of C that goes into new decomposer biomass of all the C utilized.This term is conceptually equivalent to carbon use efficiency (CUE) and can be compared to it.However, the CUE concept is in practice strongly dependent on the measurement used and its definition can be ambiguous (Geyer et al. 2016); so, the correspondence here is rather loose.For simplicity, we assumed that efficiency is independent of the quality of the C assimilated ( Ågren and Bosatta 1996).The term was extended as a normal probability distribution based on that suggested by Manzoni et al. (2012), e q,p = 0.30, to which we added a standard deviation σ e q;p ¼ e q;p Á 0:25: e q ∼ N e q;p , σ e q;p : (1) We next assumed that the decomposer growth rate per unit of C, u q , depends on the quality of the C used, and we assumed a power function of a standard decomposer growth rate u 0 as in Hyv önen et al. (1998).The term u 0 was calculated as a function of latitude as: u 0 = 0.0855 + 0.0157Á (50.6 − 0.768ÁLat) ( Ågren and Hyv önen 2003b).We extended this calculated term with a probabilistic error term ε u 0 of mean ε u 0 ;p ¼ 1 and coefficient of variation c v ¼ ðσ ε u 0 ;p =μ ε u 0 ;p Þ ¼ 0:5: The parameter β is a shape parameter determining how fast the decomposer growth rate changes with quality.It can be related to edaphic properties, e.g., the soil texture (Bosatta and Ågren 1997).This term was also expressed probabilistically as a function of former value β p = 7 and coefficient of variation c v ¼ ðσ β p =μ β p Þ ¼ 0:1: β ∼ N β p , σ β p : (3) Microbial biomass dies and then the OM returns to the substrate but with lower quality, according to η 1,q , extended statistically as a normal probability distribution based on the former value from Hyv önen et al. (1998), η 1,p = 0.36, plus a coefficient of variation c v ¼ ðσ η 1;p =μ η 1;p Þ ¼ 0:1: The decomposition model is defined for one single litter cohort as: (5 where t is the decomposition time and α is defined as: and z is defined as: This single cohort model can be used to describe one single pulse of organic material (e.g., one log), but it can also be used to describe a situation of constant inputs by overlapping multiple decomposing cohorts (each corresponding to one pulse of inputs).The initial quality of the organic material q 0 was defined probabilistically here, with average and deviation calculated from Joffre et al. (2001).We utilized as average q 0,p = 1.1 and as deviation σ q 0;p ¼ q 0;p Á 0:12: The deviation for this prior was set slightly larger than for the other terms because we considered this term to be likely more variable than the others due to the variation in the origin of the material.Since we simulated only one OM cohort, represented by the deadwood material, we did not need to integrate for the OM flux at every time step as described in Hyvönen et al. (1998).
The biomass is converted into the mass of carbon through the term f c , expressing the C content of the biomass and also introduced as a normal prior: We utilized as average f c,p = 0.5 and as deviation σ f c;p ¼ f c;p Á 0:1.

Representing delayed decomposition: introducing Tmax and Δ
We introduced the parameter T max to represent the initial diameter of the decaying material, and it influences the time needed for the colonizer microorganism to invade the wood.The model was originally based on the deduction of T max (Hyv önen et al. 1998), but since in our dataset we had the diameter measured at time 0, we utilized this value.Since we had data on measured wood density ρ at time 0, which also influences the speed at which the microorganisms colonize the wood, we used it in the model as a linear modifier of T max .To do so, we standardized the wood density data between 0.8 and 1.2 as ρ and utilized this value to directly rescale T max based on the wood density.
Tmax;0 is a deterministic term that was extended with a probabilistic error term ε T max of mean ε T max ;p ¼ 1 and coefficient of variation v www.esajournals.org In addition, we introduced in the model a delay term Δ (expressed in years), which simply delays the decomposition of the wood (Appendix S1: Fig. S7).The term Δ was extended as a uniform probability distribution between 0 and 5 yr: The term Δ is similar to T max , as they both slow down the initial phase of decomposition, but is conceptually different: T max represents a physical property of the object decaying, while Δ is related to the position of the object (lying, standing, etc.).An increase in T max generates a gentle bend of the initial shoulder of the decomposition curve, as in Appendix S1: Fig. S6, but an increase in Δ just postpones the start of the decomposition.

The calibration procedure
We performed the parameter estimation within a Markov chain Monte Carlo (MCMC) framework, implementing Bayesian statistical principles (Kruschke et al. 2012) to update the current knowledge of the model uncertainty with the information in the dataset.The method relied on statistics performed on a population of multiple model realizations, starting from our initial probabilistic knowledge of the model parameters (from now on priors), combined with the knowledge from the data to produce an updated knowledge of the probability distributions (from now on posteriors).The priors can be defined by any probability distribution derived from various sources of information.This allowed us to evaluate all the information in the data with continuous non-parametric probability distributions, retaining all the available information and synthesizing it into posterior probability distributions of predictions and parameters.
The model was written in the language Stan (Carpenter et al. 2017), a relatively recent Bayesian framework that bases the calibration on Hamiltonian MC sampling.This sampling strategy is more effective for complex spaces, as it converges faster and reduces autocorrelation compared with conventional Metropolis-Hastings samplers (Carpenter et al. 2017).The Q model, like most deterministic or semi-deterministic ecological models, tends to have more parameters than the available data can resolve (Marschmann et al. 2019).The model, therefore, displays relatively high equifinality (Beven and Binley 2014), and exploring its space can benefit greatly from this improvement in efficiency.
We performed numerical analyses and plotting from R (R Development Core Team 2019), including running the Stan model through Rstan (Stan Development Team 2019).The model fitness was evaluated as a root mean squared error (RMSE).
We ran all calibrations in independent quadruplets of MCMCs of 10,000 runs in each chain.The model gives unfeasible parameter combinations whenever z ≥ 1, as this leads to the infinite accumulation of C. Thus, we instructed the sampler to reject the proposal whenever z ≥ 1.The posterior distributions were, therefore, not influenced by unfeasible model combinations.

Parameterization diagnostics
The calibration process was monitored through a few diagnostic tools.First, we checked the trace plots of the MCMCs for selected parameters (Appendix S1: Fig. S1).The main assumption of the MCMC calibration is that the probability of each event in the chain (in our case, a model run with a certain parameter combination) depends only on the state attained in the previous node.This means that the sampling algorithm should change the value of each parameter very often, without getting stuck in local optima.In our case, all chains seemed to mix well.We also checked the auto-correlation of the MCMCs, which decreased with the increasing estimation window (Appendix S1: Fig. S2).This meant that for an increase in the number of samples, we converged toward the perfect realization of the Markov chain assumption.
We then examined the R of the whole calibration.This is a convergence criterion based on the use of multiple independent chains and measures the ratio of the average variance of each chain to the variance of the pooled draws across Markov chains (Gelman and Rubin 1992).The R value can be calculated when more than one fully independent chain are calibrated.When the sampling goes as it should and when the mixing is even, the value should converge toward 1.The conventional limit for well-mixed chains is R ≤ 1:05 (Gelman et al. 2004, Carpenter et al. 2017).In our case, R was in general below this v www.esajournals.orgthreshold (Appendix S1: Fig. S3), indicating good mixing.

Global and specific calibrations
We calibrated the model by the two separate approaches: 1. Global calibration: One single model parameterization was used to predict all the values.
No parameter was left free to vary with specific wood characteristics (i.e., different CWD classes), but one single set of values was used to predict all the data points.All the parameters were calibrated simultaneously on all the observations (Fig. 2).2. Specific calibration: The parameters q 0 , Tmax , Δ, and u q were left free to vary for specific CWD classes, while other parameters were still the same for all classes (Table 1).Specifically, Δ was assumed to vary with wood being a snag or a log, q 0 was assumed to vary depending on the tree species, and the other two parameters were assumed to vary according to the decomposition classes.The data used were the same as in the global calibration approach, and all the parameters were calibrated simultaneously on all the observations.
Global calibration: one single set of parameters and sensitivity analysis.-Wefirst assessed the model fitness response to the variation in different parameters, to gain an indication of which parameter to focus on the most.We ran an MCMC model calibration considering only one generic set of unknown parameters to describe all the data, producing a set of 5000 model runs.We then utilized this set of model runs to perform a Hornberg-Spear-Young sensitivity analysis (Beven 2008, Beven andBinley 2014).This involved splitting the population of model runs into two separate bins and then, for each parameter, assessing the difference between the distribution functions of each parameter between the two bins.
We considered the first bin ("behavioral" models Beven 2008) composed of all the parameter sets with the RMSE falling within the first five percentiles, and the second bin ("non-behavioral") composed of all the other parameters sets.We assessed the distance between the probability distributions of a certain parameter value in each bin with the Kolmogorov-Smirnov distance (Marsaglia and Wang 2003).The more the two distributions differed, the more that specific parameter was considered to contribute to the increase in fitness observed in the upper five percentiles, and, therefore, the more the model fitness was considered to depend on that parameter.The analysis indicated that the effect of the parameters was strongly dependent on the different groups (Fig. 3; Appendix S1: Fig. S5).These results were then considered to select which parameters to calibrate specifically by class.
Specific calibration: some parameters specific for different CWD classes.-In the specific calibration, we also considered 5000 model runs but let some parameters vary for different CWD classes, partly according to the sensitivity analysis and partly based on our hypotheses on processes influenced by the CWD classes.Our aim was to explain most of the variance in the dataset according to a specific decomposition parameter, and to investigate which parameters might explain the observed variation.
We considered the term q 0 dependent on tree species (birch, spruce, pine, and the spruce data from Tarasov and Birdsey [2001]).The term was defined according to Eq. 9, but with the error term dependent on the four different species: We then considered ε Tmax as dependent on decomposition classes, as this parameter summarizes the delay in decomposition due to different effects (e.g., contact with soil).The term was defined according to Eq. 10, but with the error The columns "Global" and "Specific" refers to the two calibrations, and "GL" means that the parameter is considered global, and "SP" means it is considered specific.μ is the mean and σ is the standard deviation. † We assumed also that the Δ term was dependent on whether the CWD type was a snag or a log.
Finally, we considered the term u q , which represents the growth rate of the decomposers, as specific and dependent on decomposition class.The term was defined according to Eq. 9, but with the error term dependent on the 10 different CWD classes:

Steady states
The model can be solved analytically by setting the derivatives to zero to find the C stock in an ecosystem at a steady state, assuming constant inputs: where C SS is the C stocks at steady state and I is the rate of input in terms of C (Hyv önen and Ågren 2001).Determining the steady state, although usually reached over centuries or millennia (since the system approaches the steady state, the slower the closer it gets to it), is important to test the consequences of the decay of the different organic material in a comparable way without having to assume an arbitrary time frame.Steady states are not influenced by time while, when comparing different materials or input levels, comparison results might change depending on the time perspective considered.A dynamic simulation of C stocks requires a known amount of C inputs.To compare the effect of the different organic materials in terms of C stocks at equilibrium, we assumed inputs equal to 1 MgÁha −1 Áyr −1 , a number easy to handle but plausible as possible scenario (similar values were measured for natural litterfall in unmanaged forests, Bray and Gorham 1964).These levels of input are too high for managed forests but permitted a relative comparison between the different types of CWD.The steady state is an integration of the decay represented by Eq. 5 over (very long) time, and, therefore, small differences in decomposition end up accumulating over time.

Global calibration and sensitivity analysis
The global calibration, with only one set of parameter values, could not explain much of the variance (Fig. 2).The dataset displayed a scattered cloud of observations, indicating that the decomposition is influenced by many interacting processes not considered by the global parameterization of the model.For some CWD classes, it was possible to distinguish a pattern, e.g., birchwood decomposed faster than the other tree species.However, for most CWD classes, there were v www.esajournals.orgno easily identifiable patterns, in particular concerning points with delayed decomposition which were not explained by the information in the data.
The sensitivity analysis (Fig. 3; Appendix S1: Fig. S5) identified the parameter ε Tmax , controlling the initial shoulder of the curve, as particularly important, together with the initial wood quality q 0 .Different CWD classes showed differing sensitivity to model parameters (Appendix S1: Fig. S5), with snags of spruce and pine being the most sensitive, while the fast-decaying classes (birch, spruce data from Tarasov and Birdsey [2001], and pine logs) were, in general, less sensitive.
The posterior of parameter e 0 did not differ much from the prior, since this was indicated as centered around 0.3 from Manzoni et al. (2012).

Specific calibration
The specific calibration also assumed that some parameters were common to all observations (β, η 11 , e 0 , and f c ).The posteriors for these global parameters (Table 2, Fig. 4) seemed rather similar to the priors or, in the case of e 0 , even better determined (Fig. 4).
The local parameters expressed the variance in the various CWD classes.The parameter q 0 varied with the tree species (plus the Russian spruce dataset, Fig. 5, panel A).Pine and spruce showed the lowest initial substrate quality (and, thus, slower initial decomposition), while birch showed the highest.The Russian spruce data presented high variance, but q 0 seemed in general higher than for the Finnish spruces.The Δ (delay) parameter was quite close to zero for both snags and logs (Fig. 5).The parameter Tmax , which referred to the initial diameter, had instead higher explanatory power (Fig. 5), although referring to a similar process to Δ (since  v www.esajournals.orgthey both express an initial slowdown in the decomposition rate).The posterior probability distribution for the delayed decomposition of logs was slightly higher, but still on average below 1 yr (Table 3).The error parameter ε T max varied according to the 10 CWD decomposition classes (Fig. 5, panel C).Most classes had an average ε T max value of less than 1 (Tables 2, 3), except the Russian spruce data.Short pine snags and long birch snags had the lowest Tmax .Moreover, the error parameter ε u 0 varied according to the 10 decomposition classes (Fig. 5, panel D) but was generally above 1.It was clearly below 1 only for short snags of spruce and pine (Tables 2, 3).

Calibration diagnostics
The specific calibration, being at the same time the most accurate and the main product of the study, was tested with diagnostic criteria to ensure its robustness.All diagnostic values indicated a good calibration procedure.The calibration indicated excellent mixing of the chains (Appendix S1: Fig. S2) and consequent low autocorrelation (Appendix S1: Fig. S2).The Gelman-Rubin potential scale reduction statistic R (Gelman et al. 2004) was also good for all the chains (Appendix S1: Fig. S3).

Simulation results and steady states
The parameterization generated by the global calibration missed many of the points in the datasets (Fig. 2).The local calibration captured a larger share of the variance, although some classes did not show a relevant improvement in fitness (Appendix S1: Fig. S4).Pine and spruce exhibited similar trends over time and similar parameters (Figs. 5, 6), and for both species, the pattern caused by the decomposition class was the same, with logs decaying quickly and standing snags decaying more slowly.Birch showed a similar pattern, but with faster rates.The Russian spruce data had different kinetics from the Finnish spruce data, decomposing much faster.This might be due to the different determination of time since death in these two datasets.
The different CWD decomposition classes showed clearly distinguishable kinetics (Figs. 5,6), in particular depending on tree species but also between logs or snags.The locally calibrated model was more accurate, with an RMSE of 0.153, compared with 0.169 for the globally calibrated model (also Appendix S1: Fig. S4).
The difference in the predicted remaining C between the general and local calibrations was large, with an average absolute difference of 11% at 50 yr, 10% at 80 yr, and 15% at 120 yr.These differences were larger for some specific decomposition classes, in particular for spruce snags (Appendix S1: Fig. S4).The difference was relatively large also for logs (Fig. 7).This resulted in over-or underprediction depending on the specific CWD class, while the general calibration was mostly over-predicting.
The differences in the C stocks accumulated at the steady state were substantial between the slow and fast decaying materials (Table 4).In general, the steady states were characterized by considerable uncertainty and wide boundaries, but the uncertainty was nevertheless considerably less than that during the initial phases of decomposition.While approaching the steady state, the uncertainty of all simulations decreased (Fig. 6).This seems to be due to fewer parameters influencing the steady states, i.e., Eq. 17 has considerably few parameters than the dynamic version of the model (Hyv önen et al. 1998).In particular, the OM quality converges toward similar values over time (Menichetti et al. 2019).

The model and its two parameterizations
The effects of tree species and wood position on CWD decomposition were studied by extending and parameterizing the Q decomposition model.The posteriors for the global parameters did not differ much from the values previously reported in the literature (Hyv önen et al. 1998) and from the priors.The global parameterization of the model in most cases resulted in lower C stocks at steady states compared with the specific parameterizations.The specific calibrations predicted different C stocks for each CWD class for both the steady states and predictions over time, and some cases were much lower (mostly birch) but other cases v www.esajournals.orgmuch higher than in the global calibration of the model.This indicates that there could be large errors in calculating the stand C balances if CWD classes are not considered separately, in particular, an underprediction of the C stocks for spruces and pines.This would be particularly critical for processes dependent on specific CWD classes, e.g., beetle colonization.The model we selected as a starting point represented the decomposition delay with Tmax .In our model, for calibration, we added two additional decomposition modifiers, and the term Δ is relatively similar to the term ε T max .Both have delayed the decomposition in the initial phase, but the former postpones decomposition for some time, while the latter describes an initial slowing down of decomposition (Appendix S1: Figs.S6, S7).The slow initial decomposition recorded in some of the data points was represented more often in the calibration with Tmax than with the delay term Δ, although both terms describe similarly a slowing down of the decomposition process.The slow decomposition in the early phases after tree death was also variable within the CWD classes.There are many possible causes for the slow decomposition rate and, therefore, the considerations about potential affecting factors remain speculative.One possible contributor not accounted for by the available classifiers is the sapwood/heartwood ratio, which is variable between tree species but also among individual trees.
The parameter u 0 , describing the decomposer growth rate, influenced the initial part of the curve most strongly (Appendix S1: Fig. S6), and the parameter value was higher for wood characterized by faster decomposition.This might be related to microbial community effects, with some substrates (e.g., different availability of nitrogen) influencing the general decomposers' growth rate or even favoring a particular microbial class (Bani et al. 2018).The parameter q 0 reflects the differences between the tree species, in particular between the birch and both spruce and pine.The much faster decomposition of birch compared with pine and spruce seemed to be due to intrinsic wood characteristics.Although with broad variation, the previous studies have reported comparable differences in chemical wood properties between tree species (Kahl et al. 2017).For example, Weedon et al. (2009) found large differences in decomposition kinetics related to tree species.Accordingly, in a laboratory study, (Lasota et al. 2018) were found differences in the release of organics from the wood of eight different tree species (including spruce and birch).All these results suggest that tree species have an impact on wood decomposability (represented in our model by the initial quality q 0 ).
The posterior distributions of the common parameters of the local parameterization (Fig. 4) were symmetrical.This means that the unexplained variance is normally distributed and that there are no regular biases in the model errors, which in turn suggests that the degrees of freedom introduced with the local parameter allowed the model to represent most of the key processes involved in decomposition.

Impact of CWD characteristics on decomposition
Coarse woody debris is usually left in the stand until final felling when it is often harvested, but it can be harvested earlier for utilization (Brenøe andKofman 1990, Kärhä 2011).The large variation in wood decomposition rates is heavily affected by tree species and by contact between the wood and the forest floor (Zhou et al. 2007).Once we introduced these factors in the model, the prediction error was reduced, confirming that the wood characteristics (species and position) affect both long-and short-term decomposition kinetics.This also caused differences in the predictions for longterm C stocks.
A clear effect of the position of the wood on their decomposition has also been reported.Logs lying on the ground are more easily attacked by microorganisms, decomposing faster than snags (Mäkinen et al. 2006, Harmon et al. 2020).This can be due to faster colonization by fungal hyphae, but also the local microclimate, e.g., higher moisture when the organic material is in contact with the ground (Joly et al. 2017).
The data points were scattered in all our CWD decomposition classes (Fig. 6), suggesting an impact of some local conditions or processes.This could be, e.g., a fallen tree broken at the base, but which remains for some time hanging against other trees, delaying decomposition, or isolated from the ground by the roots at the base.Our approach compared with current SOC models Current forest decomposition models do not represent the influence of local microenvironments on decomposition kinetics.In general, current CWD decomposition models, while performing well for lying organic material, tend to overestimate the decomposition of standing dead trees (Didion et al. 2014).For example, considering micro-environments, defined by temperature and moisture and tree species, affecting CWD decomposition (Zhou et al. 2007, Herrmann et al. 2015) impacts the accuracy of C balance calculations (Harmon et al. 2020).
The impact of tree species on wood decomposition is represented in other models by defining the initial quality by tree species, which in compartmental models correspond to the proportions between the pools, while in Q, the initial quality is represented by a single parameter.This information is not easily available.For example, the Yasso model (Tuomi et al. 2011) is based on the decomposition classes defined from the size of the wood residual and determined by analyzing the chemical properties of wood.The classes only represent the chemical recalcitrance of the material, which might neglect other effects related to the tree species (such as porosity or bark structure) and local parameters were considered specific to each class) by classes (minimum and maximum are calculated based on all the MCMC population).The shaded area represents the maximum and minimum of the prediction interval.
(Fig. 6.Continued) v www.esajournals.orgmicro-environments.The Century model, one of the most used compartmental decomposition models for OM and in use also for forests, has a similar approach, predicting decomposability based on the lignin content of the substrate (Parton 1996).Both models, as well as the current version of the Q model, do not consider the impact of local microclimate.More mechanistic models such as ROMUL (Chertov et al. 2001), which represents fauna effects in detail, also neglect local microenvironments (Komarov et al. 2017).The same applies to the forest DNDC (Kurbatova et al. 2008).Applications of the detailed Coup model (Svensson et al. 2008) consider kinetic modifiers only at the stand level.All these models (including Q) could be tuned with kinetic modifiers dependent on local conditions and calculated according to the parametric functions of some local ecological variable, but this would increase the complexity of the model and require local fine-scale data to run the model.
Our modeling approach, using calibrated parameters linked to specific CWD classes, contains most of the necessary information and is easier to apply.The calibration of the parameters on empirical data is more robust than in an approach based on wood chemistry, as it considers effects described by the information in the data, but which are not necessarily just chemical.The inclusion of these effects in models is important because it represents the impact of management (e.g., how wood residuals are treated, if left in place or piled, and how) on forest C balance.Moreover, the utilization of specific parameters linked to decomposition classes is simple to implement, since it does not require additional functions to model complex ecological interactions.The approach we propose can also be directly incorporated in other wood decomposition models by utilizing the decomposition equation (Eq.5).
The impact of decomposition kinetics on forest C stocks in the long term (decades) and very long term (millennia) Managed forests do not present "steady states," since management causes them to fluctuate around an average C stock.Steady states are still useful to apply to managed forests by indicating relative differences between the treatments independently of time.Steady states are conceptually similar to the mean C stock that a stand with certain management accumulates over a very long time.At a steady state, our model predicted 4-12 times more OM accumulation for spruce and pine than for birch.In general, steady states were characterized by considerable uncertainty, with wide upper and lower boundaries.In particular, microbial parameters, such as efficiency (e 0 ), lack precise quantitative estimates in the literature, and even more so β and η 11 , which are estimated by the inverse modeling.These parameters impact the steady state uncertainty (Hyvönen et al. 1998), and in our study, this impact was particularly large because of the lack of reliable data on the latter parts of the decomposition curves (Harmon et al. 2020).
The absolute C stocks we calculated depended also on the input levels, varying between the management regimes and biomes, but the relative differences between the management regimes provide indications about their influence on the C balance of a stand.If the aim is the extraction of forest energy while still increasing CWD stocks stored in the forest, birch wood is more suitable to be removed than coniferous species since it decays faster.Pine and spruce residuals and logs left in forests have a higher potential for C storage, in particular, if left standing.So, for the same mass of wood extracted from a mixed-species forest, the higher C accumulation is achieved by extracting the birch and leaving coniferous trees on site.The uncertainty we detected in the estimation of the long-term C stocks stresses the importance of longer-term decomposition observations.Without decomposition time series spanning 5-10 decades, it is difficult to constrain the model for steady states, which are typically reached after thousands of years (Hyv önen et al. 1998).Until longer-term data about the consequence of specific management become available, the only way to extrapolate from short-term data is through modeling.For this reason, representing all key processes influencing the decomposition kinetics in decomposition models is crucial when simulating the long-term effects of forest management approaches.

CONCLUSIONS
Calibration of the model based on global parameters generated less accurate predictions and explained less variance than calibration specific for the CWD classes, which captured differences in decomposition caused by CWD position and tree species.The amount of C remaining predicted after 50 yr for different classes with local parameterization was highly different from that predicted by the global parameterization (Table 5).The tree species had different decomposition rates, and the position of wood (lying on the ground or standing) affected the decomposition delay (supposedly proportional to the time needed for microbial colonization) and, thus, the forest C balance.
Considering CWD classes separately reduced the uncertainty of the predictions, but still left some unexplained variance in the decomposition kinetics within each class due to some other processes not considered in our data or model.To improve the reliability of CWD decomposition models, there is still a need for more detailed studies considering more ecological covariates, and for updates to monitoring programs and methods (e.g., the introduction of proximal sensing via three-dimensional laser scans and automatic monitoring stations).Such data are needed to explore proxies of decomposition that could be used to map and predict the unexplained variance in the kinetics.
The differences in decomposition associated with CWD classes that our model predicted affect the C balance of a stand already over a few decades.The amount of CWD in a stand can be influenced by forest management, either by harvesting certain classes or by retaining specific slowly decomposing CWD classes, and our model can be utilized to predict how much.In the future, extraction of biomass will become increasingly important, increasing the pressure on CWD harvest.Therefore, reliable modeling of the decomposition of different CWD classes will be urgently needed to decide what to harvest and not harvest, or even how to harvest (e.g., with the alternative thinning techniques such as girdling Riffell et al. 2011).Our model calibration is potentially useful for developing biomass harvesting strategies involving thinning residuals and other CWD and including them in a precise quantitative assessment of a stand C balance.This tool can be used together with highresolution data on wood biomass, such as the ones from low-altitude laser exploration (at the moment, an experimental technique in active development), to calculate a precise C balance of vast areas of the forest.In addition, the decomposition of different CWD classes left in forests can be better taken into account when considering their role in maintaining biodiversity.

Fig. 1 .
Fig. 1.In this plot, each colored point represents one data point, classified according to the 10 CWD classes considered in this study.Numbers in the legend are the number of points in each class.

Fig. 2 .
Fig. 2. The dataset plotted for the predictions of the fraction of remaining mass over time from the global model calibration (where all parameters were considered the same for all the classes).The shaded area represents the maximum and minimum of the prediction interval.

Fig. 3 .
Fig. 3.The results from the Hornberger-Spear-Young sensitivity analysis based on the global calibration (where all parameters were considered the same for all the classes), averaged for the three coarse woody debris positions in the dataset.The sensitivity measurement represents the Kolmogorov-Smirnov distance between the two bins as described in the text.

Fig. 4 .
Fig.4.The result of the specific calibration (where some parameters were considered specific to each class) for the global parameters (considered the same for all the classes also in the specific calibration).The grey dashed area represents the prior probability distribution, while the green solid area represents the posterior probability distribution.

Fig. 5 .
Fig. 5.The result of the specific calibration (where some parameters were considered specific to each class) for the local parameters (the parameters considered to vary among classes).The grey dashed area represents the prior probability distribution, while the colored solid areas represent the posterior probability distribution for each decomposition class.

v
Fig. 6.The predictions of remaining mass as a function of time from the specific calibration (where some

Fig. 7 .
Fig. 7.The difference in model predictions between the global calibration (where all parameters were considered the same for all the classes) and specific calibrations (where some parameters were considered specific for classes) for three time points (50, 85, and 120 yr).Error bars represent the standard deviation calculated over the whole MCMC.Values refer to the difference in the estimated remaining C concentration between the two calibrations (global-specific).

Table 1 .
Summary of the parameters considered in this study.

Table 2 .
Results of the specific parameterization on the common parameters.
Note: Values are calculated on all four MCMC.

Table 3 .
Results of the specific parameterization on locally calibrated parameters.
Note: Values are calculated on all four MCMC.

Table 4 .
Steady states (C stocks, Mg/ha) calculated assuming 1 MgÁha −1 Áyr −1 of input for each class, with uncertainty boundaries (95% CI, min = lower bound, max = upper bound.Values calculated on all four MCMC.Due to the nonparametric nature of the error estimation, boundaries are asymmetric).

Table 5 .
Model predictions (means) after 50, 80, and 120 yr, expressed as the fraction of the initial mass.