Predicting individual tree growth using stand-level simulation, diameter distribution, and Bayesian calibration

Abstract

Key message

We propose a methodology to develop a preliminary version of a growth model when tree-level growth data are unavailable. This modelling approach predicts individual tree growth using only one-time observations from temporary plots. A virtual dataset was generated by linking the whole stand and diameter distribution models. The individual tree model was parameterized using Bayesian calibration and a likelihood of diameter distributions.

Context

A key component of tree-level growth and yield prediction is the diameter increment model that requires at least two different points in time with individual tree measurements. In some cases, however, sufficient inventory data from remeasured permanent or semitemporary plots are unavailable or difficult to access.

Aims

The purpose of this study was to propose a three-stage approach for modelling individual tree diameter growth based on temporary plots.

Methods

The first stage is to predict stand dynamics at 5-year intervals based on stand-level resource inventory data. The second stage is to simulate diameter distribution at 5-year intervals using a Weibull function based on tree-level research inventory data. The final stage is to improve the reliability of individual tree diameter estimates by updating parameters with Bayesian calibration based on a likelihood of diameter distributions.

Results

The virtual-data-based diameter increment model provided logical patterns and reliable performances in both tree- and stand-level predictions. Although it underestimated the growth of suppressed trees compared with tree cores and remeasurements, this bias was negligible when aggregating tree-level simulations into stand-level growth predictions.

Conclusion

Our virtual-data-based modelling approach only requires one-time observations from temporary plots, but provide reliable predictions of stand- and tree-level growth when validated with tree cores and whole-stand models. This preliminary model can be updated in a Bayesian framework when growth data from tree cores or remeasurements were obtained.

1 Introduction

Efficient and accurate models for growth and yield are a fundamental tool in forest sciences, playing a key role in forest management, forest planning, ecological studies, or in fact any discipline within the field. Traditionally, large efforts have been placed to design and retrieve the necessary data, often from large databases or long-term observations in permanent plots. Several places have had these data available since long ago: there are long-term experimental plots being still surveyed in Germany that date back to the 1860s (Prodan 1965; Pretzsch 2010). As a second source, plots from national forest inventory (NFI) are a good alternative for data retrieval, although, may lack accuracy and detail as are not necessarily designed for growth and yield models. In this case, the large amount of data available may compensate these limitations.

However, what are the alternatives when the purpose of analysis is a new species or a new region with limitations in longitudinal data for single trees? For instance, despite its large role in forest resources, Brazil has only recently started a country-level NFI based on a standardized and systematized methodology, which for some time will provide a single measurement (Freitas et al. 2016; da Luz et al. 2018). The same would apply to large countries in Asia, such as Mongolia, with 11.3 million hectares of boreal forest, which has completed its first NFI recently (Altrell 2019). In Europe, despite the large experience in NFI, several countries lack permanent plots as a main source of data (Tomppo et al. 2010). In fact, according to the global forest resource assessment, 40 of the 99 tropical countries in the report did not have at least two comparable NFI measurements on permanent plots that could be the basis for individual tree growth models (FAO 2015; Romijn et al. 2015).

When sufficient data from permanent plots are not available, semitemporary or interval plots are useful to develop forest growth models. Those plots are established in stands of different development stages and measured at least in two points in time. However, local forest regulations often limit the establishment, collection, and availability of extensive plot data along time; but, while remeasurement data could be more difficult, one-time observations in resource inventory data from temporary plots are more commonly available. This type of inventory relies on temporary plot sampling that may provide databases for modelling sample whole-stand models (e.g., Vanclay 2010) and diameter distribution models (e.g., Cao 2004). The basic data units of whole-stand models are stand and site parameters. Consequently, their construction and simulation do not require tree-level information of stand dynamics. In contrast, size-class models use an entire tree class as the basic unit, and individual tree models rely on the information from each tree. This approach, although useful, presents a more restricted suitability for individual tree models.

The individual tree approach, although being more prescriptive in data needs, performs largely better to characterize growth under variable stand conditions and management practices (Weiskittel et al. 2011). For instance, concerning thinnings, a whole-stand model is limited since thinning treatments are designed to selectively remove individual trees, which leads to alteration of the residual stand structure. To address these limitations, an alternative use of temporary plots is to estimate previous stand conditions (e.g., stand basal area 5 years earlier) based on tree cores; thus, a timeline of stand attributes, along with the increment measurement from tree cores, can be used to fit individual tree models. This approach is applied, for example by several versions of the ORGANON (ORegon Growth ANalysis and projectiON) growth model (Hann et al. 2011) and many variants of the western US’ Forest Vegetation Simulator (FVS) growth model (Crookston and Dixon 2005).

Another method is to incrementally modify the explanatory variables of tree vigor and competition (Schröder et al. 2002). If explanatory variables change slowly over time, errors introduced by not having quantitative backdated measures are negligible over a short interval. However, large datasets of tree cores to construct these models are often not available in most forest inventories, and although alternative approaches have been developed to model diameter increments, these are much less accurate (Vanclay 1994).

Limited by the availability of data, most modelling efforts are quite opportunistic and often commence with any data available (Weiskittel et al. 2011). In our case, we investigate the growth and yield of Pinus tabuliformis, a dominant conifer species native to China which plays an important ecological role in soil erosion and water protection and on the regional socioeconomic development (Wu and Feng 1994), with a wide distribution area and a large altitudinal gradient (Mao and Wang 2011). In recent years, adequate designs of forest resource management, in terms of tending, thinning, timber stand management, and bush management, have been needed for restoring the degraded forest and reducing the risk of forest fires. The lack of proper individual tree level growth models due to data limitations, and the ecological and economic importance of the species was the rationale for our research efforts to provide suitable modelling approaches that could be used in forest planning related to the species. Thus, we aim to establish a modelling approach to predict individual tree diameter growth based on one-time observation data from temporary plots, with special focus on the assessment of the effects on model performance of different data types and the resulting parametric uncertainty.

2 Materials and methods

We propose a three-stage approach to generate virtual data for fitting individual tree growth models by linking stand-level and size-class simulations (Fig. 1). First, we construct a whole-stand model, and then predict stand-level dynamics. Second, we simulate diameter distributions using the Weibull function that is parameterized with observed size-class structure data. Finally, we use the simulated data to fit individual tree growth models with Bayesian calibration using a diameter distribution likelihood. All data used to parameterize models in our study were collected from temporary plots. We collected one reference dataset of tree cores and one small dataset of diameter remeasurements to validate our simulated data modelling approach.

Fig. 1

Flowchart of the main processes in the individual tree diameter growth modelling

2.1 Study area

Temporary plot data were collected in the Huoditang Forest (108° 21′–108° 39′ E, 33° 18′–33° 28′ N) and the Xunyangba Forest (103° 58′–109° 48′ E, 32° 29′–33° 13′ N), in Ningshan County, Shaanxi, China (Fig. 2). The altitude of all sampling sites varied from 1200 to 1800 m, with a slope gradient ranging from 15 to 42°. The dominant soil is a brown forest soil with an average thickness of about 50 cm. This region has a subtropical humid mountain climate. The mean annual temperature is 8–12 °C, and the mean annual precipitation is 1100 mm. The primeval forests in this region were harvested during the 1960s and 1970s. Now most of the area is covered with secondary forests or plantations.

Fig. 2

Locations of the sampling plots. The 90m-DEM (Digital Elevation Model) was downloaded from the CGIAR-CSI GeoPortal (http://srtm.csi.cgiar.org). The locations only include plots being established after 2005. Plots of other datasets were from the same compartments but lacking exact coordinates

2.2 Modelling data

The data for model parameterization were collected from 215, variable-radius plots using angle gauge sampling and 22 square plots with a fixed area of 0.04 ha (Table 1). The basal area factor of angle gauge sampling is 1 m² ha⁻¹. This dataset was provided by the forest resource inventory (1990, 2005, and 2015); most of which comes from level II forest inventory in China (Lei et al. 2009). The level II inventory is designed for forest management planning, while most plots from 1990 and 2005 lack exact coordinates. All data are from temporary plots established in Pinus tabuliformis plantations from which only one-time observations were collected (Table 2). One-quarter of plots were randomly withheld for benchmarking, while the remaining three quarters were used for model fitting. Additional data of stem analysis from 11 medium trees across various site types and conditions were combined with plot sample data to calibrate the height-age equation in the whole stand model and the height-diameter equation in the individual tree model, taking into consideration that stem analysis data and plot data perform differently in height curve estimation (Perin et al. 2013).

Table 1 Characteristics of the two types of plot data used in modelling

Table 2 Posterior probability distribution of the parameters in the whole-stand model. The posterior distribution is not analytical and it is characterized here by the mean values of the parameters, the standard deviation (SD), and the maximum a posteriori estimates (MAP). The prior was set as independent uniform distributions, of which the maxima and minima are listed

For each size class (4 cm) in each square plot, increment measures from one or two sample trees were selected randomly, taking two cores in opposite directions at the breast height. This provided 338 tree cores from 169 trees in the 22 square plots. Increment core data were processed specifically to validate the individual tree diameter increment model (Fig. 1). Therefore, individual tree diameters and stand total basal areas from 5 years earlier were back-dated based on tree cores and the successional stages of snags. Besides the tree cores, another dataset was also obtained for model validation (Fig. 1). Trees from 18 of the 22 square plots were remeasured during 2019. The other four plots were unable to locate because of inaccurate coordinate records.

2.3 Model simulations

We generated a virtual dataset for fitting an individual tree model using the following six steps: (1) constructed a whole-stand model based entirely on temporary plot data, (2) using this model, we simulated the development of stand-level attributes for 30, virtual, 1-ha plots from ages 20 to 80 years that include six site fertility classes and five initial densities (6 × 5 = 30 plots), (3) we constructed a diameter distribution model based on temporary plot data, (4) using the diameter distribution model, we simulated the diameter distributions of virtual plots at 5-year intervals, (5) based on diameter distributions and stand density, we generated individual tree diameters, and finally (6), we simulated 30 virtual 1-ha plots among sites and competition conditions, with 13 ((80–20)/5 + 1) measurement times per plot.

2.3.1 Whole-stand model

A whole-stand model that estimates stand average height, per hectare basal area, and quadratic mean diameter was built and fitted using temporary plots data. The following basic functions were hypothesized and fit:

$$ {\overline{h}}_{gc}=\exp \left({\alpha}_1+{\alpha}_2{t}^{\alpha_3}\right) $$

(1)

$$ {\overline{h}}_{30}=\exp \left({\alpha}_1+\left(\ln \overline{h}-{\alpha}_1\right)\cdotp {\left(\frac{30}{t}\right)}^{\alpha_3}\right) $$

(2)

$$ \overline{h}=\exp \left({\alpha}_1+\left(\ln {\overline{h}}_{30}-{\alpha}_1\right)\cdotp {\left(\frac{t}{30}\right)}^{\alpha_3}\right) $$

(3)

$$ SDI=N\times {\left(\frac{{\overline{d}}_g}{20}\right)}^{\alpha_4} $$

(4)

$$ G={\alpha}_5\times {{\overline{h}}_{30}}^{\alpha_6}\times \exp \left(-\frac{\alpha_7\times {\left(\frac{SDI}{1000}\right)}^{-{\alpha}_8}}{t}\right) $$

(5)

$$ {\overline{d}}_g={\alpha}_9\times {{\overline{h}}_{30}}^{\alpha_{10}}\times \exp \left(-\frac{\alpha_{11}\times {\left(\frac{SDI}{1000}\right)}^{\alpha_{12}}}{t}\right) $$

(6)

$$ N=\frac{40000\times G}{\pi \times {{\overline{d}}_g}^2} $$

(7)

where t is the stand age (years), ${\overline{h}}_{gc}$ is the average stand height in guide curve (m), and ${\overline{h}}_{30}$ is average height at the reference age (m), of which the reference age is 30 (years) for Pinus tabuliformis. $\overline{h}$ is stand height (m), SDI is stand density index (trees ha⁻¹), N is stand density (trees ha⁻¹), G is stand basal area (m² ha⁻¹), and ${\overline{d}}_g$ is quadratic mean diameter (cm).

There was a lack of records of dominant trees in previous inventory datasets. Therefore, the expected stand average height at the reference age, ${\overline{h}}_{30}$, was used to track site quality. We assumed that, similar to the widely used site index, the relationship between stand age and stand average height correlates closely with total stand production and is less dependent on thinning intensity (Eichhorn 1902). Site fertility classification can also be implemented based on average height and age relations (von Baur 1876; de Perthuis Laillevault 1803), and this site evaluation method has been widely applied in Chinese forest inventory since the 1950s. Tree height growth has been found to be impacted by stand density (Lynch 1958; Cieszewski and Bella 1993; MacFarlane et al. 2000; Flewelling et al. 2001), and the severity of this impact is related to species tolerance and site conditions (Weiskittel et al. 2011). Furthermore, the intensification of thinning from below can significantly affect the calculation of the stand average height (Assmann 1970; Pretzsch 2010). Based on this, the limitation of using average height was not considered as a severe issue in our study, because thinning was forbidden after 1998 in this region due to the National Forest Conservation Programme (NFCP) for soil and water protection (Li 2004). A guide curve of stand average height (Eq. 1) was fitted with the Bailey–Clutter model (Bailey and Clutter 1974). The equations of ${\overline{h}}_{30}$ and stand height (Eqs. 2 and 3) were then derived.

Reineke’s (1933) stand density index (SDI), which is fitted by regression of number of stems (N) and quadratic mean diameter (${\overline{d}}_g$) using fully stocked stand observations, was used to describe competition (Eq. 4). Stand total basal area (G) model (Eq. 5) and quadratic mean diameter model (Eq. 6) were modifications of the Schumacher model (Schumacher 1939), and number of stems was derived (Eq. 7).

Full-stocked stands were classified using the following three steps: (1) fit the Reineke’s (1933) equation lnN = λ₁ + λ₂ ln${\overline{d}}_g$ with all stand observations, (2) drop all points that fall below this regression line and repeat the regression, and (3) drop all points falling below the new regression line. The remaining was considered fully stocked. We also validated this method with local stock tables of “normal” stands.

This whole stand model was chosen because it can be fitted with temporary plots, although, it may be biased depending on the data distribution: the range of site fertility classes, and stand densities must be covered in all age classes. This model system was used for generating variable-density yield tables with given SDI values but not sufficient for simulating a specific real stand. Therefore, when using this whole-stand model, we set the SDI and ${\overline{h}}_{30}$ manually and then ran the ${\overline{d}}_g$ and G equations (Eqs. 5 and 6) to estimate the dynamics of the virtual stands.

2.3.2 Diameter distribution model

The Weibull distribution was used to characterize the diameter distribution. This probability density function is calculated as:

$$ f\left(x;{\beta}_1,{\beta}_2,{\beta}_3\right)=\left(\frac{\beta_3}{\beta_2}\right){\left(\frac{x-{\beta}_1}{\beta_2}\right)}^{\beta_3-1}\exp \left[-{\left(\frac{x-{\beta}_1}{\beta_2}\right)}^{\beta_3}\right] $$

(8)

where β₁, β₂, and β₃ are the location, scale, and shape parameters of Weibull distribution, respectively, and x is tree diameter at breast height. We used a moment-based parameter recovery method to calculate scale and shape parameters from quadratic mean diameter and diameter variance. The regression equation of diameter variance is:

$$ {D}_{\mathrm{var}}=\exp \left({\gamma}_1+{\gamma}_2\ln {\overline{d}}_g+{\gamma}_3\ln G\right) $$

(9)

where D_var is diameter variance. The location parameter β₁, lower limit of x, was similarly set as:

$$ {\beta}_1=\exp \left({\gamma}_4+{\gamma}_5\ln {\overline{d}}_g+{\gamma}_5\ln G\right) $$

(10)

Parameter γ_is was obtained by maximum likelihood estimator regression (Cao 2004). The β₂ and β₃were recovered by:

$$ {\beta}_2=-{\beta}_1{G}_1/{G}_2+{\left[{\left({\beta}_1/{G}_2\right)}^2\left({G_1}^2-{G}_2\right)+{D}^2/{G}_2\right]}^{0.5} $$

(11)

$$ {\beta_2}^2\left({G}_2-{G_1}^2\right)-{D}_{\mathrm{var}}=0 $$

(12)

where G_i is Γ(1 + i/β₃); Γ(∙) is the complete gamma function. The diameter distribution model was developed using only square plot data because the size-class information of variable radius plots was inadequate.

2.4 Individual tree diameter model

An individual tree model simplified from Pukkala et al. (2009) was fitted using Bayesian calibration and simulated virtual data. This model system includes a distance-independent diameter growth, survival, and height prediction models:

$$ {i}_d=\exp \left({c}_1+{c}_2{G}_{>d}+{c}_3\ln G+{c}_4\sqrt{d}+{c}_5{d}^2+{c}_6{\overline{h}}_{30}\right) $$

(13)

$$ p=\frac{1}{1+\exp \left[-\left({c}_7+{c}_8\sqrt{d}+{c}_9{G}_{>d}\right)\right]} $$

(14)

$$ h=\frac{c_{10}+{c}_{11}{\overline{h}}_{30}}{1+\left({c}_{12}/d\right)+\left({c}_{13}/{d}^2\right)} $$

(15)

where i_dis 5-year diameter increment (cm), G_>d is basal area in larger trees (m² ha⁻¹), G is stand total basal area (m² ha⁻¹), d is individual tree diameter at breast height (cm), ${\overline{h}}_{30}$ represents site fertility class (m), p is the probability to survive for the coming 5 years, h is individual tree height (m).

The virtual data simulated in this study were designed for the diameter increment and survival models. The ingrowth model was deliberately omitted because of ingrowth stochasticity in pure, even-aged Pinus tabuliformis plantations. Therefore, recruitment was assumed to be null. The height projection model is static and can be fitted separately from the diameter increment and survival equations.

2.5 Parameter estimation and uncertainty analysis

The Bayes’ theorem offers an approach for updating distributions, and can be applied repeatedly when new data are available. In its simplified notion, this theorem states:

$$ p\left(\uptheta |D\right)\propto p\left(D|\uptheta \right)\times \mathrm{p}\left(\uptheta \right) $$

(16)

where p(θ| D)and p(θ) are the posterior and prior distribution for the parameters θ, respectively, and p(D| θ) is the likelihood function, which is the probability of the data D for a given θ. For the first calibration of physiological process-based model (e.g., Van Oijen et al. 2005), prior distribution information is usually based on best guesses or published values. In this study, empirical model parameters have no physiological meaning. Thus, vague priors were used for the first calibration. Parameters were assumed as independent uniform distributions, of which the minima and maxima are shown in Tables 2, 3, and 4.

Table 3 Posterior probability distribution of the parameters in the individual tree diameter increment model. The posterior distribution is not analytical and it is characterized here by the mean values of the parameters, the standard deviation (SD) and the maximum a posteriori estimates (MAP). The prior was set as independent uniform distributions, of which the maxima and minima are listed

Table 4 Posterior probability distribution of the parameters in the individual tree survival and height prediction models. The prior was set as independent uniform distributions, of which the maxima and minima are listed

Because multiple sources of data were used in the Bayesian calibration, the likelihood function was decomposed based on the assumptions of independent errors. For example, the calibration frame of height model (Eqs. 1 and 15) was:

$$ p\left(\uptheta |D\right)\propto p\left({D}_{VR}|\uptheta \right)\times p\left({D}_{SP}|\uptheta \right)\times p\left({D}_{SA}|\uptheta \right)\times p\left(\uptheta \right) $$

(17)

where D_VR, D_SP, andD_SA are observations from variable-radius plots, square plots, and stem analysis respectively. A heavy tailed likelihood probability density function (Sivia and Skilling 2006) was used:

$$ p\left(D|\theta \right)={\prod}_{i=1}^n\frac{1}{\sigma_i\sqrt{2\pi }}\left[\frac{1-\exp \left(-{R_i}^2/2\right)}{{R_i}^2}\right] $$

(18)

where σ_iis the noise-scaling factor; a measure of the uncertainty of random error of the ith observation. R_i is the residual for the ith datum and prediction, divided by σ_i.

Data from variable-radius plots, square plots, and stem analyses were not mixed since the noise-scaling factor, σ,is different for each dataset. For each data type, a linear relation between error and prediction was established as ${\sigma}_{mn}={a}_m\cdotp {\hat{y}}_{mn}+{b}_m$, where ${\hat{y}}_{mn}$ is the nth prediction of dataset m. Parameters a and b were calibrated along with parameters θ.

A special likelihood function modified from Pukkala et al. (2011) was used to simultaneously fit individual tree diameter growth model and survival model:

$$ {\displaystyle \begin{array}{c}p\left(D|\theta \right)=p\left(E=D-M\left(\uptheta \right)\right)=\\ {}{\prod}_{k=1}^K{\prod}_{j=1}^{J_k}{\prod}_{i=1}^{I_j}\left[\upvarphi \left({G}_{jk}^D\left({d}_{ijk}\right)-{G}_{jk}^M\left({d}_{ijk}\right);0,{\sigma_{ijk}^{G^2}}\right)\cdotp \upvarphi \left({F}_{jk}^D\left({d}_{ijk}\right)-{F}_{jk}^M\left({d}_{ijk}\right);0,{\sigma_{ijk}^{F^2}}\right)\right]\end{array}} $$

(19)

where θ is the set of coefficients of the individual tree models, K is number of plots, J_k is the number of measurement intervals of plot k, I_j is number of 4-cm diameter classes in measurement interval j of plot k, φ denotes a Gaussian probability density function with a given mean and variance,${G}_{jk}^D\left({d}_{ijk}\right)$ and ${G}_{jk}^M\left({d}_{ijk}\right)$ are, respectively, cumulative basal area (m² ha⁻¹) of data and prediction at diameter d_ijk (upper limit of diameter class i) at the end of measurement interval j of plot k, and ${F}_{jk}^D\left({d}_{ijk}\right)$ and ${F}_{jk}^M\left({d}_{ijk}\right)$ are, respectively, cumulative number of trees per hectare of data and prediction at diameter d_ijk at the end of measurement interval j of plot k.

Detailed information on diameter increments and survival were unavailable from the simulated data because of an identification problem: we did not specifically know which trees die or individual tree increments during the 5-year period because trees in virtual plots were not identified. Thus, diameter increment and survival models (Eqs. 13 and 14) were fitted together using the likelihood function (Eq. 19) based on diameter distribution changes. The height prediction model (Eq. 15) was parameterized independently.

Pukkala et al. (2011) used least-squares to balance two variables at different scales (G and F in Eq. 19) by multiplying the sum of squared errors of F by 0.0001. However, this choice could be entirely arbitrary and should be justified. We converted the objective function of the optimization-based modelling into a likelihood function. Linear relations between the standard deviation and the model predictions were calculated as ${\sigma}_{ijk}^G={a}_G\cdotp {G}_{jk}^M\left({d}_{ijk}\right)+{b}_G$ and ${\sigma}_{ijk}^F={a}_F\cdotp {F}_{jk}^M\left({d}_{ijk}\right)+{b}_F$, where parameters a_G, b_G, a_F, b_F were calibrated together with model parameters. Here, two variables, G and F, will have two different linear relations. In this way, variable error increased with its magnitude and the likelihood automatically balanced the scales of the two variables.

Unlike the likelihood function for fitting whole stand model or the individual tree model, the likelihood for diameter distribution model was not based on the above assumptions of error distributions, but on the probability density function of Weibull distribution:

$$ p\left(D|\theta \right)={\prod}_{i=1}^p{\prod}_{j=1}^{n_i}\left(\frac{\beta_3}{\beta_2}\right){\left(\frac{x_{ij}-{\beta}_1}{\beta_2}\right)}^{\beta_3-1}\exp \left[-{\left(\frac{x_{ij}-{\beta}_1}{\beta_2}\right)}^{\beta_3}\right] $$

(20)

where p is the number of plots, n_i is the number of trees in the ith plot, x_ij is the diameter at breast height of tree j in the ith plot. The regression equations for D_var and β₁(Eqs. 9 and 10) were simultaneously fitted while β₂ and β₃ were recovered (Eqs. 11 and 12).

The MCMC was simulated, using differential evolution adaptive Metropolis with snooker updating (DREAMzs), which runs a few chains in parallel and explores the parameter space in an efficient way (Laloy and Vrugt 2012; Vrugt et al. 2009). We used the DREAMzs algorithm implemented in the R package BayesianTools (Hartig et al. 2017; R Development Core Team 2017). When using unanalytical joint posterior distribution as the prior for a recalibration, a multivariate sample was extracted from the previous MCMC. With this sample, a multivariate normal density was fitted and used as the prior for the new calibration.

The MCMC convergence diagnostic (Brooks and Gelman 1998; Gelman and Rubin 1992) was used to monitor the convergence in the MCMC output. The multivariate potential scale reduction factor (MPSRF) was calculated, based on two MCMC runs, each of which has three internal chains. The number of iterations was set based on the complexity of the model and data, ranging from 10⁵ to 10⁷. A large MPSRF means that the output from all chains is distinguishable, and a notable difference exists between variance and intrachain variance. In our study, convergence was diagnosed when the MPSRF was below 1.05, which is a relatively strict criterion (Brooks and Gelman 1998).

2.6 Model validation

One-quarter of plot observations were randomly selected for benchmarking, including the whole-stand, diameter distribution, and individual tree height prediction models. The increment core observations and remeasurements were respectively processed to validate the individual tree diameter increment equation. For the tree core data, we back-dated tree and stand attributes to the start of previous 5-year growth period, to create explanatory variables. The inside bark diameters from 5-year earlier were calculated for each subject tree by subtracting the 5-year inside bark increment. Bark thicknesses were estimated based on the assumption that ratios between outside and inside bark diameters were constant over a relatively short, 5-year period (e.g., Pukkala 1989). Diameter growth rate of each size class was calculated using increment cores. Tree diameters from 5-years previous were roughly estimated using estimated growth rates. From this we obtained the 5-year earlier stand basal area and basal area in larger trees. For the remeasurements, the interval years were calculated based on the number of growing seasons, considering that diameter growth mainly occurs in the early summer. Growth data of 4-, 6- or 7-year interval were linearly converted to the growth of 5 years.

We used a regression-based equivalence test developed by Robinson et al. (2005) for model benchmarking. Unlike usual parametric methods, equivalence tests rely on a null hypothesis that a model is unacceptable rather than acceptable. The equivalence test not only compares means but also for checks for similarity between individual predictions and observations. The validation strategy suggested by Robinson et al. (2005) followed the following six steps: (1) subtract mean prediction from component predictions, (2) establish regions of equivalence for the shifted intercept and slope (e.g.,±25%), (3) fit a linear regression between observations and adjusted predictions, (4) test the intercept for equality by calculating two, one-sided confidence intervals for the intercept and comparing these to the estimated equivalence regions, (5) test the slope for equality by calculating the two, one-sided confidence intervals for the slope and comparing to its equivalence regions, (6) accept or reject the hypothesis of dissimilarity based on whether the confidence interval falls within the equivalence region. The R module “equivalence” was used to perform the regression-based equivalence test (Robinson 2016). We calculated the minimum percentage that would result in rejection of the null hypothesis (Froese and Robinson 2007).

The individual tree model and the diameter increment equation (Eq. 13) were fitted separately by simulated virtual data, increment core data, and remeasurement data (Fig. 1). The performance of the core-data-based (CB) model, remeasurement-based model (RB), and virtual-data-based (VB) model was compared at both tree- and stand-levels. In addition to the equivalence test, mismatches between real growth observations and VB model predictions were quantified by decomposed mean squared deviation (MSD), which identified causes of possible large deviations between observations and predictions (Kobayashi and Salam 2000). As suggested by Gauch et al. (2003) the MSD was partitioned into three parts: squared bias (SB), nonunity slope (NU), and lack of correlation (LC). When comparing observations and predictions, squared bias means the deviation of means; nonunity slope is relevant to the slope of linear regression; lack of correlation means a scatter. In this case, we expected the mean of model predictions and the mean of observations to be the same, and the slope of observations on predictions to be 1. Thus, the data model mismatch was decomposed as:

$$ \mathrm{MSD}=\frac{\sum_{i=1}^n{\left({X}_n-{Y}_n\right)}^2}{N}=\mathrm{SB}+\mathrm{NU}+\mathrm{LC} $$

(21)

$$ \mathrm{SB}={\left(\overline{X}-\overline{Y}\right)}^2 $$

(22)

$$ \mathrm{NU}={\left(1-\beta \right)}^2\left(\sum \frac{{\left({X}_n-\overline{X}\right)}^2}{N}\right) $$

(23)

$$ \mathrm{LC}=\left(1-{r}^2\right)\left(\sum \frac{{\left({Y}_n-\overline{Y}\right)}^2}{N}\right) $$

(24)

where $\overline{X}$ and $\overline{Y}$ are the means of model predictions (X) and observations (Y), respectively, β is the slope of the least-squares regression of Y on X, r² is the square of the correlation coefficient, and N is the number of observations.

3 Results

3.1 Whole-stand and diameter distribution models

In the Bayesian calibration process, the joint posterior distribution of parameters is approximated by Markov Chain Monte Carlo simulation. The mean value, standard deviation, and MAP (maximum a posteriori probability estimates) of parameters in the whole-stand model are given in Table 2. As shown in the model of guide curves, different types of inventory data resulted in contrasting predictions and joint parametric distributions (Fig. 3). The performance of the guide curve estimated from multisource data was similar to that of the curve based on various-radius plots. However, the joint posterior distribution of a₁ and a₃ from multisource distribution was more similar to the distribution from equation based on stem analysis, because the likelihood weighted each dataset automatically based on the estimated error distributions. Based on the whole-stand model, we simulated the growth and mortality of 30, 1-ha virtual plots (Fig. 4). These plantations possessed the characteristics of both relatively high initial stand basal area (Fig. 4a) and mortality (Fig. 4b), since no thinning schemes were designed in the virtual plots.

Fig. 3

Guide curves of stand average height and joint distributions of parameters in the guide curve equation (Eq. 1) when based on different types of data. The grey area is the 95% Bayesian credible interval based on parameter uncertainty

Fig. 4

a Development of stand total basal area on different site fertility classes. b Development of stand density with different initial density

The parameters that predict diameter variance (Eq. 9) and the minimum diameter (Eq. 10) are given in Table 2. The diameter distribution changes due to growth and self-thinning in a medium site and density condition (Fig. 5). The most prominent changes were tree frequencies in small size classes, especially at young ages. The minimum diameter increased over time.

Fig. 5

Development of diameter distribution at different ages. The probability density functions were for stands with 2000 trees ha⁻¹ initial density at age 20 and ${\overline{h}}_{30}$=12 m

3.2 Individual tree diameter model

The calibrated parameters of the individual tree model are given in Table 3. Parameter vectors based on simulated virtual data, increment core method, and remeasurements were logical and general similar. Compared with the CB and RB models, the diameter increment predictions of the VB model were more sensitive to their own diameters (Fig. 6). The diameter growth rate reached its maximum at about 15–20 cm for P. tabuliformis. Simulated VB model was most sensitive to competition variable BAL. The RB model was most sensitive to the total stand basal area and site fertility. The standard deviations in Table 3 represented the parametric uncertainty; the standard deviations of parameters for lnTBA and ${\overline{h}}_{30}$ were larger than those of BAL, which suggests more uncertainty in predictions. Using the parameters’ joint posterior distribution of the VB model as the prior, we recalibrated the model using tree cores and remeasurements, which combined information of both virtual data and real growth observations in a Bayesian framework (Eq. 16). The parameter vector of the recalibrated model was somewhat a compromise among parameter vectors of the VB, CB, and RB models (Table 3).

Fig. 6

Comparison of individual tree diameter increment predictions based on simulated virtual data, increment cores, and remeasurements. (a The relationships of 5-year diameter increment over diameter, b stand basal area, c basal area of larger trees (G_>d), d and site fertility (${\overline{h}}_{30}$, average height on the reference age). The dotted, dashed, and solid lines are, respectively, virtual-data-based predictions, increment core method predictions, and remeasurement-based predictions. Lines were generated by MAP (the maximum a posterior parameter vector). The grey area represents the 95% Bayesian credible interval based on parametric uncertainty of virtual data method. For Fig. 5a, the stand total basal area is 25 m² ha⁻¹, the basal area of larger trees is 5 m² ha⁻¹, and the ${\overline{h}}_{30}$is 12 m. For Fig. 5b, the diameter is 20 cm, the basal area of larger trees is 5 m² ha⁻¹, and the ${\overline{h}}_{30}$is 12 m. For Fig. 5c, the diameter is 20 cm, the stand basal area is 35 m² ha⁻¹, and the ${\overline{h}}_{30}$is 12 m. For Fig. 5d, the diameter is 15 cm, the stand basal area is 20 m² ha⁻¹, and the basal area of larger trees is 5 m² ha⁻¹)

Survival model coefficients (Table 4) show that the survival probability increased with increasing diameter or decreasing competition (Fig. 7). When the diameter of a tree reached 30 cm, it was unlikely to die in subsequent measurement intervals, because the mortality from senescence and disturbances were not considered in the survival model. The height prediction model indicated trees of equivalent diameter are higher in more fertile sites (Fig. 8), and height uncertainty increases as the diameter increases.

Fig. 7

Dependence of the probability that a tree to survives for 5 years relative to diameter at breast height. The solid line, dashed line and dotted line represent trees on G_>d of 5, 15, and 25 m² ha⁻¹, respectively. The stand total basal area is 30 m² ha⁻¹, and the ${\overline{h}}_{30}$ is 12 m. The grey area is the 95% Bayesian credible interval based on parameter uncertainty

Fig. 8

Dependence of the height of Pinus tabuliformis on diameter at breast height. The solid line, dashed line and dotted line represent trees on good quality site (${\overline{h}}_{30}$=14 m), median site (${\overline{h}}_{30}$=12 m), and poor site (${\overline{h}}_{30}$=10 m), respectively. The grey area is the 95% Bayesian credible interval based on parameter uncertainty

3.3 Validation of different resolution models

Tests for the intercept only failed to reject the null hypothesis of dissimilarity in the individual tree diameter increment model (Table 5). However, tests for the slope failed to reject the null of quadratic mean diameter prediction function, diameter distribution function, and individual tree diameter increment model. The stand-level basal area model showed the best performance in both intercept and slope tests. Since the area of the square plot was only 0.04 ha, random errors in Weibull distributions could be large when compared with observations. The minimum rejection interval of the slope in dbh frequencies was as large as 80.25% (Table 5), which was caused by the randomness of observed diameter distribution on a limited population. As the difference between virtual data and increment CB models shows (Fig. 6), equivalence tests for individual tree diameter increment model also largely failed to reject null hypotheses of dissimilarity (Table 5).

Table 5 Summary of equivalence-based regression results. Sample size is denoted by n, and the approximate joint two one-sided 95% confidence intervals for the intercept β₀ and the slope β₁ are (${C}_{\beta}^{-},{C}_{\beta}^{+}\Big)$, the intercept interval of equivalence (${I}_{\beta}^{-},{I}_{\beta}^{+}$) for intercept and slope are, respectively, $\overline{y}\pm$25% and 1±0.25. F (8, 9,10) means frequency in specific size class, which were predicted using diameter distribution model (Eqs. 8, 9, and 10)

The virtual data-based individual tree diameter increment model underestimated the growth of the suppressed trees compared with the tree core observations (Fig. 9a) or remeasurements (Fig. 9b). The VB model was more sensitive to the competition variable G_>d than the CB and RB models (Fig. 6c and Table 3). When validated with tree core data, the MSD of the VB model was as high as 1.92 times of MSD of the CB model (Fig. 10). When validated with remeasurements, the MSD of the VB model was as high as 1.38 times of MSD of RB model. The distinction between the CB and RB model was smaller.

Fig. 9

Simulated and observed 5-year individual tree diameter increments. Comparisons between the virtual-data-based model and a tree cores or b remeasurements

Fig. 10

Decomposed mean squared deviation (MSD) between predictions and the observations of a tree cores or b remeasurements. VB = virtual-data-based individual tree diameter growth model; CB = core-data-based individual tree diameter growth model; RB = remeasurement-based individual tree diameter growth model

On the individual tree scale, this model-data mismatch means biased predictions for the highly suppressed trees (Fig. 9). However, this bias was negligible on the stand-level growth predictions (Fig. 11). The VB model shows the highest consistency with the whole-stand model.

Fig. 11

Comparison of basal area predictions between the whole-stand model and individual tree models. Initial states of 100 plots were randomly selected as input, and the stand total basal area after five years were respectively predicted by the whole stand model and three individual tree models. VB = virtual-data-based individual tree diameter growth model; CB = core-data-based individual tree diameter growth model; RB = remeasurement-based individual tree diameter growth model

4 Discussion

The principles of stand development have long been applied to model diameter growth, competition, stand structure, and diameter distribution in forests (Vanclay 1994). With remeasured data available, more complicated competition variables for modelling detailed individual tree models, such as FVS (Crookston and Dixon 2005) and ORGANON (Hann et al. 2011), are possible. Vanclay (2010) reported the feasibility of robust relationships to predict forest growth with sparse data at the stand level. At the tree level, we have applied a simple but robust competition index, G_>d, to demonstrate the sensitivity of competition effects on individual tree diameter growth. The empirical whole-stand and diameter distribution models were used as constraints to control the bias between the stand and the tree level. This makes the individual tree diameter model more reliable to track stand diameter structure and basal area development.

We used both tree core observations and tree-level remeasurements to validate the predicted individual tree diameter growth. The VB, CB, and RB models presented somewhat different predictions of diameter increment. The growth prediction from CB model was not sensitive to stand basal area, and, on the contrary, the RB model was largely sensitive to basal area, possibly due to the fact that both tree cores and remeasurements were from a limited number of plots. The large uncertainty resulting in the predictions of the RB model (Fig. 12) indicates that the remeasurement data were unable to sufficiently constrain the parameter distributions. These divergent results reflected the fact that these three methods sampled and simulated the stands in different ways. Nevertheless, the accuracy or biases were primarily dependent on database size, measurement error, and inventory sampling design rather than modelling techniques.

Fig. 12

Posterior uncertainty of predictions from three diameter increment models, by 3000 runs of each model. The diameter of the objective tree is 15 cm, the stand basal area is 30 m² ha⁻¹, basal area in larger trees is 10 m² ha⁻¹, and the ${\overline{h}}_{30}$is 12 m. Violin plots represent the density estimates of predictions from minimum to maximum values. The boxplots and white points at the median represent the quartiles of distributions. VB = virtual-data-based individual tree diameter growth model; CB = core-data-based individual tree diameter growth model; RB = remeasurement-based individual tree diameter growth model

One of the greatest problems in using data from temporary plots is that biases could be easily introduced into the model. The guide curve method assumes that average site productivity and site fertility distribution is the same in all age classes. If this is not the case, the model for stand height (Eq. 1) will be biased, affecting all other models. In practice, a temporary plot dataset with the same site fertility distributions in all age classes is almost impossible. Thus, we combined a small dataset of stem analysis in the parameterization of height prediction models by using multiplied likelihood (Eq. 17). Our dataset contains 237 temporary plots and only 11 trees with stem-analysis records. The temporary plot data were more representative, while the stem analysis data were more accurate by avoiding the biases of changing site distributions. By calibrating the values of noise-scaling factor σ, the weights of different datasets were balanced automatically.

The use of the parameter recovery estimator in a Weibull distribution linked the whole-stand and size-class models. The estimation of diameter distributions relies on stand-level variables, so variability of individuals will change consistently with the stand-level growth and mortality. The individual tree model based on virtual data will also provide consistent projections with stand-level attributes. This is because model data were simulated directly from the whole-stand model. This could be a distinct advantage when model users are concerned about accumulated error from aggregations of tree-level projections. However, systematic errors of whole-stand and diameter distribution models are also passed on to the individual tree model. In contrast, tree core–based individual tree models were developed independently of other models’ resolution. Thus, the CB model provided more reliable, tree-level projections that were well matched with tree core observations. However, the models that rely on tree cores or remeasurement may also provide unrealistic projections of stand-level attributes because changes in stand-level variables are not simulated directly but aggregated from tree-level projections (Ritchie and Hann 1997). In our case, the VB and CB models performed similarly in stand level simulations and were consistent with the whole-stand model predictions. However, the risk in the aggregation was revealed in the RB model. Several approaches can be used to modify stand-level projections by linking stand- and tree-level models. These include disaggregating and allocating stand-level growth and mortality into individual trees (e.g., Clutter and Allison 1974; Qin and Cao 2006), estimating multiresponse parameters to optimize tree-level predictions at multilevels (e.g., Zhang et al. 1997; Cao 2006) and using a composite estimator to joint estimates of stand- and tree-level models (e.g., Yue et al. 2008; Zhang et al. 2010).

With the optimization-based modelling method of Pukkala et al. (2011), we fit the individual tree model by minimizing the difference between measured and predicted diameter distributions. de Miguel et al. (2014) also found that this optimization-based method can solve tree misidentification and uneven measurement interval problems in individual tree models. Instead of the optimization algorithm, we chose Bayesian calibration for fitting models for two reasons; to quantify uncertainty of parameters and predictions, and to update models when new data became available. Considering that the predictions of the VB model were poorly matched with tree core observations, we recalibrated the model using the Bayesian calibration framework (Eq. 16). The posterior distributions of parameters in the VB model were taken as the prior information. Increment cores and remeasurements were used to calculate likelihood. The posterior distributions of parameters contained the information from both the simulated virtual stand and real growth measurements. At the end of this process, the predictions of the recalibrated model can be thought as a dynamic combination of the three datasets (Fig. 12). The basal area parameter (lnG) in the recalibrated model was greater than that in the CB model (Table 3), which means that the recalibrated model was more sensitive to the change in stand basal area. It is important to note that the posterior information of VB model was not completely converted into the new, prior information during recalibration. The posterior parameter distributions cannot be determined analytically. Therefore, we adopted the multivariate normal distribution, which was only a general description of the joint posterior distributions, to fit the density function.

We fit the individual tree diameter growth model and the individual tree survival model by linking stand-level and size-class simulations. Another important component, the ingrowth model, was deliberately omitted. This was because all data were collected from pure even-aged stands, and information of ingrowth from temporary plots was very limited. To validate our method, we used tree cores to backdate the stand condition from 5 years earlier and then fit another individual tree diameter growth model. Theoretically, an individual tree survival model could be constructed based on the backdated information. However, the exact time of tree death could not be estimated because of small age differences among trees and unclear deterioration characteristics.

More than two hundred plots were used to calibrate the whole-stand model. For other more limited cases or areas, other forms of whole-stand models or diameter distribution models could possibly be adopted. Vanclay (2010) developed a set of sample stand–level models that only contained three parameters and could be fitted with fewer data (even a single observation) for plantations. A variety of parameter recovery or parameter prediction methods are also available for simulating diameter distributions (e.g., Poudel and Cao 2013).

5 Conclusions

In this study, we proposed a three-stage modelling approach to forecast individual tree diameter growth using limited single observation data. The empirical whole stand model provides an accurate forecast of basal area development at the stand level. Links between the whole-stand model and diameter distribution model across forest sites and densities effectively enlarged the database of reliable stand structure information. This permits Bayesian calibration to model individual tree diameter without remeasured data. Our results demonstrated the individual tree diameter growth model effectively predicts competition at individual tree level; compared with empirical predictions based on incremental tree core measures, although, the simulated results seem to underestimate the diameter growth of suppressed trees.

The ultimate goal of constructing stand level models and diameter distributions is to provide adequate simulated data and appropriate likelihood functions for the calibration of individual tree models. Under the Bayesian framework, the proposed individual tree growth model can be updated when data from increment cores, interval plots or permanent plots are available, providing more accurate and robust predictions along time.