Simulating the diameter growth responses of Larix gmelini Rupr. and Betula platyphylla Suk. to biotic and abiotic factors in secondary forests in Northeast China

Abstract

Key message

The diameter growth of Dahurian larch (Larix gmelini Rupr.) and white birch (Betula platyphylla Suk.) species in secondary forest of Northeast China was not only influenced by biological factors such as tree size and stand characteristics, but also significantly affected by topographic and climatic factors such as temperature and precipitation. It is necessary to consider the abiotic factors in simulating the diameter growth.

Context

Climate change, such as global temperature rise, increased frequency of extreme weather events, and rising sea levels, has put forest ecosystems in an unstable state and has an impact on species composition, growth harvest, productivity and other functions of forests. And this impact varies in climate scenarios, regions and forest types.

Aims

To gain a comprehensive understanding of the adaptation for key species to their environment in secondary forests in Northeast China, the diameter growth responses of Dahurian larch and white birch to biotic and abiotic factors were simulated to assess the effects of climate on diameter growth.

Methods

China’s National Forest Continuous Inventory (NFCI) data from 2005 to 2015 were used to develop linear mixed-effects diameter growth models with plot-level random effects, and leave-one-out cross-validation was applied to evaluate the developed models. At the beginning of modeling, correlation analysis and best-subset regression were used to analyze the correlation between the diameter increment and the biotic and abiotic factors.

Results

(i) Sorting the categories of predictors in descending order based on the relative importance of the significant predictors, diameter growth of Dahurian larch was affected by competition, tree size, topographic conditions, stand attributes, diversity index, and climate factors, while the white birch species was affected by competition, tree size, stand attributes, climate factors, diversity index, and topographic conditions; (ii) the plot-level mixed-effects model, which achieved better fit and prediction performance than did basic linear models of individual-tree diameter growth in the cases of prediction calibration, was preferable for modeling individual-tree diameter growth; (iii) the prediction accuracy of the mixed-effects model increased gradually with increasing size of calibration sample, and the best sampling strategy was the use of nine random trees to calibrate and make predictions with the mixed-effects model for the larch and birch species; (iv) Dahurian larch was dominant in terms of interspecific competition, and the growth of this species was enhanced when it was grown with the birch.

Conclusion

In addition to biotic factors such as tree size and stand characteristics, the impact of climate on the growth of Dahurian larch and white birch should be considered in future management policies.

1 Introduction

A clear understanding of tree diameter growth not only helps to inform decision making in forest adaptation management but also helps to assess tree and stand carbon sinks. At present, individual-tree diameter growth models are effective tools for modeling growth and yield (Zhao et al. 2004; Adame et al. 2008; Cao 2022), and these models can provide detailed information about individual trees in a forest stand, support the simulation of different forest structures, and provide flexible output to evaluate a broad range of stand treatments with a higher prediction resolution than that of stand-level and size class-level models (Peng 2000; Cao 2014). However, the traditional growth model includes only some biological factors, such as tree size and stand conditions, which cannot predict growth under changing environmental conditions. The growth of trees is affected by a variety of biotic and abiotic factors (Toledo et al. 2011; Ford et al. 2017). Therefore, only a comprehensive assessment of the impact of various factors can accurately quantify the growth of trees, especially in the context of global climate change, and it is particularly important to develop growth models including biotic and abiotic factors.

Various categories of traditional variables, which are often called biotic variables and further classified as individual-level, stand-level, and competition-level variables, have been added to models (Adame et al. 2008; Pretzsch and Schütze 2009; Zhang et al. 2017). Individual-level variables have been widely used as substitutes for tree size or tree vitality. Stand-level variables have been included to reflect the differences between different stands. Competition-level variables that are related to stand density can be divided into distance-independent (nonspatial) and distance-dependent (spatial) variables (Tomé and Burkhart 1989). Distance-independent variables are more commonly used because they are easier to incorporate into models and are unlikely to provide a prediction accuracy that is much different from that obtained using distance-dependent variables (Dong et al. 2021). An increasing number of studies have revealed that traditional model variables have limitations when dealing with different forest regions and climatic regions, and it has been proposed that multiple variables related to environmental factors, such as site-growth variables, i.e., diversity, topography and climate variables, should be added to individual-tree diameter growth models (Rosenberg et al. 1983; Laubhann et al. 2009; Haase et al. 2015). Diversity reflects the mixed proportion of tree species or the richness of species, and diversity is one key difference between pure forests and mixed forests (Yang et al. 2009; del Río et al. 2019). It is also necessary to consider the abiotic factors related to tree growth. For example, aspect, slope, and elevation are recognized as important topographic variables that affect the amount and daily cycle of solar radiation received at different times of the year (Stage 1976; Fekedulegn et al. 2003). Climate variables linked to temperature (dormant and growing season), precipitation patterns, and drought intensity, which modify species growth patterns and result in productivity shifts among species, are critical to predicting tree- and stand-level growth in the context of climate change (Condés and García-Robredo 2012; Begović et al. 2020; Aldea et al. 2021). Thus, forest growth models used to inform adaptive management strategies should incorporate the sensitivity of forest dynamics to climate change (Pukkala and Kellomaki 2012; Zhang et al. 2016; Yasmeen et al. 2019).

Ordinary least squares (OLS) was the first type of model used to model individual-tree diameter growth (Newnham and Smith 1964; Hasenauer et al. 1998; Kiernan et al. 2008). However, this approach has significant deficiencies when dealing with stands with multiple tree species and age groups or with hierarchically structured data (Kuehne et al. 2016). Linear mixed-effects models have been widely used to analyze repeated sampling survey data with a nested structure (Lhotka and Loewenstein 2011; de-Miguel et al. 2013; Sanchez-Gonzalez et al. 2021). The linear mixed-effects model is composed of fixed- and random-effect parameters, and the variance-covariance structure allows for the efficient analysis of hierarchically structured data and increases the prediction accuracy of the model when the measured trees or stands are grouped into plots or regions (Ugrinowitsch et al. 2004). In addition, linear mixed-effects models may be calibrated to improve predictive ability if the values of the random parameters are predicted based on a subsample of trees measured in particular stands (Miao et al. 2021). However, the prediction accuracy may be lower than the same model fitted by OLS if samples were unavailable for random effects calibration (Xie et al. 2021).

Boreal forests are widely distributed in the eastern Daxing’an Mountains in Northeast China. These boreal forests play a crucial role in the Chinese national carbon budget and climatic system. The eastern Daxing’an Mountains have provided and continue to provide wood resources for the development of China (Ho 2006). From the 1960s to 2000s, the boreal forests in these mountains were overharvested in the absence of appropriate scientific management; thus, many overharvested secondary forests now exist in this region, and their growth is slow due to poor site conditions and frequent disturbances, including fires and snowstorm disasters (Guo et al. 2017; Chave et al. 2020; Zhu and Lo 2022). In 1998, the natural forest protection program (NFPP) was implemented in the eastern Daxing’an Mountains in Northeast China to control deforestation and forest degradation and protect upstream forest ecosystems and watersheds. After more than 20 years of protection, tree growth models that are suitable for this region are needed to estimate the developmental trends of future stands and the potential productivity and supply capacity of the land. Dahurian larch (Larix gmelinii Rupr.) and white birch (Betula platyphylla Suk.) are the two dominant species in the natural forests of the eastern Daxing’an Mountains (State Forestry and Grassland Administration 2019). Unfortunately, the forest resource inventory data in this region have not been fully developed, and no suitable growth models have been available to date for Dahurian larch and white birch in the Daxing’an Mountains.

Therefore, this study aimed to analyze and explore the diameter growth responses to biotic and abiotic factors for Dahurian larch and white birch in secondary forests of Northeast China. The detailed objectives were to (1) develop individual-tree mixed-effects models of diameter growth for Dahurian larch and white birch in a boreal overharvested secondary forest across the eastern Daxing’an Mountains in Northeast China; (2) evaluate the performances of a linear mixed-effects model and a linear model, compare them with a leave-one-out cross-validation approach and determine an appropriate sample size that considers both sampling cost and predictive accuracy; and (3) analyze the species-specific diameter growth response to different conditions, including different species compositions and climate conditions, for a better understanding of the growth patterns of the two species.

2 Materials and methods

2.1 Study area

The study area is located in the eastern Daxing’an Mountains (from 121° 12′ E to 127° 00′ E and 50° 10′ N to 53° 33′ N) (Fig. 1A). The elevation of the study area ranges from 190 to 1190 m above sea level. This area has a cold temperate continental monsoon climate, with the monthly mean annual temperature and precipitation shown in Fig. 1B and C (Zhao et al. 2016). More than 64% of the precipitation in this area falls as rain, primarily in summer (July–August). Snow typically starts in October and melts in April, lasting for an average of 6 months. The frost-free period is 85–105 days, usually spanning from May to September. The spring thaw period lasts for approximately 70 days, starting when the daily maximum air temperature exceeds 0 °C and ending with the soil thawing to a depth of 20 cm (Bai et al. 2019). The soils are mostly podzol according to the soil groups in the classification system of the Food and Agriculture Organization (Sposito 2023). The forest types in the eastern Daxing’an Mountains are mainly divided as follows: white birch forest, larch forest, larch-birch forest, deciduous broadleaf mixed forest, coniferous and broadleaf mixed forest, and coniferous mixed forest. Among these forest types, white birch forest, larch forest and larch-birch forest have the largest proportions of forest area and stand volume (State Forestry and Grassland Administration 2019). The dominant tree species are Dahurian larch and white birch; other species include Dahurian poplar (Populus davidiana Dode.), Mongolian oak (Quercus mongolica), and Mongolian pine (Pinus sylvestris var. mongolica).

Fig. 1

The geographical location of the study area and the plot distribution in the eastern Daxing’an Mountains, China (A); monthly average temperature for the climate period (2005–2015), with elevation ranging from 190 to 1190 m, in the Daxing’an Mountains (from 121° 12′ E to 127° 00′ E and from 50° 10′ N to 53° 33′ N) (B); monthly average precipitation for the climate period (2005–2015), with elevation ranging from 190 to 1190 m, in the Daxing’an Mountains (from 121° 12′ E to 127° 00′ E and 50° 10′ N to 53° 33′ N) (C). The red dots and numbers represent the forest bureaus and the number of plots nested in each forest bureau, respectively

2.2 Forest inventory data

The data in this study were collected from permanent sample plots managed by twelve forest bureaus in China’s National Forest Continuous Inventory (NFCI) in the eastern Daxing’an Mountains; the numbers of permanent sample plots for birch forest, larch forest, and larch-birch forest were 168, 212, and 240, respectively. The permanent plots were investigated in 2005–2015 at 5-year intervals. The permanent sample plots were located throughout the species’ distributions across the eastern Daxing’an Mountains; they were established at the intersections of a grid-based network (8 × 8 km) and were square, with an area of 600 ${m}^{2}$. The species of all the sample trees in each plot were identified, and the diameter at breast height (DBH) was recorded for each tree when it was greater than 5 cm. The mean height of three to five medium-sized trees, selected according to the average DBH, of each dominant species was taken as the average height of the species. The summary information of the data used in this study is listed in Table 1, and the abbreviations used in the NFCI and the calculation formulas of the variables are shown in Table 2.

Table 1 Descriptive statistics of the tree- and plot-level variables for larch and birch species in 2005–2015. The definitions of the variable abbreviations and the calculation formulas are presented in Table 2

Table 2 Abbreviations, descriptions and formulas of available variables for individual-tree diameter growth models

2.3 Climate data

To explore the effects of climate on the diameter growth of Dahurian larch and white birch, we obtained monthly and annual climate data for each plot for the entire survey interval (2005–2015) using ClimateAP v2.30 software (Wang et al. 2017). ClimateAP is a standalone MS Windows^® software application used for extracting climate data in the Asia–Pacific area, and it can provide downscaled gridded (4 × 4 km) monthly climate data for both past and future years or periods between 1950 and 2100 according to the latitude, longitude, and elevation of the plots. In our study area, the elevation ranged from 190 to 1190 m, the longitude ranged from 121° 12′ E to 127° 00′ E, and the latitude ranged from 50° 10′ N to 53° 33′ N. Climate downscaling was achieved through a combination of bilinear interpolation and dynamic local elevational adjustment (Marke et al. 2011). In this study, two growth season variables, twelve seasonal variables, and fifteen annual variables were chosen as candidate climate variables for model fitting. The input values of the climate variables were the mean values over each survey interval (5 years). The basic information on climate variables is listed in Table 3, and the corresponding climate variable names, abbreviations, and their definitions are shown in Table 4.

Table 3 Summary information of the available climate variables extracted from the ClimateAP software. The definitions of the variable abbreviations are presented in Table 4

Table 4 List of abbreviations and explanation of available climate variables for growth models. The employed values of the climate variables were the mean values over each survey interval

2.4 Model development

2.4.1 Basic model development

The individual-tree diameter growth model is a function that describes the relationship between diameter increase and controlled factors. The general form of the basic model is expressed as follows:

$$y_i=x_i^{'}\beta+\varepsilon_i$$

(1)

where ${y}_{i}$ is the ith observation of response variable, $\beta$ is the coefficients into vector, ${x}_{i}$ is the predictors of observation i into vector, ${\varepsilon }_{i}$ is the error term. In this study, both the increment of diameter and the increment of basal area were considered as response variables since these two forms have a mathematical relationship and are convenient for use in forestry (Hökkä and Groot 1999). And the final form of the response variable was determined by evaluating the fitting performance of the model that includes the different response variable, namely we used the type of response variable for which the model fit was better. The alternative predictors are generally divided into two types: (1) biotic factors (such as tree size, competition effects, stand attributes, and diversity index) and (2) abiotic factors (such as topographic conditions and climate factors), and the best combination of predictors was selected through a method of best-subset regression (Hofmann et al. 2020). To show the importance of climate factors in simulating the diameter growth, we will attempt to develop growth models with/without climate variables for both larch and birch species. In the process of selecting the predictors, a preliminary correlation analysis between the response variable and predictors was used to remove predictors with Pearson correlation coefficients below 0.1, because it would require considerable computation resources if all alternative predictors were used with a best-subset regression. In addition, significance and collinearity were considered during the process of selecting the appropriate predictors.

In many cases, both the conjoint contribution and its independent contribution of all predictors are studied. To identify the relative importance of the selected predictors for the individual-tree diameter growth model, hierarchical partitioning analysis was applied (Chevan and Sutherland 1991). In this study, the hier.part package in R software was used to perform hierarchical partitioning analysis (Nally and Walsh 2004).

2.4.2 Linear mixed-effects model

The NFCI data used in this study were typical longitudinal data from plots managed by several forestry bureaus. By identifying the nested and hierarchical structure of the data, the mixed-effects model has been widely used to explain the difference of each level (Fu et al. 2013). In this study, the measurement unit was individual trees, which were distributed in different plots nested within forestry bureaus; therefore, two levels (plot and forestry bureaus) and one level (plot or forestry) were evaluated. The general form of a two-level linear mixed-effects model can be written as (Pinheiro and Bates 2000):

$$\begin{array}{c}{y}_{ij}={X}_{ij}\beta +{Z}_{i,j}{b}_{i}+{Z}_{ij}{b}_{ij}+{\varepsilon }_{ij}\\ i=\mathrm{1,2},3,\dots ,M, j=\mathrm{1,2},3,\dots ,{M}_{i}\\ \begin{array}{c}{b}_{i}\sim N(0,{\psi }_{1})\\ {b}_{ij}\sim N(0,{\psi }_{2})\\ {\varepsilon }_{ij}\sim N(0,R)\end{array}\end{array}$$

(2)

where $i$ and $j$ represented the level of forestry bureaus and plot, and $M$ was the number of forestry bureaus and ${M}_{i}$ was the number of plots in the $i$ th forestry bureaus. ${y}_{ij}$ was response vector of ${n}_{ij}\times 1$ and n_ij was the number of trees in jth plots of ith forestry bureaus, ${X}_{ij}$ was fixed-effects model matrices of size ${n}_{ij}\times p$ where $p$ was the number of predictors, $\beta$ was the fixed-effects parameters. ${b}_{i}$ was the forestry bureaus-level random effects of length ${q}_{1}$ and ${b}_{ij}$ was the plot-level random effects of length ${q}_{2}$. ${Z}_{i,j}$ and ${Z}_{ij}$ were corresponding model matrices of size ${n}_{ij}\times {q}_{1}$ and ${n}_{ij}\times {q}_{1}$. ${\psi }_{1}$ and ${\psi }_{2}$ were the variance–covariance matrices of ${b}_{i}$ and ${b}_{ij}$, which were defaulted to the positive definite matrix structure in this study. ${\varepsilon }_{ij}$ were the within-group errors and $R$ was the corresponding the variance–covariance matrix. It should be noted that all plots used in our study were repeatedly measured three times (in 2005, 2010 and 2015); that is, for one tree, there was a time correlation between the increment of two 5-year growth periods (2005–2010, 2010–2015). This issue could also be solved through the mixed-effects model, specifically through the customized structure of $R$ for the error terms, which is expressed as follows:

$$R={\sigma }^{2}{G}^{0.5}\Gamma {G}^{0.5}$$

(3)

where ${\sigma }^{2}$ is a scaling factor of the error dispersion, equal to the residual variance in the estimated model; $G$ is a diagonal matrix that describes the heteroscedasticity variances of the within-tree error; and $\Gamma$ is a correlation structure matrix of the within-tree error. In this study, the autoregressive correlation structure of order 1 (namely, AR(1)) was applied for the correlation structures, and this structure aimed to simulate the temporal correlation between the increment of the two 5-year growth periods.

The one-level linear mixed-effects model could be expressed easily since it followed the same general pattern of Eq. (2). The maximum likelihood method was used when comparing mixed-effects models with different levels and random effects, and then the Akaike information criterion (AIC) and log likelihood (LL) methods were used to compare the fitting performance of mixed-effects models with different random effect parameters. In addition, the likelihood ratio test was used to compare the difference between one-level and two-level mixed-effects models. Finally, the restricted maximum likelihood method was used to obtain the parameter estimates of the selected mixed-effects models. All the computation of model fitting and comparison metrics was done using the nlme package in R 4.2.0 software (R Core Team 2020).

2.4.3 Calibration for model prediction

It is common to use developed models to make direct predictions with data that are outside the scope of the modeling data, but the resulting accuracy is often unsatisfactory (e.g., Miao et al. 2021; Hao et al. 2022). Therefore, a calibration for the model predictions was applied to both basic and mixed-effects models in this study. Different calibration strategies were adapted for basic and mixed-effects models, which referred to the research by Temesgen et al. (2008). If there was a subsample that includes the diameter increment information for a previous period, the model-specific predicted values could be calibrated as follows:

• (a) for the basic model:

A subject-specific calibrated coefficient (${c}^{*}$) was calculated based on Eq. (4), and the subject-specific predicted values could be obtained with Eq. (5):

$${c}^{*}=\frac{{\sum }_{i=1}^{m}{y}_{i}\cdot \left(X\widehat{\beta }\right)}{{\sum }_{i=1}^{m}{\left(X\widehat{\beta }\right)}^{2}}$$

(4)

$${{\widehat{y}}_{i}}^{*}={c}^{*}\cdot \left(X\widehat{\beta }\right)$$

(5)

where $m$ is the size of the subsample, ${y}_{i}$ is the observed value of the response variable in the subsample, $X$ is the design matrix, and $\widehat{\beta }$ is the parameter estimate in the basic model.

• (b) for the mixed-effects model:

Different from the basic model, the subject-specific predicted values of the mixed-effects model were obtained through both fixed effect and random effect parts, where the random effect part had to be calibrated with a subsample (Yang et al. 2009). We used the estimated best linear predictors (EBLUPs) method to estimate the random effect for a specific subject (Mehtatalo and Lappi 2020):

$${\widehat{u}}_{i}=\widehat{D}{Z}_{i}^{T}({Z}_{i}\widehat{D}{Z}_{i}^{T}+{\widehat{R}}_{i}{)}^{-1}{\widehat{e}}_{i}$$

(6)

where ${\widehat{u}}_{i}$ is the estimated random effect parameter vector; $\widehat{D}$ is the variance-covariance matrix of random effect parameters, which was replaced by the estimates obtained from model fitting; ${Z}_{i}$ is the matrix of the partial derivatives of the model function corresponding to the random effect parameters; ${\widehat{R}}_{i}$ is the variance-covariance matrix of the error term, which can be recalculated according to Eq. (3); and ${\widehat{e}}_{i}$ is the bias, which is defined as the difference between the observed increment of diameter growth and the predicted increment of diameter growth by the fixed effects parameters.

Without local observations of response variables, namely, the null subsample, the direct predicted values without calibration for basic models and the fixed part of the mixed-effects model in prediction were both used as controls to compare the effects of calibration on prediction.

2.4.4 Assessment of model fitting and validation

The adjusted coefficient of determination (${R}_{a}^{2}$, Eq. (7)) and root mean square error ($RMSE$, Eq. (8)) were used to assess the fitting performance of each model developed in this study. ${R}_{a}^{2}$ reflects the degree to which the selected predictors explain the variation in the response variable, and RMSE represents the fitted bias; therefore, a larger ${R}_{a}^{2}$ and a smaller RMSE indicate a better corresponding model.

$${R}_{a}^{2}=1-\left(\frac{n-1}{n-p}\right)\frac{{\sum }_{i=1}^{n}{\left({y}_{i}-{\widehat{y}}_{i}\right)}^{2}}{{\sum }_{i=1}^{n}{\left({y}_{i}-{\overline{y} }_{i}\right)}^{2}}$$

(7)

$$RMSE=\sqrt{\sum\limits_{i=1}^{n}{\left({y}_{i}-{\widehat{y}}_{i}\right)}^{2}/\left(n-p\right)}$$

(8)

Generally, model validation with independent data is necessary and important to evaluate the applicability and reliability of the developed model. The most difficult aspect of independent validation is the acquisition of independent data. At present, the common practice is to use the leave-one-out cross-validation (LOOCV) method, which contributes to optimizing the investigated data for model construction and can be used to simulate subsampling data (Kearns and Ron 1997). After obtaining the predicted values of all the developed models based on the description above, four statistics—mean error ($ME$), mean absolute error ($MAE$), mean percentage error ($MPE(\mathrm{\%})$), and mean absolute percentage error ($MAPE\left(\mathrm{\%}\right)$)—were calculated with the model bias generated from the validation process to assess and compare the predictive performances of the models. The formulas were as follows:

$$ME=\sum_{i=1}^{n}\frac{{y}_{o,i}-{\widehat{y}}_{p,i}}{n}$$

(9)

$$MAE=\sum_{i=1}^{n}\left|\frac{{y}_{o,i}-{\widehat{y}}_{p,i}}{n}\right|$$

(10)

$$MPE(\mathrm{\%})=\frac{{\sum }_{i=1}^{n}\left({y}_{o,i}-{\widehat{y}}_{p,i}\right)/{\overline{y} }_{o}}{n}\times 100$$

(11)

$$MAPE(\mathrm{\%})=\frac{{\sum }_{i=1}^{n}\left|{y}_{o,i}-{\widehat{y}}_{p,i}\right|/{y}_{o,i}}{n}\times 100$$

(12)

where ${y}_{o, i}$ is the diameter observation of the $i$ th tree, ${\widehat{y}}_{p,i}$ is the $i$ th predicted value of the developed models obtained by the LOOCV method, ${\overline{y} }_{o}$ is the mean value of all the observations, and $n$ is the number of all trees. Please note that the four statistics presented in Eqs. (9) ~ (12) were calculated using observations and model predictions in the original scale, specifically the observed and predicted diameter transformed from the model predictions. This was done in order to directly evaluate the magnitude of the bias between the measured and predicted diameter.

Available evidence has shown that prediction accuracy can be improved when the model predictions are calibrated by a subsample (Temesgen et al. 2008). In addition, the subsample size can influence the calibrated effect in model prediction for the same model (Hao et al. 2022). In this sense, subsamples with different sample sizes in a plot, which were obtained by random sampling without replacement, were used for prediction calibration. Specifically, various l (l = 1, 2, 3, …, 15) sample sizes were evaluated through bias statistics. In addition, the process of acquiring the subsample was repeated 500 times in each size to avoid random bias. In this analysis, we expected to find an appropriate subsample size that comprehensively considered the model prediction accuracy and measurement cost.

3 Results

3.1 Basic model

When the predictors were kept the same, the fitting performance of the corresponding model was better when the increment of basal area was used as the response variable. Therefore, the logarithm of 5-year growth in squared diameter ($ln\left({DBH}_{2}^{2}-{DBH}_{1}^{2}+1\right)$) was chosen as the final response variable. However, the final predictors used in the growth model for Dahurian larch and white birch were determined by collinearity and significance analyses according to the results of full subset regression. The final forms of the basic model were as follows (GMwoc: growth model without climate variables; GMwc: growth model with climate variables):

Dahurian larch:

$$\begin{array}{cc}\mathrm{GMwoc}:& ln\left({DBH}_{2}^{2}-{DBH}_{1}^{2}+1\right)={\beta }_{0}+{\beta }_{1}{DBH}_{N}+{\beta }_{2}{BALD}_{1}+{\beta }_{3}{BALD}_{2}+{\beta }_{4}lnG+{\beta }_{5}ln{G}_{LR}+{\beta }_{6}ALT+\varepsilon \end{array}$$

(13)

$$\begin{array}{cc}\mathrm{GMwc}:& ln\left({DBH}_{2}^{2}-{DBH}_{1}^{2}+1\right)={\beta }_{0}+{\beta }_{1}{DBH}_{N}+{\beta }_{2}{BALD}_{1}+{\beta }_{3}{BALD}_{2}+{\beta }_{4}lnG+{\beta }_{5}{lnG}_{LR}+{\beta }_{6}ALT+{\beta }_{7}TD+{\beta }_{8}PPT\_SON+\varepsilon \end{array}$$

(14)

White birch:

$$\begin{array}{cc}\mathrm{GMwoc}:& ln\left({DBH}_{2}^{2}-{DBH}_{1}^{2}+1\right)={\beta }_{0}+{\beta }_{1}{DBH}_{L}+{\beta }_{2}{BALD}_{1}+{\beta }_{3}lnG+{\beta }_{4}{G}_{BR}+{\beta }_{5}ALT+\varepsilon \end{array}$$

(15)

$$\begin{array}{cc}\mathrm{GMwc}:& ln\left({DBH}_{2}^{2}-{DBH}_{1}^{2}+1\right)={\beta }_{0}+{\beta }_{1}{DBH}_{L}+{\beta }_{2}{BALD}_{1}+{\beta }_{3}lnG+{\beta }_{4}{G}_{BR}+{\beta }_{5}ALT+{\beta }_{6}Tmax\_MAM+{\beta }_{7}PPT\_DJF+\varepsilon \end{array}$$

(16)

where ${\beta }_{0}$~${\beta }_{8}$ are model parameters, ${DBH}_{1}$ and ${DBH}_{2}$ are the initial and final tree diameters in the period, ${DBH}_{N}$ is the reciprocal of the initial tree diameter, ${DBH}_{L}$ is the logarithmic initial tree diameter, ${BALD}_{1}$ is the ratio of the intraspecific basal area of larger trees to ${DBH}_{1}$, ${BALD}_{2}$ is the ratio of the interspecific basal area of larger trees to ${DBH}_{1}$, $G$ is the logarithmic stand basal area, ${G}_{LR}$ is the basal area proportion of larch, ${G}_{BR}$ is the basal area proportion of birch, $ALT$ is the elevation, $TD$ is the difference between the mean warmest month temperature and mean coldest month temperature, Tmax_MAM is the maximum temperature for spring, and PPT_SON and PPT_DJF are the mean precipitation for autumn and winter, respectively. In addition, the response variable was increased by 1 to avoid the same initial and final DBH in the period. The parameter estimates for Eqs. (13)–(16) are presented in Table 5, and they were all significant (p < 0.05). In addition, the fitting performances of GMwoc and GMwc are provided in Table 5, where GMwc performed better than GMwoc; specifically, GMwc had a higher ${R}_{a}^{2}$ (0.3594 for larch and 0.4040 for birch) and a lower RMSE (0.9527 for larch and 0.9390 for birch), which indicated that the introduction of climate variables improved the model.

Table 5 Parameter estimates and fitting indices of the best basic models for the larch and birch species

No systematic trend was observed when residuals were plotted against the corresponding model fitted values (Fig. 2), which suggested that the diameter increment model was able to adequately describe the variabilities of response variables. Figure 3 shows the proportion of variation explained by the respective predictors to the response variable based on the hierarchical partitioning analysis of GMwc. Intraspecific competition contributed the most in GMwc, followed by initial DBH, interspecific competition, elevation, basal area, basal area proportion of larch, temperature, and precipitation for larch, and we found that the proportion of variation explained by competition and tree size was more than 70% for larch. For birch, both intraspecific competition and tree size contributed the most (more than 70%) to the corresponding GMwc, followed by basal area, temperature, basal area proportion of birch, elevation, and precipitation. Therefore, topography and climate had little effect on diameter growth, although the corresponding variables were statistically significant.

Fig. 2

Residuals plotted against fitted values of developed models in the diameter growth model

Fig. 3

Contribution ratios of predictors to the diameter increment growth model by species

The effects of all predictors except for ${DBH}_{1}$ on the diameter growth for larch and birch were simulated and analyzed by the controlled variable method (Fig. 4). In the process of simulation, mean values were used for the analyzed predictors except for the predictor of interest, which was set as the minimum, mean, and maximum values in modeling data, and therefore, the influences of different predictors on the increment of diameter growth were clearly shown. Consistent with the results of hierarchical partitioning analysis, Fig. 4 also supported that competition, especially intraspecific competition, had the largest contribution to the variation in diameter growth, followed by elevation, basal area, basal area proportion of larch, temperature, and precipitation for larch. Additionally, intraspecific competition provided the largest contribution to the diameter growth model for birch, followed by basal area, temperature, basal area proportion of birch, elevation, and precipitation. For both larch and birch, the simulation shown in Fig. 4 indicated that the diameter increment of a 5-year period decreased with the increasing competition effect described by ${BALD}_{1}$, ${BALD}_{2}$, and $LnG$. The diameter increment of a 5-year period also decreased as the elevation and precipitation of autumn and winter increased. Interestingly, the effect of stand composition described by the species-specific basal area proportion on the 5-year diameter increment of the corresponding species was completely opposite, namely, the basal area proportion of larch had a negative effect on the 5-year diameter increment of larch, and the basal area proportion of birch had a positive effect on the 5-year diameter increment of birch.

Fig. 4

Responses of the DBH increment curve to changes in different covariates for larch and birch trees. A, B, and C represent cases in which the control variables are at the minimum, median, and maximum, respectively

3.2 Mixed-effects models

The stepwise addition of the random effects of different levels (plot and forestry bureaus) to Eqs. (13)–(16) resulted in the best combination of random effects, according to the AIC and the likelihood ratio test (LRT). The LRT result supported that the performance of the plot-level mixed-effects model was better than that of the forestry bureau-level mixed-effects model and was not significantly different from that of the two-level (plot nested within forestry bureaus) mixed-effects model. Therefore, the plot-level mixed-effects diameter growth models were developed based on GMwoc and GMwc, namely, the mixed-effect model without climate variables (MGMwoc) and with climate variables (MGMwc). The final formulas were as follows:

Dahurian larch:

$$\begin{array}{cc}\mathrm{MGMwoc}:& ln\left({DBH}_{2}^{2}-{DBH}_{1}^{2}+1\right)=({\beta }_{0}+{u}_{0})+({\beta }_{1}{+{u}_{1})DBH}_{N}+{\beta }_{2}{BALD}_{1}+{\beta }_{3}{BALD}_{2}+{(\beta }_{4}+{u}_{2})LnG+{\beta }_{5}{lnG}_{LR}+{\beta }_{6}ALT+\varepsilon \end{array}$$

(17)

$$\begin{array}{cc}\mathrm{MGMwc}:& ln\left({DBH}_{2}^{2}-{DBH}_{1}^{2}+1\right)={\beta }_{0}+{(\beta }_{1}+{u}_{0}){DBH}_{N}+{\beta }_{2}{BALD}_{1}+{\beta }_{3}{BALD}_{2}+\left({\beta }_{4}+{u}_{1}\right)+LnG({\beta }_{5}+{u}_{2}){lnG}_{LR}+{\beta }_{6}ALT+{\beta }_{7}TD+{\beta }_{8}PPT\_SON+\varepsilon \end{array}$$

(18)

White birch:

$$\begin{array}{cc}\mathrm{MEMwoc}:& ln\left({DBH}_{2}^{2}-{DBH}_{1}^{2}+1\right)=({\beta }_{0}+{u}_{0})+({\beta }_{1}+{u}_{1}){DBH}_{L}+{\beta }_{2}{BALD}_{1}+({\beta }_{3}+{u}_{2})LnG+{\beta }_{4}{G}_{BR}+{\beta }_{5}ALT+\varepsilon \end{array}$$

(19)

$$\begin{array}{cc}\mathrm{MGMwc}:& ln\left({DBH}_{2}^{2}-{DBH}_{1}^{2}+1\right)=({\beta }_{0}+{u}_{0})+({\beta }_{1}+{u}_{1}){DBH}_{L}+({\beta }_{2}+{u}_{2}){BALD}_{1}+{\beta }_{3}LnG+{\beta }_{4}{G}_{BR}+{\beta }_{5}ALT+{\beta }_{6}Tmax\_MAM+{\beta }_{7}PPT\_DJF+\varepsilon \end{array}$$

(20)

where ${u}_{0}$, ${u}_{1}$, and ${u}_{2}$ are random effect parameters, and the others are defined above. The fitting results of the mixed-effects models are listed in Table 5, where the parameter estimates were all significant and the signs of the parameters were consistent with the basic models. The ${R}_{a}^{2}$ of MGMwoc was 59.6 and 41.0% higher than those of GMwoc for larch and birch, respectively, as were the MGMwc and GMwc, which indicated that the introduction of the plot-level random effect improved the basic model significantly. However, the performance of MEMwoc was slightly better than that of MGMwc based on the ${R}_{a}^{2}$ and RMSE values. In addition, the correlation coefficient in the $R$ matrix (Eq. (3)) was approximately 0.4 for the two species, which suggested that there was a temporal correlation between the increments of the two 5-year growth periods.

3.3 Analysis of model validation

The bias statistics that assess the prediction accuracy of diameter for the developed models are listed in Table 6, and the results suggest that the basic and mixed-effects models we developed slightly underestimated the DBH increments for the two species since the $ME$ and $MPE$ values were greater than 0. However, they were approximately unbiased models in the prediction because $ME$ and $MPE$ were close to 0 according to Table 6. The $MAE$ s were found to be smaller than 0.4 cm, and the corresponding $MAPE$ s were also less than 4% for the two species, indicating that the basic and mixed-effects models performed adequately well for predictions related to an individual. In addition, the introduction of plot-level random effects and climate variables provided more accurate predictions, which were similar to the fitting results of the models.

Table 6 The bias statistics of diameter for the diameter increment models of larch and birch

Figure 5 further demonstrates the superiority of the MGMwc model in all diameter classes for the two species. The $MAE$ increased gradually with increasing diameter classes, but the maximum was not more than 0.7 cm. GMwoc performed worse than GMwc in each grade, which indicated that climate variables are important in modeling diameter growth. Compared with the addition of climate variables, the introduction of plot-level random effects improved the model prediction more significantly in all diameter classes.

Fig. 5

Plot of the absolute bias across different diameter classes

3.4 Relationship between prediction calibration and prediction accuracy

The prediction results obtained with the models developed under different sampling sizes are shown in Fig. 6. After calibration, GMwoc, GMwc, and MEMwc had similar prediction accuracies, which were all better than the prediction accuracy of MGMwoc. The absolute bias between the measured $\mathrm{DBH}$ and predicted $\mathrm{DBH}$ gradually decreased as the number of samples increased; however, the rate of decrease tended to be stable when the sample size was more than nine trees. In Fig. 6, the prediction calibration actually increased the prediction bias at the start of sampling (sample size was smaller than 5), and then, calibration gradually improved the prediction effect of the basic model. For the mixed-effects model, the prediction accuracy of the calibrated models was still improved compared with that of the uncalibrated models and was better than that of the basic model under the same conditions, even if only one tree was used in model calibration. In addition, MGMwoc performed slightly worse than MGMwc in prediction calibration with small samples; otherwise, the opposite was true. Overall, the prediction errors tended to be basically stable when the sample size was greater than nine trees, and to limit costs while achieving relatively high efficiency, we suggest sampling nine trees per plot for the prediction of individual-tree diameter growth in birch and larch with the developed models in this study.

Fig. 6

Relationship between prediction biases and the sample sizes used to calibrate the growth models of larch and birch. The red lines represent the prediction bias in the scenario in which the calibration data are not available

4 Discussion

4.1 Responses to biotic and abiotic factors

Individual-tree diameter growth models are frequently estimated using tree-, competition-, site-, and climate-level variables, which are usually easy to obtain via field investigations and software such as ClimateAP (Biging and Dobbertin 1995; Fekedulegn et al. 2003; Kuehne et al. 2016; Saud et al. 2019). Based on the fitting results of GMwc, the variable ${\mathrm{BALD}}_{1}$ was the most important tree-level variable, which referred specifically to intraspecific competition and indicated the social status of individual conspecific trees; it was calculated as $\mathrm{BAL}$ divided by the logarithmic diameter of the tree. The negative correlation between ${\mathrm{BALD}}_{1}$ and tree growth showed that competition within tree species reduced tree growth rates. This finding was supported by the findings of Wykoff et al. (1982); since there was less light available for smaller trees, ${\mathrm{BALD}}_{1}$ may be a surrogate for light measurements. ${\mathrm{BALD}}_{2}$ was also considered to be an effective evaluation index for competition among different species, which influenced the diameter growth of larch in our study. The coefficient of ${\mathrm{BALD}}_{2}$ was slightly greater than that of ${\mathrm{BALD}}_{1}$, as shown in Table 5, indicating that the growth of larch was more sensitive to changes in interspecific competition. Therefore, mixing with birch can promote the diameter growth of larch in a scenario in which larch has an advantage in interspecific competition. In addition, tree size was an important tree-level variable, which may reflect the growth time in mixed forests, and the positive value of this parameter indicated that the diameter growth accelerated with increasing ${\text{DBH}}$ in the initial period. $\mathrm{lnG}$ represents the stand attributes, especially when substituted for the stand density variable, and $\mathrm{G}$ is the stand basal area of all species and was negatively correlated with growth. Increased crowding and reduced growing space resulted in lower growth rates.

The variables ${G}_{LR}$ and ${G}_{BR}$ represent the basal area proportions of larch and birch species in the stand and were used in this study to represent species diversity or mixture (Riofrío et al. 2019; Bottero et al. 2021). The parameter estimate for ${G}_{LR}$ was negative and had the opposite sign as that for ${G}_{BR}$. This difference showed that the growth rate of the birch species was suppressed, while the growth of the larch species was promoted when the two species were mixed, which was the same as the analysis of the effect of competition. Figure 4 also shows the differences in growth under different proportions of these two species, with larch having a more obvious gain. The larch species, a more important managed tree species, showed negative growth responses to intraspecific neighbors, but the effects were counterbalanced at the stand level by a corresponding increase in interspecific trees. The principles of interaction among mixtures of species are complex; they involve nutrient cycling, photosynthesis, and soil physical and chemical properties, and they remain mostly unknown. Larch and birch have similar growth traits, such as light demands and preferences for fertile soil, and they can survive in regions with severe environments, such as high mountain sites or along the upper parts of valleys (Kloeppel et al. 1998). During interspecific interactions, larch competes effectively with birch due to its higher resistance to harsh climatic conditions. Larch surpasses other tree species in water use efficiency and can survive at a semidesert level of precipitation (Berg and Chapin 1994). Moreover, larch can survive in the permafrost zone because its deciduous leaf characteristic and dense bark protect stems from winter desiccation and snow abrasion (Mao et al. 2010).

The parameter estimates of variable $\mathrm{ALT}$ were negative for both larch and birch species, showing that tree growth was restrained with increasing elevation. This pattern was consistent with the results of a study on stand-level volume growth (Wang et al. 2021). In the process of adding the climate variables, the best-performing variables in this study were $\mathrm{TD}$ and $\mathrm{PPT\_SON}$ for larch and $\mathrm{Tmax\_MAM}$ and $\mathrm{PPT\_DJF}$ for birch. The effects of precipitation and temperature on diameter growth were not uniform for the two species, namely, precipitation inhibited the diameter growth of larch and birch trees, while temperature had a negative effect on the diameter growth of larch and a positive effect on the diameter growth of birch. However, both larch and birch species were more sensitive to the climate variables related to temperature than to those related to precipitation (Fig. 4). TD was calculated as the difference between the warmest and coldest months, which affected the aboveground conditions and soil conditions (Zhang et al. 2019). In general, the lower the maximum temperature in the coldest month in winter is, the higher the ${\text{TD}}$, which strongly affects winter ground conditions. Larch rarely grows in winter because its leaves are stripped, and most physiological activity tends to stall; however, the higher temperature in the warmest month promotes the growth of larch. The maximum temperature in spring was more important than that in the other seasons for birch, and it may play an important catalytic role in photosynthesis for this species (Levani and Eggertsson 2008). In addition, a higher maximum temperature in spring may lead to a reduced frozen depth, which allows for greater tree root expansion and higher root activity in the future (Alvarez‐Uria and Körner 2007). The precipitation in autumn and winter had a negative effect on the growth of birch and larch, respectively. The region of this study in China with the highest latitude and the temperature in autumn and winter is less than 0℃ (Fig. 1) which indicated that the precipitation (snow) in cold seasons would slow the diameter growth of larch and birch. The same result was presented in studies of the common juniper (Juniperus communis L.) growth in the Alpine tundra (Carrer et al. 2019). Some other studies had opposite results that snow promote the growth by providing protection against frosts and melting snow as water source for next growth season (Rixen et al. 2010; Hallinger et al. 2010). However, here the snowpack could shorten cambial activity by delaying the onset or anticipating the end of the growing season the extensive duration of snow cover delays the onset of the vegetative period, inducing a negative effect on growth (Francon et al. 2017). And the negative effect is likely not direct and it might represent the collateral consequence of temperature on snow characteristics (Carrer et al. 2019).

4.2 Model calibration

In the mixed-effects model, the fixed-effects parameters were global to all groups, whereas the random effects parameters were plot specific. One new observation was that the mixed-effects model should be calibrated, which involves predicting the random effect parameters for the group, i.e., individual-tree diameter growth. New samples were used to predict the random effect parameters by estimating the best linear predictors (Mehtatalo and Lappi 2020), and the effects of different sample sizes on the prediction accuracy for the base models (GMwoc and GMwc) and mixed-effects models (MGMwoc and MGMwc) were evaluated. To ensure the randomness of sampling, the sampling of random trees was conducted to simulate new observations to calculate the random effects parameter to calibrate the models. Figure 6 shows that the model performance improved as the number of trees sampled increased, which was supported by other studies (Bronisz and Mehtatalo 2020a, b; Miao et al. 2021). However, the prediction accuracy tended to be basically stable when the sample size was greater than nine trees. When the sample size was small, the calibrated base model was less accurate than the uncorrected model, and the prediction accuracy of the mixed-effects model after calibration was always higher than that of the uncalibrated scenarios. Therefore, the prediction accuracy of the fixed effects (the scenario in which the calibration data are not available) for the mixed-effects model was lower than that of the basic model if the premeasured samples could not be obtained. In this case, it is not recommended to use the mixed-effects model to predict the future growth of trees. The mixed-effects model was more suitable for predicting the growth of trees at the present stage since it was easy to calibrate the samples. Because the measurement of DBH is relatively simple, nine trees were randomly selected for each plot for the calibration of random effects. Compared with the scale of the whole forest, this approach can not only lower the cost of investigation but also ensure the prediction accuracy of the mixed-effects model, which further enhances the practicability of the mixed-effect model (Bronisz and Mehtatalo 2020a, b).

4.3 Prospect and notes

In China, nine national forest resources continuous inventories have been completed, and a large amount of dynamic forest resource data has been generated for both planted and natural forests. However, such a long-term and systematic NFCI dataset is typically only used for the aggregation and statistics of forest resource information, without fully exploring its deeper potential. For instance, it can be utilized for conducting large-scale research on forestry-related scientific questions or guiding the sustainable utilization and management of forest resources. Therefore, further exploration and utilization of such data are needed. Exploring the diameter growth mechanisms based on biotic and abiotic factors using these continuous and systematic NFCI data can not only help estimate the forest productivity and carbon balance on the regional scale, but also provide certain indicators for the evaluation of the structure and function of forest ecosystems. In this sense, a case study for larch and birch in Daxing’an Mountains, Northeast China, was provided here, and it has the potential to contribute not only a practical model for supporting sustainable forest management in secondary forests in the eastern Daxing’an Mountains, but also as an indicator for analyzing tree diameter growth in a wider range of forest stand types.

In the past, diameter growth models mainly considered biotic factors such as tree size and stand characteristics since these biometric factors were the main drivers linked to diameter growth (Wykoff 1990; Lessard et al. 2001; Andreassen and Tomter 2003). However, such models have certain limitations, such as they cannot evaluate the tree growth in the context of climate change. In this study, the response of diameter growth of larch and birch to tree size, stand characteristics, elevation, precipitation, and temperature were simulated. While obtaining general conclusions such as the positive correlation between tree size and diameter growth, and the negative correlation between competition and stand basal area with diameter growth, we also quantified the effects of elevation, temperature, and precipitation on diameter growth for larch and birch. Our study showed that larch and birch diameter growth are more sensitive to temperature than to precipitation, and less diameter growth occurs at lower temperatures and with more precipitation in winter. Furthermore, a better understanding of the growth mechanisms and patterns of trees and forests can be gained, thereby improving existing prediction models and frameworks, and more accurately predicting the future growth of trees and forests in a more scientific way.

Some biomass models for larch and birch in this region were previously developed at Northeast Forestry University, China (Dong et al. 2018). In this sense, accuracy diameter growth models would be more meaningful to contribute to research on the carbon cycle and changes in carbon dynamics in study area. Moreover, considering the influence of climate factors in diameter increment prediction could improve the accuracy of diameter estimation and further increase the prediction accuracy of the carbon storage and carbon sequestration capacities of forest ecosystems in the eastern Daxing’an Mountains, Northeast China. In addition, if our growth models are used in other regions, caution should be taken because different environmental and tree growth conditions may yield different relationships between tree growth and variables; thus, the models could have larger prediction errors.

5 Conclusion

Individual growth linear mixed-effects models were developed for Dahurian larch and white birch species in a boreal, overharvested, secondary forest in the eastern Daxing’an Mountains, China. The growth of trees increased with increasing $\mathrm{DBH}$ and decreased with increasing ${\mathrm{BALD}}_{1}$, ${\mathrm{BALD}}_{2}$, $\mathrm{lnG}$, ${\text{ln}}{\text{G}}_{\text{LR}}$, $\mathrm{ALT}$, ${\text{TD}}$, and PPT_SON for larch. For birch, tree growth also increased with increasing $\mathrm{DBH}$, ${\text{G}}_{\text{BR}}$, and Tmax_MAM and decreased with increasing ${\mathrm{BALD}}_{1}$, $\mathrm{lnG}$, $\mathrm{ALT}$, and PPT_DJF. The subject-specific predictions were obtained by calibrating the mixed-effects model using plot-level random effects based on which zero to fifteen randomly sampled trees of both species were sampled per plot. Our results indicated that (1) the mixed-effects model had obvious advantages over the OLS method; (2) to calibrate the mixed-effects models, it was best to include at least nine random additional individual tree diameter measurements to predict the growth of larch and birch; (3) larch species were dominant in interspecific competition when mixed with birch; and (4) temperature and precipitation, especially temperature, had important growth-promoting effects on both the larch and the birch species. The impact of species mixture and climate on the growth of Dahurian larch and white birch must be considered in future management policies. Quantifying the response of diameter growth to biotic and abiotic factors enables the monitoring of biological and ecological processes within relevant forest ecosystems, facilitating a more precise assessment of present and future forest resources in the area. Furthermore, it also can provide a scientific basis for formulating effective forest management strategies, thus narrowing the gap between forest ecological value and forestland management by balancing the ecological equilibrium of forests with the rational utilization of forest resources simultaneously.