Bayesian inference of biomass growth characteristics for sugi (C. japonica) and hinoki (C. obtusa) forests in self-thinned and managed stands

Forests are an important sink for atmospheric carbon and could release that carbon upon deforestation and degradation. Knowing stand biomass dynamic of evergreen forests has become necessary to improve current biomass production models. The different growth processes of managed forests compared to self-managed forests imply an adaptation of biomass prediction models. In this paper we model through three models the biomass growth of two tree species (Japanese cedar, Japanese cypress) at stand level whether they are managed or not (self-thinning). One of them is named self-thinned model which uses a specific self-thinning parameter α and adapted to self-managed forests and an other model is named thinned model adapted to managed forests. The latter is compared to a Mitscherlich model. The self-thinned model takes into account the light competition between trees relying on easily observable parameters (e.g. stand density). A Bayesian inference was carried out to determine parameters values according to a large database collected. In managed forest, Bayesian inference results showed obviously a lack of identifiability of Mitscherlich model parameters and a strong evidence for the thinned model in comparison to Mitscherlich model. In self-thinning forest, the results of Bayesian inference are in accordance with the self-thinning 3/2 rule (α=1.4). Structural dependence between stand density and stand yield in self-thinned model allows to qualifying the expression of biological time as a function of physical time and better qualify growth and mortality rate. Relative mortality rate is 2.5 times more important than relative growth rate after about 40 years old. Stand density and stand yield can be expressed as function of biological time, showing that yield is independent of initial density. This paper addressed stand biomass dynamic models of evergreen forests in order to improve biomass growth dynamic assessment at regional scale relying on easily observable parameters. These models can be used to dynamically estimate forest biomass and more generally estimate the carbon balance and could contribute to a better understanding of climate change factors.


Background
Forests are an important sink for atmospheric carbon and could release that carbon upon deforestation and degradation (Lim et al. 2013). Therefore the dynamics of biomass growth are important to study especially in relation to future environmental changes (Fang et al. 2014;Lehtonen et al. 2004). The forest ecosystem can also become a source of contamination for several ecosystems in the context of post-accident (watershed, leaching leading to contamination in aquatic and agricultural environments (Fukuyama et al. 2005;Karki and Shibano 2006), groundwater contamination (Bugai et al. 1996), fire problem reemitting particles or gases in the atmosphere (Evangeliou et al. 2016), ...).
Two-thirds of Japan's land area is covered with forests, with a total forested area of 25 million hectares. Approximately 40% of these forests are artificially planted forests and the major planted species are Sugi (cedar), Hinoki (cypress), and Karamatsu (larch) (Forestry 2012). In Fukushima prefecture there are 83 cedar forests and 42 cypress forests from the 190 private forests, from which 16 cedar forests from the 25 public forests (Ministry of Agriculture 2012). This paper aims to extract an application methodology for an even-aged, mono specific population corresponding to the cases of stands used for timber. We focused on coniferous stands present in Japan (on Japanese cedar (Sugi) and Japanese cypress (Hinoki)) (Yoshihara et al. 2014;Yoshihara et al. 2016;Ogawa et al. 2016;Aoki et al. 2017).
The determination of the biomass growth driving at the forest stand level animated different researchers and modelers (Bartelink et al. 1997;Canham et al. 2004;Will et al. 2002). Several studies reporting the construction of yield tables were provided in literature, often based on a diameter growth model (Matsumoto 1997;Nagahama and Kondoh 2006;Shimada and Mie-ken 2010;Nakajima et al. 2014).
The notion of competition between trees within a population density to maximize their access to resources (light, soil nutrients, ...) is an essential factor in stand growth. The initial density of a forest stand has a large impact on the yield of the logs in the plot. Relationships between yield and stand density will naturally result from tree growth assumptions.
Moreover, a forest can evolve either in autonomy that means without human intervention on its production of biomass (e.g. by pruning, cutting thinning, weeding and cutting the undergrowth species before closing the canopy, ...), it is then said to be auto-managed, or it is managed by foresters. When forest is self-managed, a selfthinning phase occurs naturally. This phenomenon allows a stand to self-regulate and leads to a naturally optimized production of biomass. We then obtain so-called self-thinning lines in a log-log plot of stand density yield that indicate maximum forest yield. The method to determine the maximum density, line of self-thinning was discussed and a self-thinning model was developed using observed long-term data in (Fengri and Lingbin 1995). When the forest is harvested by foresters, thinning is scheduled regularly to achieve higher returns to the short term. A forest growth model is an abstraction of the natural dynamics of the tree or stand that may include growth, mortality, and other changes in stand composition and structure. Classically predictive forest growth models are of three types: • Population model which are characterized by global variables of the forest stand (stand density, basal area of the stand); • Distribution model which refers to particular species or classes of trees within a forest and uses relationships based on the driving dendrometric variables such as diameter at breast height (DBH), dominant height, ...; • Individual-tree model which consists of predicting stand growth and yield using individual trees as a base unit.
These different types of growth models position themselves in relation to each other in terms of the resolution of the results (Cao 2014). An individual-tree model is finer at the tree level but leads to stand-level errors. Conversely, a population model provides information on the whole population but has a low resolution. The population models respond more to the problem of prediction over huge spatial domains. Thus the notion of stand density as biomass growth explanatory variable answers this need of granularity, and allows to overcome the heterogeneity of smaller metric scales (individual level).
This article also aims to verify whether the process leading to the dependence of mortality and self-thinning on population density can advantageously explain the nonmanaged forest stands observed data through the application of an empirical self-thinned model. For managed forest, we compare two types of empirical model by a Bayesian approach with a likelihood criteria: one called thinned model and the other one was a Mitscherlich model.
To achieve this goal, a Bayesian inference was carried out with the software SAMCAT (Nicoulaud-Gouin et al. 2016) to determine parameters values according to large collected data.

Methods
In this section we present in detail the forest growth models we have chosen according to its resolution granularity corresponding to our objectives.
Non managed forests are linked to the competition effect by the density of their population as we will see in the next Competition-density effect paragraph.
We describe the stand density, the stand yield and its compartments biomass (branch, leaf and stem) at a time t according to known characteristics in Density, yield and biomass paragraph.
The data will be analyzed in the paragraph Data and the Bayesian inference method allowing to calibrate the different parameters of the models on the data is detailed in paragraph Bayesian inference and model validation.

Competition-density effect
Self-thinning in overcrowded pure stands under cultivated and natural conditions was studied (Yoda et al. 1963;Shinozaki and Kira 1956), and it was found that for each initial stand density a corresponding rate of maximal productivity was associated. That means that competition among trees and self-thinning were the most important ecological process in tree population dynamics. Authors derived the reciprocal equation of the competition-density (C-D) effect (Eq. 1) (Shinozaki and Kira 1956) for even-aged pure plant populations using the general logistic growth equation and the 3/2 power law of self-thinning (Yoda et al. 1963) (Eq. 2).
where y was the total biomass per unit area, and N was stand density.
where the coefficient α tends to 3/2 for all species. A theoretical equation (Eq. 3) which was derived by unifying the density effect and the 3/2 power law of self-thinning was proposed (Hagihara 1996): where Y was the stand yield, N 0 the initial stand density, α and β were species specific coefficients. The main result of this formulation was that the density dependent mortality induced for any population an exponential decrease of its constituents after a sufficient lapse of biological time which is the integral of intrinsic growth rate. This result will developed in Results of self-thinned model paragraph.

Density, yield and biomass
We propose a forest dynamic model for managed and non-managed forest which is summarized in the Fig. 1.
This diagram presents the process of modeling biomass. Depending on whether there is thinning or not, two development branches are chosen: No thinning, the stand density noted N sth is calculated and then the yield noted Y sth based on it is also calculated. If there is thinning, the stand density noted N th is calculated as well as the yield noted Y th . From the stand yield Y derived from Y sth or Y th the biomass of the trunk W S , branches W B and leaves W L can be derived. The equations of all models are gathered in this figure too.

Density models
Without forest management, stand density time series followed self-thinning law taking account of the natural competition between trees (Ogawa 2007). The following empirical relation which is a logistic function was therefore used to consider these phenomena (Ogawa 2017) (Eq. 4): where N 0 was the initial stand density (t = 0) and the coefficients ω sth (w.d.), δ sth (w.d.) and γ sth (y −1 ) were fitted by Bayesian inference.
In the other case of managed forest, thinning out was done with different practical scenario, the forecast of stand density N (tree·ha −1 ) could be modelled by a Mitscherlich law. A German soil scientist, E. A. Mitscherlich, developed an equation that relates growth to the supply of growth factors. He observed that the growth response was proportional to the limiting element. Mitscherlich's law states that "the increase in any crop produced by a unit increase in a deficient factor is proportional to the decrease in that factor from the maximum. " The response is curvilinear, not linear (Eq. 5): The maximum possible density obtained by an unlimited supply of time is the product of the initial density N 0 (tree·ha −1 ) and the parameter ω th (w.d.). γ th (y −1 ) is a constant of proportionality. These parameters were fitted by Bayesian inference.

Yield models
Without forest management, stand yield (m 3 · ha −1 ) time series were dependent on stand density by a natural selfthinning relationship developed at first for Japanese cedar then for Japanese cypress (Ogawa 2005), which is a generalization of Eq. 3 proposed by Hagihara (Eq. 6): where Y max (m 3 · ha −1 ) was a maximum yield, α (w.d.), β (w.d.) and n (w.d.) were characteristic coefficients of the self-thinning curve. The formulation was initially calibrated with a series of observations of a natural Sugi plantation in the Nagoya University Experimental Forest at Inabu, located about 55 km east of Nagoya, Aichi Prefecture, central Japan, over a period of 20 years. Parameters values were in this article fitted by Bayesian inference. This Eq. 6 associated to the Eq. 4 concerns the self-thinned model. In the case of managed forest, thinning out was done with different practical scenario, the forecast of forest growth could be modelled by deriving the pipe model from the modelling of DBH, H and N. Actually, Inoue (2006) proposed a model describing the relationship between the form-factors for stem volume and those for stem surface area in coniferous species, especially for Japanese cedar and Japanese cypress. It based on the Kunze's equation which led to the following equation of volume, as it is developped in supplementary data (Eq. 7): We can find this result again by considering v as a cone volume which equation is π 3 (d/200) 2 ×h with the diameter d proportional to the height h.
If the height of a tree is modeled by a Mitscherlich equation with three degrees of freedom, then we can model the yield at the stand level Y (m 3 · ha −1 )(Eq. 8) as follows by having combined the Mitscherlich equation with the equation 7 by multiplying the density N (Eq. 8): where Y max th (m 3 ), ω yth (w.d.) and γ yth (y −1 ) were fitted by Bayesian inference. Some author advise Mitscherlich equation for modelling stand yield (Fukuda et al. 2003) (Eq. 9): where Y max thm (m 3 ), ω yth m (w.d.) and γ yth m (y −1 ) were specific coefficients. The two models named thinned model (Eq. 8 associated to Eq. 5) and Mitscherlich model (Eq. 9 associated to Eq. 5) were compared to determinate the best one with regards the data.

Biomass
Stand stem biomass was deduced from stand yield by wood density ρ wood (dry kg per m 3 ) (Eq. 10): Values of basic wood density were fixed to their average values (Fujiwara et al. 2007) (314 kg·m −3 for Sugi and 407 kg·m −3 for Hinoki). Allometric relationships established by Nishizono et al. were used to model stand branch biomass (Eq. 11) and stand leaf biomass (Nishizono et al. 2005) (Eq. 12):

Data
From the compiled database on Cryptomeria japonica (Sugi) and Chamaecyparis obtusa (Hinoki) (Usoltsev 2013), we extended it with several others publications (In supplementary data 1: Table S1). In total 311 data were collected. Among these data, 234 were coming from Sugi specy and 142 from Hinoki specy. 40 of Hinoki were selfmanaged and 25 of Sugi were self-managed. In managed forest the avaliable data for density were 156 for Sugi and 57 for Hinoki, for stand volume they were 209 for Sugi and 102 for Hinoki, for stem biomass they were 102 for Sugi and 58 for Hinoki, for branch biomass they were 99 for Sugi and 54 for Hinoki and for leaf biomass they were 105 for Sugi and 54 for Hinoki. In managed forest the available data for density were 25 for Sugi and 40 for Hinoki, for stand volume they were 19 for Sugi and 37 for Hinoki, for stem, branch and leaf biomass they were 6 for Hinoki, and no available data for Sugi. These data all concern Japaneses forests except for few data of South Korea. The distribution of Sugi stands in the Japan database is fairly well represented except for some areas north of Sendai such as Akita and Iwate, the Koti Peninsula and the northwestern coastal margins as shown on the map showing the stands locations according to their density (top plot of Fig. S1 in supplementary data 1). On the other hand, Hinoki data are more disparate and cover a small part of Japan (bottom plot of Fig. S1 in supplementary data 1).
Managed forests have been identified and separated from self-thinned forests. Main statistics are gathered in Table 1. These data were used to perform Bayesian inference with the software SAMCAT (Nicoulaud-Gouin et al. 2016), giving a density probability of all model parameters and goodness-of-fit as explained above.

Bayesian inference and model validation
Bayesian inference was conducted from data gathered on Hinoki and Sugi forest without distinguishing species, in this manner: A statistical model f (t; θ) is defined for the observations D(t), and the parameters θ (Eq. 13): where is the random error unexplained by the model f. θ is assumed to be random. It has a prior probability distribution P(θ), reflecting the existing knowledge before the calibration experiment, that can be updated by the observed data D through application of the Bayes' rule (Jeffreys 1961) (Eq. 14): Likelihood was modelled with a normal distribution law centered on the data with the following error model as variance. For the error modelling a proportional error has been chosen for stand density, yield, branch, and leaf biomass, as the data were well fitted with a lognormal distribution. The parameters to be identified concern both the dynamic biomass model θ and the error model (σ ). The components of θ were assigned log-uniform distributions that cover several orders of magnitude and reflect the lack of knowledge of dynamic biomass model parameters (Tables 2, 3). By testing the calibration with several variation domains of the parameters we tested the sensitivity of prior values. As to conclude, often the model gave posterior distribution of some parameters bounded by the low bound or the high bound of the interval of variation. That led to enlarge some intervals in order to obtain good posterior distribution of parameters.
Since the range of density data for managed forests is from 204 trees per hectare to 29500 trees per hectare, we chose an initial density N 0 range that allows us to reach these densities, i.e. [1000,30000]. The ω th parameter being the percentage of the initial density obtained after an infinite time, was chosen so that the product with the initial density is the maximum possible stand density after an infinite time. For the γ th proportionality factor giving the decreasing velocity of density towards the final density we decided to enclose the values found by Fukuda et al. (2003) for this factor around 0.1 for the biomasses.
Since the range of density data for non managed forests is from 406 trees per hectare to 9800 trees per hectare, we chose an initial density N 0 range that allows us to reach these densities, i.e. [2000,20000] and the ω sth parameter so that the product with the initial density is the maximum possible density of the stand after an infinite time. As the proportionality factor is a time decreasing function of γ sth and δ sth , we chose these previous parameter in according to their values found in Table 1 of (Ogawa 2018). For the coefficients of the yield Y according to the Mitscherlich formulation (Eq. 9), we established the bounds of the domains of variation by referring partially to the Table 3 of regression coefficients for the Mitscherlich formula relating the accumulation of wood volume to stand age for Sugi and Hinoki forest types published in Fukuda et al. (2003): The asymptotic value Y max showed values between 575 and 6188, hence the search interval of [500, 6000] chosen. For the parameter ω yth m the values found in Fukuda seem strange to us because they are all greater than 1, but we can express omega as a function of the initial value of efficiency and the asymptotic value (Eq. 15).
If we consider the initial value as zero, we should find 1 as a value, hence our search domain of [0.5, 1] chosen. For gamma we have widened the lower bound of the search interval because successive calibration tests indicated a poor identifiability of the parameter for domains that were too small.
The same reasoning has led to the variation domain of the coefficients for yield Y according to the Thinned formulation (Eq. 8).
For the coefficients of the yield Y of self-managed stands (Eq. 6), we established the ranges of variation from the values found in (Ogawa 2005), α Y = 3/2, n Y = 5.01, β Y = 1.2 , and Y max = (1/141.10 −6 ) 0.833 × 6143 0.5 = 20. Table 1 in Nishizono's publication (Nishizono et al. 2005) shows the values of the regression constants for the branch and leaf biomass models. We chose expanded domains around these average values.
Since the range of leaf biomass values is [5; 27], a range for the W L max coefficient of [10; 30] seems reasonable.
An inverse gamma distribution was chosen for the unknown standard deviation σ of the Gaussian error.
Two performance criteria quantifying the deviation of model median predictionsŷ 1 , · · · ,ŷ n from observations y 1 , · · · , y n were also used. Geometric Reliability Index (GRI) is a dissimilarity measure, quantifying the accuracy factor of predictions around observed values. It is defined by Jachner et al. (2007) (Eq. 16): and optimal values correspond to 1. The efficiency factor EF (Nash and Sutcliffe 1970) is a validation measure, ranging from −∞ to 1 (Eq. 17): where prediction power increases as EF tends to ideal value 1, corresponding to exact model predictions.
Marginal likelihood P(Y |M) for each model M was used to select the most plausible model. This decision criterion is the probability of the observations Y according to model M (Eq. 18): By analogy with classical likelihood ratio tests, a model M has a strong (resp. very strong) evidence if its marginal log-likelihood difference is larger than 3 (resp. 5) in comparison with alternative models (Kass and Raftery 1995). Marginal likelihoods were calculated numerically with SAMCAT by the Laplace-Metropolis approximation (Kass and Raftery 1995).

Results
After analysing the database, results of Mitscherlich and thinned models are given. Then we focus on the model for self-thinning forests (self-thinned model), looking at the 3/2 power law. The yield rate and best fit as well as the correlations between parameters are also analyzed.

Database analyses
An overview of the different variables volume V (in m 3 · ha −1 ), density N (in tree·ha −1 ) in log scale, age and DBH (in cm) shows ( Fig. S2 in supplementary data 1) a correspondence between high DBH, high yields, an age above 50 years and a fairly high stocking density. The high stand densities concern the early development stage below 10 years of age. Some DBH data are not available as shown by the gray circles in the graph. 76 on 160 Hinoki data were all filled for the main variables (Age, DBH, Volume, Density and Height) and 83 on 248 Sugi data were also all filled. 40 data were self-thinned Hinoki stand and 25 data were self-thinned Sugi stand. A non-hierarchical cluster partitioning by K-means method was conducted on the data set of the two species Sugi and Hinoki reduced to a pool where all important variables (Age, DBH, Volume, Density and Height) are filled in to assess if they could be grouped together in order to have unique parameterization in models of yield versus stand density. The Fig.  S3 in supplementary data 1 shows the normalized variable V V max versus the normalized variable N N max for the two species noticed C o for Hinoki and C j for Sugi. Five groups were identified taking account the different variables (Age, DBH, Volume, Density and Height). The fifth group reduced to a Sugi stand is clearly an outlier because of the high stand density it has relative to all the other stands. The second group also stands out with a high density. Then we have three other groups in which the two species are distributed and which show similarities in behavior between species. These two species have been gathered together in order to have an unique parametrization for the models. This implies that our models are independent of species.

Results of mitscherlich and thinned models
On the same data set Mitscherlich model has larger proportional error variances than our model for yield and stand density ( The major interest of our model has been to obtain fine results for the Y max th parameter with a median of 1.91 and [1.36, 2.78] credible interval, compared to the Y max thm parameter with a median of 5370 and [4220, 6000] credible interval. It is due to the fact that the two parameters do not correspond to the same physical expression. The effective yields when age is infinity could be expressed as: Because of the senescence of the forest, it is more realistic to have in median Y t→∞ th = 1897 m 3 ·ha −1 compared to the median value for Mitscherlich model Y t→∞ thm = 5370 m 3 ·ha −1 .
The branches and leaves biomass models are equivalent with respect to the parameters values and errors variance.
The posterior probability distributions of Y max th and γ yth are quite well established showing good identifiability, with regards to the posterior probability distributions of Y max thm , ω yth m and γ yth m .
This result could be forecasted due to biological phenomenon which it is explained based on geometrical consideration: the height is a one-dimensional model and the biomass is a three-dimensional model. Therefore, as height modelling is also a Mitscherlich model the application of this type of model cannot work suitably for yield dynamic growth and the proposal model named thinned model based on this three-dimensional aspect reflects more the reality of yield growth (Ogawa 2019).

Results of self-thinned model
Without forest management, the Ogawa model (Ogawa 2005) was chosen for the stand yield (m 3 ·ha −1 ). Although the self-thinned model was initially calibrated on a single Sugi stand, the pooling of Sugi and Hinoki data to a total of 65 data improves the robustness of the parameter values. With the Bayesian calibration results, we will analyze the 3/2 power law of self-thinning.
Changes in self-thinning slope (Xue et al. 1997) with increasing stand age, showed that the self-thinning slope approaches the self-thinning line of 3/2 which is near the median value of α = 1.4 we found with a credible interval of [0.9, 1.7]. Therefore, results on this database respects the 3/2 power law of self-thinning. The credible interval Table 4 Optimal parameters values of Mitscherlich model (Eqs. 5,9,10,11,12) and thinned model (Eqs. 5,8,10,11,12)  The median and lower and upper credible interval (Highest Posterior Density interval) were gathered takes into account the phenomena of aggregation and clumping which affect the competition process (Weiner 1985). Indeed, spatial pattern can influence self-thinning. With a same stand density, the proximity of neighbors in a clumped stand increases competition between plants compared to regularly spaced stand. The first part of the self-thinned yield model indicates that trajectory gradually approaches and eventually moves along the following self-thinning line at the final stage of stand development (Ogawa 2005). But the n and β parameters allow deviating from this line between the beginning of stand growth and before the end of the stand development. Self-thinning process in even-aged stands consists of creating gaps due to mortality and filling some of these by surviving trees (Zeide 1987).

Yield velocity among models
The Fig. S4 in supplementary data 1 reflects the velocity of yield with regards to maximum of yield for different modelling with parameters median values (Tables 5, 6): for self-thinned and thinned models, the velocity is high at the beginning of tree life to achieve 80% of their maximum at about 80 years or 150 years. Unlike to Mitscherlich model or the model consisting in a proportion of height cube Y Y max ∝ H 3 or the model consisting of a proportion of DBH square times height Y Y max ∝ DBH 2 × H , in which

Goodness-of-fit
Managed observed datas of Sugi were quite well explained by the calibrated models (as well Mitscherlich or thinned models) with a more dispersive fitting with Mitscherlich model (Figs. 2 and 3). Managed observed datas of Hinoki were largely more dispersive far from the plume prediction. There is little significant difference between the Mitscherlich model and the thinned model for both the GRI and EF indices for modeling population density or yield (Table 6). Nevertheless the marginal loglikelihood discriminates between these two models (with a difference of 11.2 on marginal loglikelihood, thinned model has a strong evidence in comparison to Mitscherlich model). This finding confirms the ability of Bayesian inference to identify the correct models of population density, yield and error. One can also see a low goodness-of-fit since for a very good quality of prediction one expects a GRI and an EF close to 1. This low quality of prediction is due to the fact that the data are very largely variable with a great range of values. The self-thinned model is better than the Mitscherlich or thinned models with regards to this goodness-of-fit GRI and EF.

Correlation between parameters
The Mitscherlich model has no correlation between its stand density and yield parameters (Fig. 4). This is because the models of yield and stand density are completely independent of each other. They depend only on the stand age. Two strong correlations are observed, one positive (0.94) between N 0 and γ N on the one hand and negative (-0.97) between Y max and γ Y on the other hand. A lesser correlation (0.57) between ω N and γ N is also observed. No correlation is observed between ω Y and γ Y . Fig. 2 Plots of Mitscherlich dynamic growth model with stand density, stand yield, stand branch biomass and stand leaf biomass curves for the two species of coniferous (Japanese cedar or Sugi, and Japanese cypress or Hinoki) and associated data (from Table  S1 in supplementary data 1). There are the predictions with median parameters, with mean parameters and the plume prediction with 1000 samples of parameters within the posterior distribution coming from Bayesian inference (grid of prediction) Fig. 3 Plots of thinned dynamic growth model with stand density, stand yield, stand branch biomass and stand leaf biomass curves for the two species of coniferous (Japanese cedar or Sugi, and Japanese cypress or Hinoki) and associated data (from Table S1 in supplementary datarefMOESM1). There are the predictions with median parameters, with mean parameters and the plume prediction with 1000 samples of parameters within the posterior distribution coming from Bayesian inference (grid of prediction) Fig. 4 Plots of self-thinned dynamic growth model with stand density, stand yield for the two species of coniferous (Japanese cedar or Sugi, and Japanese cypress or Hinoki) and associated data (from Table S1 in supplementary data 1). There are the predictions with median parameters, with mean parameters and the plume prediction with 1000 samples of parameters within the posterior distribution coming from Bayesian inference (grid of prediction) In contrast to the Mitscherlich model the thinned model exhibits non-negligible correlations between the parameters of density and yield models. For example between ω N and ω Y with a negative correlation (-0.47); between Y max and γ N with a negative correlation (-0.58); between γ N and γ Y with a positive correlation (0.61). We observe the same types of strong intra-model correlations between N 0 and γ N (0.80) and between Y max and γ Y (-0.85). Therefore, the higher the initial density, the higher the exponential decay rate of density, and the higher the Y max , the lower the exponential decay rates γ Y and γ N too. On the other hand, there is a correlation between ω Y and γ Y (0.44) whereas none is observed with the Mitscherlich model. A double and inverse correlation between ω N and N 0 compared to Mitscherlich model, a weaker correlation between ω N and γ N than for Mitscherlich model.
For the self-thinned model, a strong negative correlation between Y max and the self-thinning parameter α is detected (-0.86). The quite strong positive correlation between ω N and the self-thinning parameter α is observed too (0.6). The initial density has weak correlation with all other parameters of stand density or yield. As we saw in the precedent paragraph yield is independent of initial density. The negative quite weak correlation between n Y and intial density N 0 (-0.30) is due to the negative correlation between n Y and δ N (-0.43) and the positive correlation between δ N and N 0 (0.27). We observe strong intra-model correlations for yield model and density model.

Discussions
In this paragraph, we explain why we did not use biomass expansion factors, we discuss the two models adapted to managed forests (the Mitscherlich model and the thinned model), looking carefully at the identifiability of the parameters and the posterior distributions. Concerning the self-thinned model we discuss relative growth and mortality rates, physical and biological times. Then, characteristics of DBH and height models are discussed. Finally, comparison with other models is looked at.

Biomass expansion factors
Biomass expansion factors (Lim et al. 2013)(Eq. 19 and Eq. 20) were not applied. Their formulation indicates that the leaf biomass ratio ( W L W ) on an average tree is decreasing over time, because the biomass expansion factors vary on age class with decreasing values. By Nishizono formulation (Eq. 12) the leaf biomass is limited by the parameter W L max and this whatever the density of the stand. This seems to be in a more realistic physical sense than allocation fractions or BEFs that do not take into account the maximum leaf biomass capacity that a stand of a given species can produce on the basal area. From canopy closure, leaf biomass growth is limited, although log biomass may continue to increase more or less strongly depending on initial stand density. In other words, the fact that leaf biomass per unit area remains constant in closed forests of the same species results from the limited intensity of light, although other factors necessary for photosynthesis such as water, temperature ... may differ from one stand to another. We can say that light is a limiting factor to leaf growth (Tadaki and Kawasaki 1966). If biomass allocation or expansion fraction formulations are kept, the calculation of biomass and leaf area index will be largely distorted and will not reflect the reality for stands of high initial density. Indeed, when the stand biomass becomes maximal or time is over the canopy closure time, foliar biomass approaches an upper limit and the stand leaf biomass can be expressed as (Eq. 21): Therefore, the stand leaf biomass Y F (t) decreases exponentially as stand age increases. A simple allometric relationship between mean leaf mass and mean tree mass does not reflect the age-related changes in forest stand leaf biomass (Ogawa 2018).

Comparison of mitscherlich and thinned models Advantage of bayesian analysis
The advantage of Bayesian analysis over frequentist analysis is the possibility of comparing models without necessarily being nested, thanks to the application of marginal log-likelihood. Bayesian methods need much smaller sample size than frequentist methods to estimate parameters accurately (Zapata-Cuartas et al. 2012). The hierarchical Bayesian framework provides robust approach for estimating models parameters (Price et al. 2009). Moreover, Bayesian analysis directly provides probability densities for the parameters of the models that take into account the data. In addition, the parametric uncertainty obtained allows to know if a parameter is well identifiable.

Identifiability
Raw data used to calibrate Mitscherlich formulas (Eq. 9) in the initial publication (Fukuda et al. 2003) were issue from national stand density control diagram (forestry agency of Japan 1981a; forestry agency of Japan 1981b), and were separately pooled into 10 Sugi and 5 Hinoki regional datasets. However, the variations in R 2 values of Mitscherlich formulas for the accumulation of wood volume per hectare over time (between 0.476 and 0.706 for the two species) were showed without scrutinizing the identifiability of the parameters of the formulas (Table 3 of (Fukuda et al. 2003)). Actually, we have tried to fit this type of equation with regards to our data and it was obvious that the parameters were not well identifiable compared with the suggested model in this paper (Eq. 8 ) for the thinned model (Eq. 9 in supplementary data 1).
The plot of these sensitivity functions, locally calculated at a likely parameter values with the managed forests data observations showed the lack of parameters identifiability for the Mitscherlich model compared to thinned model ( Fig. S5 in supplementary data 1) (Nicoulaud-Gouin et al. 2016). The correlation matrix of the scaled sensitivity matrix had conditioning indices η 1 = 1, η 2 = 10.89 and η 3 = 151.35 for Mitscherlich model and η 1 = 1, η 2 = 3.44 and η 3 = 25.42 for thinned model which indicate possible problems of parameter unstability, parameter redundancy and poor precision of estimation (Seber and Wild 1989) (Table 7).
Parameter identifiability was critical for all parameters Y max thm , ω yth m , γ yth m in Mitscherlich model and for Y max th and γ yth in thinned model (Table 7). Parameter uncertainty was maximum in the direction of the eigenvector associated to conditioning index η 3 (corresponding to the smallest eigenvalue), for which the contributions of Y max thm , ω yth m and γ yth m were high (0.97, 1 and 0.99) for Mitscherlich model and for which the contributions of Y max th and γ yth were high (0.99 and 0.99) for thinned model.

Relative growth and mortality rate
We can express the quotient of relative growth rate on relative mortality rate, and this quotient depends on time and particularly on N N 0 (Eq. 22). Therefore, this quotient is independent on initial stand density N 0 (see equations 10, Table 7 Matrix of the variance decomposition and conditioning indices for Mitscherlich model and thinned model: Collinear parameters were identified in the variance decomposition matrix deduced from the correlation matrix of X (Belsley et al. 1980), whose lines correspond to conditioning indices, and columns correspond to parameters

Conditioning Variance proportion
Mitscherlich model where RGR is the relative growth rate with regard physical time and m is the relative mortality rate with regard physical time.
At the beginning of life stage (about five years old) the relative growth rate is 60 times the relative mortality rate. Stand yield growth is greater than trees mortality. At a infinity time the relative mortality rate is 2.5 times more important than relative growth rate if we consider median values for model parameters (Eq. 22). Indeed, the limit is quickly reached at around 40-50 years old. This can also be seen in the figures of the density and yield curves for the two species Sugi and Hinoki (Fig. 5) which show an inflection around 40 years old. Unfortunately we have not data between 40 years and 140 years to confirm this theory.

Physical and biological time
We can also express physical time in function of stand density (Eq. 23): For this model and the median parameters values of Bayesian inference according to database (Table 5), the valid maximal time is about 180 years with a limit of stand density of 716 tree·ha −1 .
The growth rate of mean tree biomass decreases exponentially with the age and this phenomenon defines a biological time, a specific time of trees species, different from physical time and such that time is slowing down in relation to physical time as the tree evolves as it ages.
We can deduce the biological time where λ(t) is relative to an intrinsic growth rate in the generalized logistic curve (Eq. 16 in supplementary data 1). Therefore, we have a sigmoid curve for the biological time.
The limit of the biological time when the physical time tends towards infinity is lim Extreme Surface Estimator (ESE) and Extremum Distance Estimator (EDE) methods were used to find the inflection point of the convex/concave sigmoid curve (Fig. S6 in supplementary data 1). We find an inflection point between 32 and 43 years.
The stand density can be expressed with biological time as (Eq. 20 in supplementary data 1).
The derivative of stand density with regards physical time, gives us velocity of the stand density decreasing (Eq. 20 in supplementary data 1).
Up to 21 years of rapid loss of stand density can be observed and then fall asleep in old stands. This biological time allows linking self-thinning parameter α with stand density dynamic. The stand yield becomes independent of initial stand density (Eq. 22 in supplementary data 1): In parallel with stand density, the rate of decreasing of stand yield is very high in young stands and decreases until the age of 40 to increase more reasonably until it stabilizes in old stands (Eq. 23 in supplementary data 1).

Comparison with others models
Some models finely characterize the generic parameters of the Richards or Gompertz type models in terms of quantities such as dominant height, basal area, etc., by expressing them as a function of different indexes (site index, last thinning index, time elapsed since the last thinning, etc.) (Candy 1989). These details certainly make it possible to better qualify the models of the managed forests and to reduce the parametric uncertainty, but they require the knowledge of specific quantities for each stand, which are not necessarily accessible, and they also result in a considerable number of modeling parameters.
We did not distinguish between stands with high thinning and those with low thinning, and the models (thinned or Mitscherlich) do not take into account explicitly this thinning weight nor the type of thinning carried out. This weight can be appreciated by an index of relative space between trees by the relation (Eq. 24) and the type of thinning by the change in diameter distribution (Eq. 25).

RS
where N thn , and N tot are respectively removed and total number of stems and G thn and G tot are respectively removed and total basal area.
In fact the thinned model could be expressed with the relative space index, if we model the yield with a product of dominant height model and DBH model instead using cube of Mitscherlich law (Eq. 8).
However, in the study conducted, information on thinning practices is not always available, so it would have been difficult to assign a weight and type of thinning to the different stands in the database to establish a clear classification that would improve the predictive quality of the models.
Our models don't reflect spatial structure of biomass. We did not test general linear models (GLM) of every possible first-order combination of environmental variables like latitue, longitude, elevation, rainfall, ... in order to assess what extent environmental effects can explain biomass spatial variation as it is made in (Guitet et al. 2015).

Conclusions
This paper addressed stand biomass dynamic models of evergreen forests in order to improve biomass growth dynamic assessment at regional scale relying on easily observable parameters (e.g. stand density). Two types of population were considered, self-thinned and managed one, and two species (Japanese cedar and Japanese cypress) were looked at. The models were performed on a large database with Bayesian inference.
For managed forests two types of modelling approaches were compared Mitscherlich one and a new approach called thinned model based on geometric considerations. The results showed an obvious lack of identifiability of Mitscherlich model parameters by the analysis of sensitivity functions and conditioning indices and with the posterior density of input parameters too. By Bayesian inference, marginal loglikelihood showed that thinned model has a strong evidence in comparison to Mitscherlich model too.
The study of the relationships between the parameters of the presented models allowed to highlight the characteristics of the biomass dynamics at different stages of the forest evolution: These models, unlike biomass expansion factor models, provided an ecological view by estimating yield rate dynamics across life stages, growth rates and mortality rates.
Results of self-thinned model brought evidence of light self-thinning law with the 3/2 rule with a α parameter (α = 1.4) quite close to 3/2. Quotient of relative growth rate on relative mortality rate were given as function of physical time and relative mortality rate is 2.5 times more important than relative growth rate after about 40 years old. The growth rate of mean tree biomass can be expressed as a generalized sigmoid curve relative to biological time. Stand density and stand yield can be expressed as function of biological time, showing that yield is independant of initial density.
Managed forest models (thinning or Mitscherlich one) do not allow the expression of the relative growth rate quotient by the relative mortality rate, or the growth rate of an average tree. This is because these models do not incorporate a structural dependence between stand density and stand yield. It is therefore not easy to express biological time as a function of physical time, unlike the self-thinning forest model, which has a structural dependence between stand density and stand yield.
The implementation of these dynamic models with few parameters and easy-to-measure explanatory variables such as stand density will allow them to be integrated into larger forest ecosystem modeling systems for predictive purposes.
The interest of having at one's disposal a parameterization of dynamic biomass models as a function of accessible physical quantities such as stand density, goes beyond the problem of integrating these models into predictive systems within forest ecosystems. These models can be used to dynamically estimate the carbon balance and contribute to a better understanding of climate change factors (Schulze et al. 2021). Dynamic models can also quantify the carbon stored in Sugi and Hinoki forests in specific years when no data are available (Aguirre et al. 2021).