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

Valérie Nicoulaud-Gouin, Marc Andre Gonze, Pierre Hurtevent, Philippe Calmon

<figure><img src="image_1.png" /></figure>

## To cite this version:

Valérie Nicoulaud-Gouin, Marc Andre Gonze, Pierre Hurtevent, Philippe Calmon. Bayesian inference of biomass growth characteristics for sugi (C. japonica) and hinoki (C. obtusa) forests in self-thinned and managed stands. Forest Ecosystems, 2021, 8 (75), pp.1-18. 10.1186/s40663-021-00354-4. hal-03993902

HAL Id: hal-03993902

https://hal.science/hal-03993902v1

Submitted on 17 Feb 2023

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
RESEARCH

Open Access

<figure><img src="image_3.png" /></figure>

Check for
updates

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

Valérie Nicoulaud-Gouin* , Marc-André Gonze, Pierre Hurtevent and Phillippe Calmon

**Abstract**

**Background:** 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.

**Methods:** 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.

**Results:** 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.

**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. 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.

**Keywords:** Dynamic, Biomass, Self-thinning, Sugi, Hinoki, Competition-Density rule

<figure><img src="image_4.png" /></figure>

<footnote>*Correspondence: valerie.nicoulaud_gouin@irsn.fr
Institut de Radioprotection et de Sûreté Nucléaire, St Paul-Lez-Durance, France</footnote>
## 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 re-emitting 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 models (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 self-thinning 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 non-managed 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).

$$ \frac{1}{y} = A + \frac{B}{N} \qquad (1) $$

where y was the total biomass per unit area, and N was stand density.

$$ y = K \times N^{1-\alpha} \qquad (2) $$

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):

$$ Y = K \times N^{1-\alpha} \times \left(1 - \frac{N}{N_0}\right)^{\beta} \qquad (3) $$

where Y was the stand yield, N₀ 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 Nsth is calculated and then the yield noted Ysth based on it is also calculated. If there is thinning, the stand density noted Nth is calculated as well as the yield noted Yth. From the stand yield Y derived from Ysth or Yth the biomass of the trunk WS, branches WB and leaves WL 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):

$$ N = N_0 \left( \omega_{sth} + \frac{(1 - \omega_{sth})(1 + \delta_{sth})}{1 + \delta_{sth} \times e^{\gamma_{sth} \times t}} \right) \qquad (4) $$

where N₀ was the initial stand density (t = 0) and the coefficients ωsth (w.d.), δsth (w.d.) and γsth (y⁻¹) 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⁻¹) 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):

$$ N = N_0 (\omega_{th} + (1 - \omega_{th}) e^{-\gamma_{th} \times t}) \qquad (5) $$

The maximum possible density obtained by an unlimited supply of time is the product of the initial density N₀ (tree·ha⁻¹) and the parameter ωth (w.d.). γth (y⁻¹) is a constant of proportionality. These parameters were fitted by Bayesian inference.

### Yield models

Without forest management, stand yield (m³·ha⁻¹) time series were dependent on stand density by a natural self-thinning 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):

$$ Y = Y_{max} \left( \frac{N}{N_0} \right)^{1-\alpha} \left( 1 - \left( \frac{N}{N_0} \right)^n \right)^{\beta} \qquad (6) $$
where $Y_{max_{thm}}$ (m³), $\omega_{ythm}$ (w.d.) and $\gamma_{ythm}$ (y⁻¹) 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 determine the best one with regards the data.

## Biomass

Stand stem biomass was deduced from stand yield by wood density $\rho_{wood}$ (dry kg per m³) (Eq. 10):

$$W_T = \rho_{wood} \times Y \qquad (10)$$

Values of basic wood density were fixed to their average values (Fujiwara et al. 2007) (314 kg·m⁻³ for Sugi and 407 kg·m⁻³ 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):

$$W_B = a_B \times W_T^{b_B} \qquad (11)$$

$$W_L = \frac{1}{\frac{1}{a_L \times W_T^{b_L}} + \frac{1}{W_{Lmax}}} \qquad (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 self-managed 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 avaliable 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 Japanese 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 north-western 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; \theta)$ is defined for the observations $D(t)$, and the parameters $\theta$ (Eq. 13):

$$D(t) = f(t; \theta) + \epsilon \qquad (13)$$

where $\epsilon$ is the random error unexplained by the model $f$.

$\theta$ is assumed to be random. It has a prior probability distribution $P(\theta)$, 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):

$$P(\theta|D) \propto P(D|\theta) \cdot P(\theta) \qquad (14)$$

Likelihood was modelled with a normal distribution law centered on the data with the following error model as

<table><caption>**Table 1** Statistic data issued from Table S1 in supplementary data 1</caption><thead><tr><th></th><th></th><th></th><th></th><th></th><th></th><th></th><th></th><th></th></tr></thead><tbody><tr><td colspan="9"><strong>Thinned forests</strong></td></tr><tr><td>Japanese cedar</td><td>Age</td><td>DBH</td><td><i>H</i></td><td><i>N</i></td><td><i>V</i></td><td><i>W</i><sub>B</sub></td><td><i>W</i><sub>L</sub></td><td></td></tr><tr><td>Min</td><td>5</td><td>1</td><td>3</td><td>204</td><td>12</td><td>0.5</td><td>6</td><td></td></tr><tr><td>Mean</td><td>41</td><td>23</td><td>16</td><td>2188</td><td>447</td><td>9</td><td>20</td><td></td></tr><tr><td>Max</td><td>153</td><td>65</td><td>40</td><td>29500</td><td>1726</td><td>30</td><td>35</td><td></td></tr><tr><td>Japanese cypress</td><td>Age</td><td>DBH</td><td><i>H</i></td><td><i>N</i></td><td><i>V</i></td><td><i>W</i><sub>B</sub></td><td><i>W</i><sub>L</sub></td><td></td></tr><tr><td>Min</td><td>17</td><td>8</td><td>5</td><td>350</td><td>74</td><td>7</td><td>5</td><td></td></tr><tr><td>Mean</td><td>38</td><td>18</td><td>13</td><td>1883</td><td>302</td><td>16</td><td>15</td><td></td></tr><tr><td>Max</td><td>91</td><td>34</td><td>22</td><td>4170</td><td>560</td><td>32</td><td>27</td><td></td></tr><tr><td colspan="9"><strong>Unthinned forests</strong></td></tr><tr><td>Japanese cedar</td><td>Age</td><td>DBH</td><td><i>H</i></td><td><i>N</i></td><td><i>V</i></td><td><i>W</i><sub>B</sub></td><td><i>W</i><sub>L</sub></td><td></td></tr><tr><td>Min</td><td>4</td><td>2</td><td>35</td><td>406</td><td>376</td><td>(-)</td><td>(-)</td><td></td></tr><tr><td>Mean</td><td>44</td><td>26</td><td>35.8</td><td>3200</td><td>679</td><td>(-)</td><td>(-)</td><td></td></tr><tr><td>Max</td><td>153</td><td>46</td><td>36</td><td>5020</td><td>1200</td><td>(-)</td><td>(-)</td><td></td></tr><tr><td>Japanese cypress</td><td>Age</td><td>DBH</td><td><i>H</i></td><td><i>N</i></td><td><i>V</i></td><td><i>W</i><sub>B</sub></td><td><i>W</i><sub>L</sub></td><td></td></tr><tr><td>Min</td><td>3</td><td>7</td><td>7</td><td>900</td><td>0</td><td>7</td><td>8</td><td></td></tr><tr><td>Mean</td><td>30</td><td>11</td><td>11</td><td>6100</td><td>199</td><td>13</td><td>14</td><td></td></tr><tr><td>Max</td><td>105</td><td>20</td><td>15</td><td>9800</td><td>393</td><td>19</td><td>22</td><td></td></tr></tbody></table>

Minimum, maximum and mean for the main variables as Age, DBH, *H*, *V*, *W*<sub>B</sub>, *W*<sub>L</sub> are gathered for the two species Japanese cedar and Japanese cypress and distinguishing managed forests (Thinned ones) and unmanaged forests (unthinning ones)
<table><caption>**Table 2** Prior distributions for parameters of Mitscherlich model (Eqs. 5, 9, 10, 11, 12) and thinned model (Eqs. 5, 8, 10, 11, 12)</caption><thead><tr><th colspan="4">Mitscherlich Model</th><th colspan="4">Thinned Model</th></tr><tr><th>Model parameters</th><th>Law</th><th>Law parameter 1</th><th>Law parameter 2</th><th>Model parameters</th><th>Law</th><th>Law parameter 1</th><th>Law parameter 2</th></tr></thead><tbody><tr><td>Y<sub>max<sub>thm</sub></sub></td><td>LogUniform</td><td>500</td><td>6000</td><td>Y<sub>max<sub>th</sub></sub></td><td>LogUniform</td><td>1</td><td>50</td></tr><tr><td>ω<sub>ythm</sub></td><td>LogUniform</td><td>0.5</td><td>1</td><td>ω<sub>yth</sub></td><td>LogUniform</td><td>0.5</td><td>1</td></tr><tr><td>γ<sub>ythm</sub></td><td>LogUniform</td><td>1e-12</td><td>0.01</td><td>γ<sub>yth</sub></td><td>LogUniform</td><td>0.001</td><td>0.1</td></tr><tr><td>σ<sub>ythm</sub></td><td>InvGamma</td><td>3</td><td>0.05</td><td>σ<sub>yth</sub></td><td>InvGamma</td><td>3</td><td>0.05</td></tr><tr><td>N<sub>0</sub></td><td>LogUniform</td><td>1000</td><td>30000</td><td>N<sub>0</sub></td><td>LogUniform</td><td>1000</td><td>30000</td></tr><tr><td>ω<sub>N</sub></td><td>LogUniform</td><td>0.01</td><td>1</td><td>ω<sub>N</sub></td><td>LogUniform</td><td>0.01</td><td>1</td></tr><tr><td>γ<sub>N</sub></td><td>LogUniform</td><td>0.01</td><td>1.5</td><td>γ<sub>N</sub></td><td>LogUniform</td><td>0.01</td><td>1.5</td></tr><tr><td>σ<sub>N</sub></td><td>InvGamma</td><td>3</td><td>0.05</td><td>σ<sub>N</sub></td><td>InvGamma</td><td>3</td><td>0.05</td></tr><tr><td>a<sub>B</sub></td><td>LogUniform</td><td>1e-5</td><td>1</td><td>a<sub>B</sub></td><td>LogUniform</td><td>0.05</td><td>1</td></tr><tr><td>b<sub>B</sub></td><td>LogUniform</td><td>0.1</td><td>2</td><td>b<sub>B</sub></td><td>LogUniform</td><td>0.1</td><td>2</td></tr><tr><td>σ<sub>B</sub></td><td>InvGamma</td><td>3</td><td>0.05</td><td>σ<sub>B</sub></td><td>InvGamma</td><td>3</td><td>0.05</td></tr><tr><td>a<sub>L</sub></td><td>LogUniform</td><td>1e-12</td><td>1</td><td>a<sub>L</sub></td><td>LogUniform</td><td>0.1</td><td>3</td></tr><tr><td>b<sub>L</sub></td><td>LogUniform</td><td>0.1</td><td>10</td><td>b<sub>L</sub></td><td>LogUniform</td><td>0.1</td><td>5</td></tr><tr><td>W<sub>L<sub>max</sub></sub></td><td>LogUniform</td><td>10</td><td>30</td><td>W<sub>L<sub>max</sub></sub></td><td>LogUniform</td><td>10</td><td>30</td></tr><tr><td>σ<sub>WL</sub></td><td>InvGamma</td><td>3</td><td>0.05</td><td>σ<sub>WL</sub></td><td>InvGamma</td><td>3</td><td>0.05</td></tr></tbody></table>

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 $\theta$ and the error model ($\sigma$). The components of $\theta$ 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 $\omega_{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 $\gamma_{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 $\omega_{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 $\gamma_{sth}$ and $\delta_{sth}$, we chose these previous

<table><caption>**Table 3** Prior distributions for parameters of self-thinned model (Eqs. 4, 6, 10, 11, 12)</caption><thead><tr><th>Model parameters</th><th>Law</th><th>Law parameter 1</th><th>Law parameter 2</th></tr></thead><tbody><tr><td>Y<sub>max<sub>thm</sub></sub></td><td>LogUniform</td><td>5</td><td>5000</td></tr><tr><td>α</td><td>LogUniform</td><td>0.8</td><td>2</td></tr><tr><td>n<sub>γ</sub></td><td>LogUniform</td><td>0.001</td><td>8</td></tr><tr><td>β<sub>γ</sub></td><td>LogUniform</td><td>0.2</td><td>6</td></tr><tr><td>σ<sub>γ</sub></td><td>InvGamma</td><td>3</td><td>0.05</td></tr><tr><td>N<sub>0</sub></td><td>LogUniform</td><td>2000</td><td>20000</td></tr><tr><td>ω<sub>N</sub></td><td>LogUniform</td><td>0.001</td><td>0.5</td></tr><tr><td>δ<sub>N</sub></td><td>LogUniform</td><td>0.001</td><td>0.5</td></tr><tr><td>γ<sub>N</sub></td><td>LogUniform</td><td>0.01</td><td>0.3</td></tr><tr><td>σ<sub>N</sub></td><td>InvGamma</td><td>3</td><td>0.05</td></tr><tr><td>a<sub>B</sub></td><td>LogUniform</td><td>0.5</td><td>1</td></tr><tr><td>b<sub>B</sub></td><td>LogUniform</td><td>0.1</td><td>2</td></tr><tr><td>σ<sub>B</sub></td><td>InvGamma</td><td>3</td><td>0.05</td></tr><tr><td>a<sub>L</sub></td><td>LogUniform</td><td>0.01</td><td>3</td></tr><tr><td>b<sub>L</sub></td><td>LogUniform</td><td>0.5</td><td>3</td></tr><tr><td>W<sub>L<sub>max</sub></sub></td><td>LogUniform</td><td>10</td><td>30</td></tr><tr><td>σ<sub>WL</sub></td><td>InvGamma</td><td>3</td><td>0.05</td></tr></tbody></table>
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 $\omega_{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).

$$ \omega_{yth_m} = \frac{Y_{max} - Y_0}{Y_{max}} \qquad (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), $\alpha_Y = 3/2$, $n_Y = 5.01$, $\beta_Y = 1.2$, and $Y_{max} = (1/141.10^{-6})^{0.833} \times 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_{Lmax}$ coefficient of [10; 30] seems reasonable.

An inverse gamma distribution was chosen for the unknown standard deviation $\sigma$ of the Gaussian error.

Two performance criteria quantifying the deviation of model median predictions $\hat{y}_1, \dots, \hat{y}_n$ from observations $y_1, \dots, 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):

$$ GRI = \frac{1 + \sqrt{\frac{1}{n} \sum_i \left( \frac{\hat{y}_i - y_i}{y_i + y_i} \right)^2}}{1 - \sqrt{\frac{1}{n} \sum_i \left( \frac{\hat{y}_i - y_i}{y_i + y_i} \right)^2}} \qquad (16) $$

and optimal values correspond to 1. The efficiency factor EF (Nash and Sutcliffe 1970) is a validation measure, ranging from $-\infty$ to 1 (Eq. 17):

$$ EF = 1 - \frac{\sum_i (\hat{y}_i - y_i)^2}{\sum_t (\hat{y}_i - \bar{y})^2} \qquad (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):

$$ P(Y|M) = \int P(Y|\theta)P(\theta)d\theta \qquad (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³·ha⁻¹), density N (in tree·ha⁻¹) 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 $\frac{V}{V_{max}}$ versus the normalized variable $\frac{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 (Table 4):

* for the yield (Eqs. 8 and 9) there is a $\sigma_{Y_{thm}}$ equal to 41.4 with [33.8, 50.1] as credible interval against a $\sigma_{Y_{th}}$ equal to 36.6 with [29.6, 44.2] as credible interval;

* for the stand density (Eq. 5) there is a $\sigma_{N_{thm}}$ equal to 604 with [508, 717] as credible interval against a $\sigma_{N_{th}}$ equal to 598 with [505, 707] as credible interval.

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: $Y_{t \to \infty_{th}} = Y_{max_{th}} \times N_0 \times \omega_{th}$ and $Y_{t \to \infty_{thm}} = Y_{max_{thm}}$. Because of the senescence of the forest, it is more realistic to have in median $Y_{t \to \infty_{th}} = 1897 \text{ m}^3 \cdot \text{ha}^{-1}$ compared to the median value for Mitscherlich model $Y_{t \to \infty_{thm}} = 5370 \text{ m}^3 \cdot \text{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 $\gamma_{yth}$ are quite well established showing good identifiability, with regards to the posterior probability distributions of $Y_{max_{thm}}$, $\omega_{ythm}$ and $\gamma_{ythm}$.

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 ($\text{m}^3 \cdot \text{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 $\alpha = 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><caption>**Table 4** Optimal parameters values of Mitscherlich model (Eqs. 5, 9, 10, 11, 12) and thinned model (Eqs. 5, 8, 10, 11, 12)</caption><thead><tr><th colspan="3">Mitscherlich Model</th><th colspan="3">Thinned Model</th></tr><tr><th>Model parameters</th><th>Lower CI</th><th>Median</th><th>Upper CI</th><th>Model parameters</th><th>Lower CI</th><th>Median</th><th>Upper CI</th></tr></thead><tbody><tr><td>Y<sub>max<sub>thm</sub></sub></td><td>4220</td><td>5370</td><td>6000</td><td>Y<sub>max<sub>th</sub></sub></td><td>1.36</td><td>1.96</td><td>2.78</td></tr><tr><td>ω<sub>ythm</sub></td><td>0.99</td><td>1</td><td>1</td><td>ω<sub>yth</sub></td><td>0.95</td><td>0.97</td><td>0.99</td></tr><tr><td>γ<sub>ythm</sub></td><td>1.610<sup>-3</sup></td><td>1.910<sup>-3</sup></td><td>2.510<sup>-3</sup></td><td>γ<sub>yth</sub></td><td>1.410<sup>-2</sup></td><td>1.810<sup>-2</sup></td><td>2.110<sup>-2</sup></td></tr><tr><td>σ<sub>ythm</sub></td><td>33.8</td><td>41.4</td><td>50.1</td><td>σ<sub>yth</sub></td><td>29.6</td><td>36.6</td><td>44.2</td></tr><tr><td>N<sub>0</sub></td><td>1590</td><td>2070</td><td>2630</td><td>N<sub>0</sub></td><td>15000</td><td>18600</td><td>22200</td></tr><tr><td>ω<sub>N</sub></td><td>3.810<sup>-2</sup></td><td>4.710<sup>-2</sup></td><td>5.710<sup>-2</sup></td><td>ω<sub>N</sub></td><td>4.310<sup>-2</sup></td><td>5.210<sup>-2</sup></td><td>6.210<sup>-2</sup></td></tr><tr><td>γ<sub>N</sub></td><td>0.09</td><td>0.11</td><td>0.13</td><td>γ<sub>N</sub></td><td>0.08</td><td>0.10</td><td>0.11</td></tr><tr><td>σ<sub>N</sub></td><td>508</td><td>604</td><td>717</td><td>σ<sub>N</sub></td><td>505</td><td>598</td><td>707</td></tr><tr><td>a<sub>B</sub></td><td>0.06</td><td>0.11</td><td>0.19</td><td>a<sub>B</sub></td><td>6.210<sup>-2</sup></td><td>0.14</td><td>0.24</td></tr><tr><td>b<sub>B</sub></td><td>0.89</td><td>1.01</td><td>1.11</td><td>b<sub>B</sub></td><td>0.85</td><td>0.98</td><td>1.13</td></tr><tr><td>σ<sub>B</sub></td><td>0.88</td><td>1.12</td><td>1.41</td><td>σ<sub>B</sub></td><td>0.92</td><td>1.18</td><td>1.4</td></tr><tr><td>a<sub>L</sub></td><td>1.510<sup>-6</sup></td><td>0.03</td><td>0.72</td><td>a<sub>L</sub></td><td>0.1</td><td>0.96</td><td>2.6</td></tr><tr><td>b<sub>L</sub></td><td>2.42</td><td>5.8</td><td>9.5</td><td>b<sub>L</sub></td><td>2.05</td><td>3.39</td><td>4.96</td></tr><tr><td>W<sub>Lmax</sub></td><td>17.6</td><td>18.5</td><td>19.4</td><td>W<sub>Lmax</sub></td><td>17.6</td><td>18.5</td><td>19.4</td></tr><tr><td>σ<sub>WL</sub></td><td>1.42</td><td>1.8</td><td>2.22</td><td>σ<sub>WL</sub></td><td>1.4</td><td>1.8</td><td>2.2</td></tr></tbody></table>

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 ($\left(\frac{Y}{Y_{max}} \propto H^3\right)$ or the model consisting of a proportion of DBH square times height ($\left(\frac{Y}{Y_{max}} \propto DBH^2 \times H\right)$), in which

<table><caption>**Table 5** Optimal parameters values of self-thinned model (Eqs. 4, 6, 10, 11, 12)</caption><thead><tr><th>Model parameters</th><th>Lower CI</th><th>Median</th><th>Upper CI</th></tr></thead><tbody><tr><td>Y<sub>max</sub><sub>thm</sub></td><td>185</td><td>366</td><td>1090</td></tr><tr><td>α</td><td>0.9</td><td>1.4</td><td>1.7</td></tr><tr><td>n<sub>Y</sub></td><td>0.8</td><td>4</td><td>7</td></tr><tr><td>β<sub>Y</sub></td><td>1.8</td><td>3.7</td><td>5.8</td></tr><tr><td>σ<sub>Y</sub></td><td>36</td><td>56</td><td>82</td></tr><tr><td>N<sub>0</sub></td><td>8500</td><td>9940</td><td>11400</td></tr><tr><td>ω<sub>N</sub></td><td>3.510<sup>-2</sup></td><td>7.210<sup>-2</sup></td><td>0.12</td></tr><tr><td>δ<sub>N</sub></td><td>5.510<sup>-2</sup></td><td>0.25</td><td>0.48</td></tr><tr><td>γ<sub>N</sub></td><td>4.710<sup>-2</sup></td><td>6.610<sup>-2</sup></td><td>9.710<sup>-2</sup></td></tr><tr><td>σ<sub>N</sub></td><td>238</td><td>349</td><td>490</td></tr><tr><td>a<sub>B</sub></td><td>0.05</td><td>0.33</td><td>0.88</td></tr><tr><td>b<sub>B</sub></td><td>0.49</td><td>0.71</td><td>1.04</td></tr><tr><td>σ<sub>B</sub></td><td>0.2</td><td>0.6</td><td>1.4</td></tr><tr><td>a<sub>L</sub></td><td>0.01</td><td>0.52</td><td>2.4</td></tr><tr><td>b<sub>L</sub></td><td>0.57</td><td>1.52</td><td>2.7</td></tr><tr><td>W<sub>Lmax</sub></td><td>11.5</td><td>15.5</td><td>25.3</td></tr><tr><td>σ<sub>W</sub><sub>L</sub></td><td>0.26</td><td>0.76</td><td>1.8</td></tr></tbody></table>

The median and lower and upper credible interval (Highest Posterior Density interval) were gathered

<table><caption>**Table 6** With the Bayesian inference, goodness of fit and marginal likelihood of the three models Mitscherlich model, thinned model and self-thinned model</caption><thead><tr><th colspan="4">N</th></tr><tr><th>Model</th><th>GRI</th><th>EF</th><th>ML</th></tr></thead><tbody><tr><td>Mitscherlich model</td><td>2.29</td><td>0.30</td><td>-4940.9</td></tr><tr><td>thinned model</td><td>2.28</td><td>0.31</td><td>-4929.6</td></tr><tr><td>self-thinned model</td><td>1.52</td><td>0.65</td><td>-1006.5</td></tr><tr><td colspan="4">Y</td></tr><tr><td>Mitscherlich model</td><td>4.53</td><td>0.28</td><td>-4940.9</td></tr><tr><td>thinned model</td><td>4.52</td><td>0.36</td><td>-4929.7</td></tr><tr><td>self-thinned model</td><td>2.50</td><td>-0.05</td><td>-1006.8</td></tr></tbody></table>

GRI: geometric reliability index, EF: efficiency factor and ML: marginal loglikelihood

the velocity does not overtake 30% of their maximum at about 200 years.

### 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<sub>0</sub> and γ<sub>N</sub> on the one hand and negative (-0.97) between Y<sub>max</sub> and γ<sub>Y</sub> on the other hand. A lesser correlation (0.57) between ω<sub>N</sub> and γ<sub>N</sub> is also observed. No correlation is observed between ω<sub>Y</sub> and γ<sub>Y</sub>.
<figure><img src="image_5.png" /><figcaption>**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)</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>**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 data ref MOESM1). 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)</figcaption></figure>
<figure><img src="image_7.png" /></figure>

In contrast to the Mitscherlich model the thinned model
exhibits non-negligible correlations between the parame-
ters of density and yield models. For example between ωN
and ωY with a negative correlation (-0.47); between Ymax
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 N0
and γN (0.80) and between Ymax and γY (-0.85). Therefore,
the higher the initial density, the higher the exponen-
tial decay rate of density, and the higher the Ymax, 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 Mitscher-
lich model. A double and inverse correlation between
ωN and N0 compared to Mitscherlich model, a weaker
correlation between ωN and γN than for Mitscherlich
model.

For the self-thinned model, a strong negative corre-
lation between Ymax and the self-thinning parameter α
is detected (-0.86). The quite strong positive correla-
tion between ωN and the self-thinning parameter α is

observed too (0.6). The initial density has weak correla-
tion with all other parameters of stand density or yield.
As we saw in the precedent paragraph yield is indepen-
dent of initial density. The negative quite weak correla-
tion between nY and intial density N0 (-0.30) is due to
the negative correlation between nY and δN (-0.43) and
the positive correlation between δN and N0 (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 mod-
els adapted to managed forests (the Mitscherlich model
and the thinned model), looking carefully at the iden-
tifiability of the parameters and the posterior distribu-
tions. Concerning the self-thinned model we discuss rel-
ative 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

$$W_B = BEF_B \times W_T \qquad (19)$$

$$W_L = BEF_L \times W_T \qquad (20)$$

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_{Lmax}$ 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):

$$Y_F = W_{Lmax} \times N = W_{Lmax} \times N_0 \left( \omega_{sth} + \frac{(1 - \omega_{sth})(1 + \delta_{sth})}{1 + \delta_{sth} \times e^{\gamma_{sth} \times t}} \right) \qquad (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). The stand yield $Y$ function was derived to obtain analytical expressions of ($\frac{\partial Y}{\partial Y_{max_{thm}}}$, $\frac{\partial Y}{\partial \omega_{thm}}$ and $\frac{\partial Y}{\partial \gamma_{thm}}$) for the Mitscherlich model (Eq. 8 in supplementary data 1) and ($\frac{\partial Y}{\partial Y_{max_{th}}}$, $\frac{\partial Y}{\partial \omega_{th}}$ and $\frac{\partial Y}{\partial \gamma_{th}}$) 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 $\eta_1 = 1$, $\eta_2 = 10.89$ and $\eta_3 = 151.35$ for Mitscherlich model and $\eta_1 = 1$, $\eta_2 = 3.44$ and $\eta_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}}$, $\omega_{thm}$, $\gamma_{thm}$ in Mitscherlich model and for $Y_{max_{th}}$ and $\gamma_{th}$ in thinned model (Table 7). Parameter uncertainty was maximum in the direction of the eigenvector associated to conditioning index $\eta_3$ (corresponding to the smallest eigenvalue), for which the contributions of $Y_{max_{thm}}$, $\omega_{thm}$ and $\gamma_{thm}$ were high (0.97, 1 and 0.99) for Mitscherlich model and for which the contributions of $Y_{max_{th}}$ and $\gamma_{th}$ were high (0.99 and 0.99) for thinned model.

### Self-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><caption>**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</caption><thead><tr><th>Conditioning</th><th colspan="3">Variance proportion</th></tr></thead><tbody><tr><td colspan="4">Mitscherlich model</td></tr><tr><td>Indices</td><td>S<sub>γmax<sub>thm</sub></sub></td><td>S<sub>ω<sub>γthm</sub></sub></td><td>S<sub>γ<sub>γthm</sub></sub></td></tr><tr><td>η<sub>1</sub> = 1</td><td>0</td><td>0</td><td>0</td></tr><tr><td>η<sub>2</sub> = 10.89</td><td>0.031</td><td>0</td><td>0.008</td></tr><tr><td>η<sub>3</sub> = 151.35</td><td>0.97</td><td>1</td><td>0.99</td></tr><tr><td colspan="4">Thinned model</td></tr><tr><td>Indices</td><td>S<sub>γmax<sub>th</sub></sub></td><td>S<sub>ω<sub>γth</sub></sub></td><td>S<sub>γ<sub>γth</sub></sub></td></tr><tr><td>η<sub>1</sub> = 1</td><td>0.001</td><td>0.011</td><td>0.001</td></tr><tr><td>η<sub>2</sub> = 3.44</td><td>0.013</td><td>0.296</td><td>0.002</td></tr><tr><td>η<sub>3</sub> = 25.42</td><td>0.99</td><td>0.69</td><td>0.99</td></tr></tbody></table>

For a large conditioning index ηᵢ, parameters are collinear if their contributions to the variance are high

11 and 12 in supplementary data 1).

$$ \lim_{t \to \infty} \frac{RGR}{m} = \lim_{t \to \infty} -\frac{dlog(Y)}{dlog(N)} = (\alpha - 1) + \beta \times n \times \frac{\omega^n}{1 - \omega^n} \approx 0.4 \quad (22) $$

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):

$$ t = \frac{1}{\gamma} \log \left( \frac{1}{\delta} \left( \frac{(1 - \omega)(1 + \delta)N_0}{N - \omega N_0} - 1 \right) \right) \quad (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⁻¹.

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 $τ(t) = \int_0^t \lambda(t)dt$ 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_{t \to \infty} τ = -αlog(ω) \approx 3.68$. The 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
<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /><figcaption>**Fig. 5** Correlation matrix of Mitscherlich, thinned and self-thinned model parameters (respectively at top, middle and bottom)</figcaption></figure>
thinning by the change in diameter distribution (Eq. 25).

$$RS = \frac{\sqrt{10000/N}}{H} \qquad (24)$$

$$SG = \frac{N_{thm}/N_{tot}}{G_{thm}/G_{tot}} \qquad (25)$$

where $N_{thm}$, and $N_{tot}$ are respectively removed and total number of stems and $G_{thm}$ 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 latitude, 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 $\alpha$ parameter ($\alpha$ = 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 independent 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).

## Supplementary Information

The online version contains supplementary material available at https://doi.org/10.1186/s40663-021-00354-4.

Additional file 1: Supplementary data 1

### Acknowledgments

This work has been done thanks to the State financial support managed by the Agence Nationale de la Recherche, allocated in "Investissements d'Avenir" framework programme under reference AMORAD,ANR-11-RNSR-0002.

### Authors' contributions

Conception and design of the study: VNG. Acquisition of data: All the authors. Drafting the manuscript: VNG. Revising the manuscript critically for important intellectual content: VNG, MAG, PH. All the authors read and approved the final manuscript

### Declarations

### Competing interests

The authors declare that they have no competing interests.
Received: 18 June 2021 Accepted: 19 October 2021

Published online: 28 December 2021

## References

Aguirre A, del Río M, Ruiz-Peinado R, Condès S (2021) Stand-level biomass models for predicting C stock for the main Spanish pine species. For Ecosyst 8:29. https://doi.org/10.1186/s40663-021-00308-w

Aoki D, Asai R, Tomioka R, Matsushita T, Asakura H, Tabuchi M, Fukushima K (2017) Translocation of 133Cs administered to Cryptomeria japonica wood. Sci Total Environ 584(585):88-95

Bartelink H. H, Kramer K, Mohren G. M. J (1997) Applicability of the radiation-use efficiency concept for simulating growth of forest stands. Agric For Meteorol 88(1):169-179

Belsley D. A, Kuh E, Welsch R. E (1980) Regression Diagnostics: Identifying Influential Data and Sources of Collinearity, chapitre Detecting and Assessing Collinearity. John Wiley & Sons, Inc.

Bugai D. A, Waters R. D, Dzhepo S. P, Skal'skij A. S (1996) Risks from radionuclide migration to groundwater in the Chernobyl 30-km zone. Health Phys 71(1):9-19

Candy S. G (1989) Growth and yield models for Pinus radiata in Tasmania NZ. J For Sci 19:112-133

Canham C. D, LePage P. T, Dave Coates K (2004) A neighborhood analysis of canopy tree competition: effects of shading versus crowding. Can J For Res 34(4):778-787

Cao Q. V (2014) Linking individual-tree and whole-stand models for forest growth and yield prediction. For Ecosyst 1:1-8

Evangeliou N, Balkanski Y, Cozic A, Hao W. M, Mouillot F, Thonicke K, Paugam R, Zibtsev S, Mousseau T. A, Wang R, Poulter B, Petkov A, Yue C, Cadule P, Koffi B, Kaiser J. W, Moller A. P (2016) Fire evolution in the radioactive forests of Ukraine and Belarus: future risks for the population and the environment. Ecol Monogr 85(1):49-72

Fang J, Kato T, Guo Z, Yang Y, Hu H, Shen H, Zhao X, Kishimoto-Mo AW, Tang Y, Houghton RA (2014) Evidence for environmentally enhanced forest growth. Proc Natl Acad Sci U S A 111(26):9527-32. https://doi.org/10.1073/pnas.1402333111

Fengri L, Lingbin M (1995) Maximum density line and model of self-thinning for plantation. J Northeast For Univ 6:1-7. https://doi.org/10.1007/BF02842944

Forestry A (2012) Annual Report on Forest and Forestry in Japan 2012 Fiscal year 2011. Rapp Tech For Agency

forestry agency of Japan (1981a) Manual of Stand Density Control Diagram for Hinoki Plantations. Japan Forest Technology Association, Tokyo

forestry agency of Japan (1981b) Manual of Stand Density Control Diagram for Sugi Plantations. Japan Forest Technology Association, Tokyo

Fujiwara T, Yamashita K, Kuroda K (2007) Basic densities as a parameter for estimating the amount of carbon removal by forests and their variation. Bull FFPRI 6(4(405)):215-226

Fukuda M, Iehara T, Matsumoto M (2003) Carbon stock estimates for sugi and hinoki forests in Japan. For Ecol Manag 184(1(3)):1-16

Fukuyama T, Takenaka C, Onda Y (2005) 137Cs loss via soil erosion from a mountainous headwater catchment in central Japan. Sci Total Environ 350(1):238-247

Guitet S, Hérault B, Molto Q, Brunaux O, Couteron P (2015) Spatial Structure of Above-Ground Biomass Limits Accuracy of Carbon Mapping in Rainforest but Large Scale Forest Inventories Can Help to Overcome. PLOS ONE 10:1-22. Public Library of Science

Hagihara A (1996) Theoretical studies on the time-trajectory of mean plant weight and density in self-thinning populations. Bull Nagoya Univ For 15:51-67

Inoue A (2006) A model for the relationship between form-factors for stem volume and those for stem surface area in coniferous species. J For Res 11(4):289-294

Jachner S, van den Boogaart K. G., Petzoldt T (2007) Statistical methods for the qualitative assessment of dynamic models with time delay (r package qualv). J Stat Soft 22:1-30

Jeffreys H (1961) Theory of Probability. Oxford Univ. Press, Oxford

Karki K. B, Shibano H (2006) Soil loss in a forested watershed underlain by deeply weathered granite: comparison of observations to predictions of a GIS-based USLE. Bull Tokyo Univ For 115:1-36

Kass R. E, Raftery A. E (1995) Bayes factors. J Am Stat Assoc 430:773-795

Lehtonen A, Mäkipää R, Heikkinen J, Sievänen R, Liski J (2004) Biomass expansion factors (BEFs) for Scots pine, Norway spruce and birch according to stand age for boreal forests. For Ecol Manag 188:211-224

Lim H, Lee K-H, Lee K. H, Park I. H (2013) Biomass expansion factors and allometric equations in an age sequence for Japanese cedar (Cryptomeria japonica) in southern Korea. J For Res 18(4):316-322

Ministry of Agriculture, Forestry and Fisheries of Japan (2012) The 87th Statistical Yearbook of Ministry of Agriculture, Forestry and Fisheries [in japanese]

Matsumoto M (1997) Construction of yield tables for sugi (Cryptomeria japonica). In Kumamoto district using LYCS. J For Plann 3:55-62

Nagahama T, Kondoh H (2006) Improved Yield Prediction Information for Sugi (Cryotomeria japonica) Plantations Applied to Long-rotation Management in Kagoshima Prefecture. J Japan For Soc 88(2):71-78

Nakajima T, Yang YC, Matsumoto M, Zheng CL, Tatsuhara S, Lin WL, Shiraishi N (2014) Estimation of diameter growth parameters for Cryptomeria Plantations in Taiwan Using the Local Yield Table Construction System. NCHUIP:14-21

Nash J. E, Sutcliffe F (1970) River flow forecasting through conceptual models. Part I. A discussion of principles. J. Hydrol 10:282-290

Nicoulaud-Gouin V, Garcia-Sanchez L, Giacalone M, Attard J. C, Martin-Garin A, Bois F. Y (2016) Identifiability of sorption parameters in stirred flow-through reactor experiments and their identification with a Bayesian approach. J Environ Radioact 162(163):328-339

Nishizono T, Iehara T, Kuboyama H, Fukuda M (2005) A forest biomass yield table based on an empirical model. J For Res 10(3):211-220

Ogawa H, Hirano Y, Igei S, Yokota K, Arai S, Ito H, Kumata A, Yoshida H (2016) Changes in the distribution of radiocesium in the wood of Japanese cedar trees from 2011 to 2013. J Environ Radioact 161:51-57. Special Issue: Radioactive contamination in forest ecosystems: From Chernobyl to Fukushima

Ogawa K (2005) Time-trajectory of mean phytomass and density during a course of self-thinning in a sugi (Cryptomeria japonica D. Don) plantation. For Ecol Manag 214(1):104-110

Ogawa K (2007) Self-thinning and Stand Dynamics of even-aged pure Stands of Sugi (Cryptomeria japonica) and Kinoki Cypress (Chamaecyparis obtusa). In: Scaggs AK (ed). New research in forest ecology. Nova Science Publishers, New York

Ogawa K (2017) Modeling age-related leaf biomass changes in forest stands under the assumptions of the self-thinning law. Trees 31:165-172. https://doi.org/10.1007/s00468-016-1465-7

Ogawa K (2018) Mathematical consideration of the age-related decline in leaf biomass in forest stands under the self-thinning law. Ecol Model 372:64-69

Ogawa K (2019) Scaling relations based on the geometric and metabolic theories in woody plant species: A review. Perspect Plant Ecol Evol Syst 40:125480

Price C. A, Ogle K, White E. P, Weitz J. S (2009) Evaluating scaling models in biology using hierarchical Bayesian approaches. Ecol Lett 12:641-651. https://doi.org/10.1111/j.1461-0248.2009.01316.x

Schulze E. D, Valentini R, Bouriaud O (2021) The role of net ecosystem productivity and of inventories in climate change research: the need for "net ecosystem productivity with harvest". NEPH For Ecosyst 8:15. https://doi.org/10.1186/s40663-021-00294-z

Seber GAF, Wild CJ (1989) Nonlinear regression. John Wiley & Sons

Shimada H, Mie-ken (2010) Construction of yield tables for sugi (Cryptomeria japonica) and hinoki (Chamaecyparis obtusa) plantations applied to long-rotation management in Mie [Japan] Prefecture theory of the C-D effect. For Res Inst Tsu (Japan) 2:1-28

Shinozaki K, Kira T (1956) Intraspecific competition among higher plants. VII. Logistic theory of the C-D effect. J Inst Polytechnics 12:69-82

Tadaki Y, Kawasaki Y (1966) Studies on the Production Structure of Forest. IX. Primary Productivity of a Young Cryptomeria Plantation with Excessively High Stand Density. J Japan For Soc 48:55-61

Usoltsev VA (2013) Forest biomass and primary production database for Eurasia. CD-version. V.A. Usoltsev; [sci. ed. S.G. Shiyatov]; Ural State Forest Engineering University, Yekaterinburg

Weiner J (1985) Size Hierarchies in Experimental Populations of Annual Plants. Ecology 66(3):743-752

Will RE, Munger GT, Zhang Y, Borders BE (2002) Effects of annual fertilization and complete competition control on current annual increment, foliar development, and growth efficiency of different aged Pinus taeda stands. Can J For Res 32(10):1728-1740

Xue L, Hagihara A, Tadaki Y, Du X (1997) Self-thinning in Prince Rupprecht's Larch (Larix principis-rupprechtii Mayr) Stands. J For Res 2(4):221-226
Yoda K, Kira T, Ogawa H, Hozami K (1963) Self thinning in overcrowded pure stands under cultivated and natural conditions. J Biol Osaka City Univ 14:107–129

Yoshihara T, Matsumura H, nosuke Hashida S, Nakaya K (2016) Radiocesium contamination in living and dead foliar parts of Japanese cedar during 2011–2015. J Environ Radioact 164:291–299

Yoshihara T, Matsumura H, Tsuzaki M, Wakamatsu T, Kobayashi T, nosuke Hashida S, Nagaoka T, Goto F (2014) Changes in radiocesium contamination from Fukushima in foliar parts of 10 common tree species in Japan between 2011 and 2013. J Environ Radioact 138:220–226

Zapata-Cuartas M, Sierra C. A, Alleman L (2012) Probability distribution of allometric coefficients and Bayesian estimation of aboveground tree biomass, Forest Ecology and Management. ISSN 0378-1127 277:173–179.
https://doi.org/10.1016/j.foreco.2012.04.030

Submit your manuscript to a **SpringerOpen**
journal and benefit from:

► Convenient online submission

► Rigorous peer review

► Open access: articles freely available online

► High visibility within the field

► Retaining the copyright to your article

Submit your next manuscript at ► springeropen.com