Using Functional Traits to Improve Estimates of Height–Diameter Allometry in a Temperate Mixed Forest

Abstract

Accurate estimates of tree height (H) are critical for forest productivity and carbon stock assessments. Based on an extensive dataset, we developed a set of generalized mixed-effects height–DBH (H–D) models in a typical natural mixed forest in Northeastern China, adding species functional traits to the H–D base model. Functional traits encompass diverse leaf economic spectrum features as well as maximum tree height and wood density, which characterize the ability of a plant to acquire resources and resist external disturbances. Beyond this, we defined expanded variables at different levels and combined them to form a new model, which provided satisfactory estimates. The results show that functional traits can significantly affect the H–D ratio and improve estimations of allometric relationships. Generalized mixed-effects models with multilevel combinations of expanded variables could improve the prediction accuracy of tree height. There was an 82.42% improvement in the accuracy of carbon stock estimates for the studied zone using our model predictions. This study introduces commonly used functional traits into the H–D model, providing an important reference for forest growth and harvest models.

1. Introduction

In the context of global climate change and national biodiversity conservation, forests have received increasing attention due to their productivity and role as a carbon stock reservoir [1,2,3]. However, the direct estimation of forest biomass over large areas at the regional or stand scale is impractical and necessitates using refined statistical methods [4,5,6]. Allometric models have proven effective in providing accurate estimates of forest volume or biomass for carbon stocks based on extensive data measurements [7,8]. Typically, allometric models require the diameter at breast height (DBH) and tree height (H) as fundamental explanatory factors [9]. The DBH can be easily obtained through field surveys, while the tree height is more difficult to measure, and constructing accurate height–diameter (H–DBH) models will improve the accuracy of tree height predictions [5,10]. In addition, H–D models are not only used as part of estimating biomass and carbon stocks but are also utilized by foresters and ecologists to assess a tree’s access to light resources, competition with neighboring trees, or stand quality [11,12,13,14].

Despite their gradual decline, natural forests still account for 93% of the global forest area [15]. Protecting existing natural forests will continue to be the most stable and economical way to stabilize climate change in the future [16]. This implies the importance of natural forests as they contribute to forest productivity and the carbon increment [17,18]. It is necessary to obtain the tree height of natural forest stands in regional or local areas in time to understand their growth dynamics and health status [19,20]. However, compared with plantation forests, natural forests cannot easily be used to construct suitable H–D models due to their complex horizontal and vertical structures, numerous tree species, and heterogeneous ages [21,22,23]. According to previous studies, two kinds of natural mixed forest models are usually used—models for single tree species and models for whole tree species. Although these two types have been verified in many practical applications [23,24,25,26], there are still some drawbacks. Modeling for single tree species cannot be used for all tree species because of the small sample sizes of some species or the multiplicity of species in stands [9]. In contrast, modeling for a whole tree species can lead to a lack of accuracy because the characteristics of the tree species are not considered [23]. Therefore, a new idea is urgently needed to link the two together to achieve an accurate and simple H–D model.

Functional traits are important attributes of plants that respond and adapt to external biotic and abiotic resources or stresses [27,28,29], and they can express individual growth and system function both at the individual and stand levels [30,31,32]. Most studies usually use the average trait values of tree species to quantify the phenological traits of specific species, not only because average traits are easily accessible, but more importantly, the prediction of individual growth from traits measured at the species level has been shown to be reasonable [33]. Thus, in a sense, the bottom-up functional traits of plants can cover the characteristics of individual organs of tree species, such as the leaf photosynthetic capacity, defense and stress tolerance, hydraulic transport efficiency, and soil nutrient uptake capacity [34,35]. Although research on functional traits has been well established and extended to be more systematic and networked [36,37], they are not yet widely promoted in the field of forest growth and forest management. In addition, the effect of these traits on the allometric growth relationship has hardly been reported; filling this gap could provide an important reference for the trade-off between plant growth and survival.

With the continuous development of science and technology, modern H–D models could input more information to accurately serve specific forest conservation and production efforts. The generalized H–D model is an improvement from the traditional H–D model that can add more details to basic tree height and diameter relationships to further improve the prediction accuracy [38]. Jilin Province in China is rich in natural forest resources, with a natural forest area of about 6.09 million ha and a natural forest stock of 900 million m³ [39]. Therefore, the establishment of an accurate and applicable generalized mixed-effects H–D model is important for biodiversity conservation and the carbon balance in this region. In this study, we used ten tree species in Northeast China, including broadleaf and coniferous trees, with the following objectives: (I) to test whether adding functional traits to the base model improves the model’s performance; (II) to compare generalized H–D models at different levels (stand, species, and single tree) and explore the main factors affecting H–D models; (III) to combine the influencing factors at different levels and develop the best mixed-effects model for natural mixed coniferous forests in the region; and (IV) to evaluate whether using this improved model is more accurate for predicting aboveground carbon stocks in a similar region. This modeling approach will provide an important reference for growth models of natural mixed forests in the future.

2. Materials and Methods

2.1. Study Site

The study was conducted at Jingouling Experimental Forest Farm (130°5′–130°20′ E, 43°17′–43°25′ N) of the Wang Qing Forest Bureau in Jilin Province, China. The vegetation type is a mixed broadleaf/conifer forest in Changbai Mountains, and the temperate climate is characterized by alpine monsoons. The annual precipitation ranges from 600 to 700 mm, and the annual average relative humidity is 66%. The wet season is from June to August, and a relatively dry season lasts from September to May. The mean annual temperature is 3.9 °C, and the mean monthly maximum and minimum temperatures are 22 and −32 °C in July and January, respectively. The soil is brown forest soil with a rootable depth ranging from 20 to 100 cm. The dominant tree species in this area are Abies fabri, Tilia amurensis, Pinus koraiensis, Acer pictum, and Picea asperata.

2.2. Data Measurements

Data used in this study were collected in two zones. Zone I was established in 1987 with 104 permanent sample plots (20 m × 20 m), and Zone II was established with 29 plots of the same size as Zone I in 1988. In this study, Zone I was used as the training dataset to fit and compare the models in order to ensure that the best model we developed could be commonly applied to other similar types of new stands. Zone II was used as the validation dataset to evaluate the sensitivity and accuracy of the model. The distribution of Zone I and Zone II is presented in Figure 1.

All trees (DBH ≥ 5.0 cm) in the permanent sample plots were identified and measured in 2020. We only used tree species with more than 30 samples in the training dataset, as this ensures the accuracy of the modeling evaluation of the later division of the tree species. A total of 3301 pairs of height–diameter measurements were selected from 10 key species: 2602 from Zone I (training dataset) and 699 from Zone II (validation dataset). The primary information of each plot and tree species is detailed in

Table 1. The primary information of ten major tree species from the training dataset and validation dataset. D is the diameter at breast height; H is the tree height.

Species	Number	D (cm)			H (m)
		MaxD	MeanD	MinD	MaxH	MeanH	MinH
Training dataset
Abies fabri	689	42.5	14.0	5.0	26.2	10.7	2.1
Acer pictum	308	53.6	15.3	5.1	24.7	11.6	2.3
Betula costata	200	75.0	18.2	6.0	30.8	15.9	6.2
Betula platyphylla	72	39.7	21.5	8.9	26.2	17.4	8.4
Fraxinus mandshurica	46	36.9	18.9	5.3	27.8	14.9	4.1
Larix gmelinii	38	45.2	25.2	11.0	25.4	18.3	10.6
Picea asperata	308	67.1	23.7	5.0	28.4	14.1	3.2
Pinus koraiensis	315	78.5	24.2	5.0	28.4	13.9	2.0
Tilia amurensis	539	56.5	16.6	5.1	26.3	12.0	3.1
Ulmus davidiana Planch. var. japonica	87	52.0	18.9	5.2	26.1	13.2	4.6
Validation dataset
Abies fabri	64	36.5	14.4	5.1	26.8	10.0	3.0
Acer pictum	36	41.7	19.2	5.6	22.3	12.0	4.7
Betula costata	112	37.1	18.8	7.8	24.8	14.6	5.7
Betula platyphylla	6	30.7	23.1	15.9	18.7	15.3	12.4
Fraxinus mandshurica	26	38.3	22.2	6.8	22.8	14.3	5.1
Larix gmelinii	175	48.2	24.1	7.4	25.6	16.9	4.1
Picea asperata	77	50.8	26.0	5.5	28.3	15.4	3.2
Pinus koraiensis	75	62.7	22.6	6.3	24.6	12.8	4.2
Tilia amurensis	116	42.8	13.4	5.1	22.3	10.1	4.1
Ulmus davidiana Planch. var. japonica	12	56.5	26.4	8.4	20.4	14.4	7.2

All species Latin from FRPS (www.iplant.cn/frps).

and Table S2.

2.3. Functional Traits

The aim of the study was to test whether functional traits of tree species can be used as an additional variable for the H–D model. Seven functional traits that have been proven to be closely associated with plant growth and development were measured [40,41]: leaf area (LA), leaf thickness (LT), leaf dry matter content (LDMC), specific leaf area (SLA), leaf nitrogen (LN), wood density (WD), and Hmax (maximum tree height). The specific measurements of these functional traits are shown in Appendix A.

In addition to considering individual traits, we also tested the effect of combined traits on the model by using principal component analysis (PCA) to reduce the dimensionality of these seven functional traits. The PCA axis 1 (PC1) was used to represent the combined traits, explaining 54.4% of the trait variation among ten species (Figure 2).

2.4. Other Expanded Variables

The functional traits above represent expanded variables at the species level. In addition, we consider stand level and individual level variables. For stand level, we measured the elevation, slope, and aspect of each sample plot. These primary topographic variables were used to calculate two metrics that are closely related to stand structure and development [42,43]:

$$\mathrm{S}\mathrm{I}\mathrm{C}=\tan \left(\mathrm{s}\mathrm{l}\mathrm{o}\mathrm{p}\mathrm{e}\right)\times \mathrm{c}\mathrm{o}\mathrm{s}\left(\mathrm{a}\mathrm{s}\mathrm{p}\mathrm{e}\mathrm{c}\mathrm{t}\right)$$

(1)

$$\mathrm{C}\mathrm{E}=\cos \left(\mathrm{a}\mathrm{s}\mathrm{p}\mathrm{e}\mathrm{c}\mathrm{t}\right)\times \ln \left(\mathrm{e}\mathrm{l}\mathrm{e}\mathrm{v}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}\right)$$

(2)

Competition is described using the basal area (BA) of all trees in sample plots. Stand quality is an important indicator for evaluating ecosystem functioning, production, and regeneration. Here, we used the maximum height of the dominant tree species at each sample plot to quantify stand quality. At the individual level, we used basal area in large trees (BAL) of plots to target individual trees. The details of all extended variables are shown in

Table 2. Description of the expanded variables used in the generalized H–D model.

Measurement Level	Application Situation	Variable Type	Variable Name	Specific Description
Plot level	Modeling for dividing species or all trees	Topography	SIC	tan(slope) × cos(aspect)
			CE	cos(aspect) × ln(elevation)
		Competition	BA	Basal area
		Stand quality	DMH	Maximum tree height of Dominant species
Species level	Modeling for all trees	Functional trait	LA	Leaf area
			LT	Leaf thickness
			LDMC	Leaf dry matter content
			SLA	Specific leaf area
			LN	Leaf nitrogen
			WD	Wood density
			Hmax	Maximum tree height
Individual level	Modeling for dividing species or all trees	Competition	BAL	Basal area in larger trees

.

2.5. Base H–D Model Selection

We used ten nonlinear models to describe the basic allometric relationship between tree height and diameter. These empirical models have been broadly used in previous research [44,45,46] and do not have heteroscedasticity or autocorrelation problems in this study. The specific model forms and sources are presented in

Table 3. Information of H–D base models in this study. H is the tree height (m); D is the tree diameter at breast height (cm); and β₁ and β₂ are estimated model parameters.

Number	Equation	References
BM.1	$H=1.3+{\left[\frac{D}{\left({\beta }_1+{\beta }_2\times D\right)}\right]}^3$	Näslund (1936) [47]
BM.2	$H=1.3+{\beta }_1\times e^{\left(\frac{{\beta }_2}{D}\right)}$	Schumacher (1939) [48]
BM.3	$H=1.3+\frac{D^2}{{\left({\beta }_1+{\beta }_2\times D\right)}^2}$	Meyer (1940) [49]
BM.4	$H=1.3+{\beta }_1\times \left[1-e^{\left(-{\beta }_2\times D\right)}\right]$	Meyer (1940) [49]
BM.5	$H=1.3+{\beta }_1\times {\left[1-e^{\left(-{\beta }_2\times D\right)}\right]}^3$	Bertalanffy (1957) [50]
BM.6	$H=1.3+{\beta }_1{\times \left[\frac{D}{\left(1+D\right)}\right]}^{{\beta }_2}$	Curtis (1967) [51]
BM.7	$H=1.3+{\beta }_1{\times D}^{{\beta }_2}$	Stage (1975) [52]
BM.8	$H=1.3+{\beta }_1\times \frac{D}{\left({\beta }_2+D\right)}$	Bates and Watts (1980) [53]
BM.9	$H=1.3+e^{\left[{\beta }_1+\left(\frac{{\beta }_2}{D+1}\right)\right]}$	Wykoff et al. (1982) [20]
BM.10	$H=1.3+\frac{{\beta }_1\times D}{\left(D+1\right)}+{\beta }_2\times D$	Larson (1986) [54]

. Note that we only used base models with two parameters that can be easily applied in practice.

We first considered the base model for each tree species separately. This is a regular modeling approach used in many previous studies [10,25,55]. This was used as an important reference for whether our trait-based model could improve the traditional model. Next, we filtered the best base model for all trees, regardless of species. This model was used to add expanded variables to develop the generalized model.

2.6. Generalized H–D Model at Different Levels

First, to determine the effect of functional traits on allometric growth, we established the univariate linear relationship between functional traits and height–diameter ratio. Then, to include traits and information at the stand and individual tree levels, the best base model was extended to the generalized model. We used a reparameterization approach, which assumes that the expanded variables have a linear relationship with the parameters of the base model [46].

We considered the single parameter and both parameters separately to ensure a meaningful comparison of the results. The base model BM.1 is shown as an example as follows:

$$\mathrm{Base}\mathrm{model}:H=1.3+{\left[\frac{D}{\left({\beta }_1+{\beta }_2\times D\right)}\right]}^3$$

(3)

$$\mathrm{Reparameterization}\mathrm{for}{\beta }_1:H=1.3+{\left[\frac{D}{\left({\beta }_1+{\beta }_{e1}\times a+{\beta }_2\times D\right)}\right]}^3$$

(4)

$$\mathrm{Reparameterization}\mathrm{for}{\beta }_2:H=1.3+{\left[\frac{D}{\left({\beta }_1+{(\beta }_2+{\beta }_{e2}\times a)\times D\right)}\right]}^3$$

(5)

$$\mathrm{Reparameterization}\mathrm{for}{\beta }_1\mathrm{and}{\beta }_2:H=1.3+{\left[\frac{D}{\left({\beta }_1+{\beta }_{e1}\times a+{(\beta }_2+{\beta }_{e2}\times a)\times D\right)}\right]}^3$$

(6)

where $a$ denotes all expanded variables of different levels; and ${\beta }_1$, ${\beta }_2$, ${\beta }_{e1}$, and ${\beta }_{e2}$ are the estimated model parameters. Finally, we tried to combine expanded variables with different features and levels to build a generalized H–D model that contains more information.

2.7. Mixed-Effects H–D Model Building

The best-performing generalized H–D models were selected for the purpose of developing a mixed-effects model using the training dataset. All permanent plots were used as a random factor to add to the two original parameters. The following formulation was used [56]:

$$\begin{gathered} H_{ijk}=f\left({\overrightarrow{\widehat{\beta }}}_i,D_{ijk}\right)+{\epsilon }_{ijk},{\epsilon }_{ijk}~N\left(0,{\sigma }_{\epsilon }^2\right), \\ i=1,\dots ,I,j=1,\dots ,J_i,k=1,\dots ,K_{ij} \end{gathered}$$

(7)

where $H_{ijk}$ is the height of the $k$th tree from the $j$th species in $i$th plot, $f$ is the best generalized model function of the relationship between H and DBH, and ${\epsilon }_{ijk}$ is a normally distributed independent within-plots error term. ${\overrightarrow{\widehat{\beta }}}_i$ represents two parameters of the model (${\beta }_1$ and ${\beta }_2$). In this study, the $p$th element of the vector ${\overrightarrow{\widehat{\beta }}}_i$ was modeled as follows:

$${\widehat{\beta }}_{ip}={\beta }_p+{\theta }_{ip},p=\mathrm{1,2}{\theta }_{ip}~N\left(0,{\phi }_p\right),$$

(8)

where ${\widehat{\beta }}_p$ are fixed effects, and ${\theta }_{ip}$ indicates the random effect of the $i$th plot. The ${\theta }_{ip}$ follows a normal distribution with mean of 0 and variance–covariance matrix of ${\phi }_p$.

For each model, we checked the residuals for heteroskedasticity problems. The best mixed-effects model was used to estimate carbon stocks for the validation dataset. This is a key step to ensure that our models can be widely applied to any new forest stand. The specific method was as follows. First, the real biomass (tree height from our observations), biomass according to the traditional method (tree height from the base H–D model divided into tree species), and biomass according to our improved method (tree height from the best-generalized mixed-effects model) were obtained using the binary aboveground biomass (AGB) allometric growth equation based on tree height and diameter at breast height, respectively. Then, the total aboveground carbon stock of all trees was derived from the carbon coefficients of each tree species (biomass allometric models and carbon content coefficients for ten species are shown in Table S2). Finally, our model was further evaluated by comparing the difference between the two methods and the carbon stocks obtained from the actual observations.

2.8. Model Comparison

The fitting results of the base models and the generalized H–D models of the training dataset (Zone I) were compared using the root mean square error (RMSE), the mean square error for 10-fold cross-validations (CV-MSE), the Akaike information criterion (AIC), and the Bayesian information criterion (BIC).

For the RMSE, the fitting results of the models were examined by randomly splitting the training dataset again for model training (80%) and then testing the data 200 times (20%). The specific calculation process was as follows:

$$\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}=\sqrt{\frac{1}{M}\sum _{m=1}^M\frac{\sum _i^N{(H_i-{\widehat{H}}_i)}^2}{N}},$$

(9)

where RMSE is the estimated root mean square error; $H_i$ represents the observation of the dependent variable from the test data; ${\widehat{H}}_i$ indicates the fitted value of the dependent variable using the model parameters obtained from the training data to fit the test data; $M$ denotes the number of resamplings (200); and $N$ is the number of trees of a given data set per resampling $m$.

CV-MSE is a metric used to evaluate 10-fold cross-validations. The definition is as follows:

$$\mathrm{C}\mathrm{V}-\mathrm{M}\mathrm{S}\mathrm{E}=\sqrt{\frac{\sum _i^S{\left(H_i-{\widehat{H}}_i\right)}^2}{10\times S}},$$

(10)

where CV-MSE is the normalized mean square error of the cross-validation; $H_i$ represents the $i$th observation of the dependent variable from the one remaining dataset; ${\widehat{H}}_i$ represents the fitted value of the dependent variable using the model obtained from the training data to fit the test data; and $S$ is the number of every fold of data divided by 10-fold cross-validations.

AIC and BIC are metrics for evaluating the simplicity and accuracy of the model:

$$\mathrm{A}\mathrm{I}\mathrm{C}=2q-2\mathrm{l}\mathrm{o}\mathrm{g}\mathrm{L}\mathrm{i}\mathrm{k},$$

(11)

$$\mathrm{B}\mathrm{I}\mathrm{C}=\ln \left(N\right)q-2\ln \left(\mathrm{l}\mathrm{o}\mathrm{g}\mathrm{L}\mathrm{i}\mathrm{k}\right),$$

(12)

where $q$ is the number of parameters; and N is the number of samples. logLik is the maximum likelihood under the model.

For the validation of the generalized models, the mean square error (MSE) was used to compare the variation between the best base H–D models of differentiated and undifferentiated tree species and the generalized H–D models with all expanded variables. The specific MSE was calculated as follows:

$$\mathrm{M}\mathrm{S}\mathrm{E}=\frac{1}{N}\sum _i^N{(H_i-{\widehat{H}}_i)}^2,$$

(13)

where MSE is the estimated mean square error; $H_i$ represents the observation of the dependent variable in the validation dataset (Zone II); ${\widehat{H}}_i$ indicates the fitted value of the dependent variable using the model obtained from the training dataset; and $N$ is the number of trees. For the mixed-effects models, the model evaluations and comparisons used logLik, AIC, and BIC. More importantly, we used the generalized mixed-effects H–D model for Zone II to predict the carbon.

All the data processing in this study was performed using the nlme package of R software (version 4.1.3; R Core Team, 2022; Vienna, Austria).

3. Results

3.1. Base Model

The results of the base H–D model showed that for all trees (without dividing the tree species), BM.9 was the best with the lowest RMSE, CV-MSE, AIC, and BIC (

Table 4. Fitting results of the base empirical model.

Model	RMSE	CV-MSE	AIC	BIC
BM.1	2.8805	2.8651	12,757.48	12,775.05
BM.2	2.8834	2.8690	12,765.07	12,782.64
BM.3	2.8876	2.8718	12,769.31	12,786.88
BM.4	2.9004	2.8844	12,791.41	12,808.97
BM.5	2.9554	2.9443	12,898.41	12,915.98
BM.6	2.8800	2.8653	12,758.24	12,775.81
BM.7	3.0826	3.0643	13,102.85	13,120.42
BM.8	2.9355	2.9187	12,852.45	12,870.02
BM.9	2.8785	2.8635	12,754.90	12,772.47 1
BM.10	3.2619	3.2426	13,393.53	13,411.10

¹ Bold indicates the best performing model.

). BM.9 was used as the initial model for generalized and mixed-effects model building. The specific equation is as follows:

$$H=1.3+e^{\left[3.165+\frac{-11.3}{D+1})\right]}$$

(14)

The best base models for the ten tree species were inconsistent when dividing the tree species to be modeled individually (Table S3). After testing and comparing the validation dataset, the modeling for whole species was better than that for single species (Figure 3).

3.2. Trait-Based Model Test Results

Figure 4 shows the univariate linear relationship between each functional trait and the H–D ratio. All functional traits can significantly affect the H–D ratio (p < 0.05). Specifically, only the relationship between LT and LDMC was negative, while the other relationships were positive.

The superiority of the model was improved when the traits were considered in the base model, both for single properties and mixed traits (PC) (Table S4). In addition, there were differences in the same trait variable for different reparameterizations, and these also appeared in the results for other expanded variables (Table S4).

The generalized H–D model in the best-reparametrized form for each trait showed a higher accuracy than that of the base model in the validation dataset (Figure 3). Note that, despite the excellent performance of Hmax in the training dataset among all the traits, its validation, even less than the estimation of tree height using the base model. Therefore, in the later multilevel combinations of the expanded variables, mixed trait indicators (PC) were used to represent the expanded variables at the species level. The specific form of a generalized model with PC is as follows:

$$H=1.3+e^{\left[3.1502-0.0188\times PC+\frac{(-10.9549+0.5868\times PC)}{D+1}\right]}$$

(15)

3.3. Different Levels of Generalized H–D Comparison and Combination

The model was further improved by adding a single expanded variable at the stand level or individual level, except for the SIC, which showed no convergence problem (Table S4). Among them, the stand quality at the stand level showed the best fitting result, and only it could result in better predictions than PC among all the expanded variables at every level.

When combining the traits with other variables, the effect of different combinations varied (

Table 5. Fitting and validation results for all generalized H–D models. The forms used here are the best reparametrized forms for each expansion variable. n.s. indicates that the parameters of the model were not significant (p < 0.05).

Model Type	Expanded Variable	RMSE	NMSE	AIC	BIC	MSE
Base model	Null	2.8785	2.8635	12,754.90	127,72.47	10.3711
Generalized model for one expanded variable	CE	2.8571	2.8584	12,749.62	12,773.04	10.4462
	BA	2.8427	2.8424	12,721.31	12,744.74	10.3414
	DMH	2.7225	2.7239	12,493.77	12,523.05	9.8701
	PC	2.8348	2.8370	12,709.06	12,738.34	10.1739
	BAL	2.8525	2.8525	12,739.78	12,763.20	10.4187
Generalized model for mixed extra variable	Based on level (DMH + PC + BAL)	n.s.	n.s.	n.s.	n.s.	n.s.
	Based on feature 1 (CE + BA + DMH + PC + BAL)	2.6704	2.6188	12,398.10	12,456.66	9.7457
	PC + CE	2.8268	2.8037	12,702.12	12,737.25	10.2478
	PC + BA	2.8111	2.7905	12,672.88	12,708.02	10.1766
	PC + DMH	2.6965	2.6576	12,440.17	12,481.16	9.7707
	PC + BAL	2.8214	2.7993	12,691.55	12,726.69	10.2587

¹ Bold indicates the best performing model.

). Among them, the feature-based combination form (considering all the expanded variables) showed the best accuracy, which was also the case in the validation dataset. The specific equation is as follows:

$$H=1.3+e^{\left[2.4712+0.0286\times DMH-0.0169\times PC+\frac{\left(-7.8520-0.1101\times DMH+0.5677\times PC-0.1402\times CE+3.6273\times BA-3.7214\times BAL\right)}{D+1}\right]}$$

(16)

Compared to the generalized model with one expanded variable, PC showed a better fit when combined with any other single expanded variable.

3.4. Mixed-Effects H–D Development and Examination

After building the mixed-effects model with plots as random effects, the final form of the generalized mixed-effects model is as follows:

$$\begin{gathered} H=1.3+e^{\left[2.4893+{\theta }_1+0.0316\times DMH-0.03\times PC+\frac{\left(-7.8879+{\theta }_2-0.1996\times DMH+0.8028\times PC-0.1099\times CE+0.9703\times BA+0.0844\times BAL\right)}{D+1}\right]}+\epsilon \\ \epsilon ~N\left(\mathrm{0,5.5781}\right) \end{gathered}$$

(17)

where ${\theta }_1$ and ${\theta }_2$ are the random-effects parameters for ${\beta }_1$ and ${\beta }_2$ separately; and $\epsilon$ is the within-plots error term. The residuals of this model do not have weak heteroskedasticity (Figure 5).

The observed value of the aboveground carbon stock in Zone II is 63,365.55 Kg. When a generalized mixed-effects model is used for carbon stock estimation in this zone, the result is 64,031.85 Kg (an error of 1.0515%). The prediction accuracy improves by 82.42% from the traditional approach (the base model for single tree species) (with an estimated value of 67,154.95 Kg and an error of 5.9802%).

4. Discussion

The aim of this study was to develop a generalized mixed-effects H–D model that included functional traits in a temperate mixed forest. For our comparison, we chose BM.9 as the base model, which has a higher prediction accuracy in the natural forest than even the base model divided into tree species (Table 4). BM.9 has a logistic S-shaped growth trend and is consistent with the conclusions of previous studies on the allometric growth relationship between tree height and diameter [57], which has also been demonstrated in other similar H–D models [58,59]. Moreover, BM.9 was used as the base model for estimating tree height in the northern part of Idaho, USA, within the same latitude range [38]. Therefore, in this study, all generalized models were extended based on BM.9.

We first validated the effect of functional traits on allometric growth, and then developed a set of multilevel generalized mixed-effects models incorporating functional traits and validated their applicability in non-modeling zones. This will provide an essential reference for applying ecological and plant physiological mechanisms in forest management.

4.1. Modeling Thoughts and Methods for Natural Mixed Forests

In this study, both modeling ideas for natural mixed forests were evaluated. Surprisingly, the test results in the new stand (Zone II) showed that the modeling method without a division into tree species had a smaller MSE compared with that of the modeling method for single tree species (Figure 3), which implies that the division of tree species for separate models was not appropriate for our natural mixed forest.

There are two possible reasons for this result: first, for some tree species, the allometric growth relationship between tree height and diameter is not stable in complex natural mixed forests, and large variations may arise within the same diameter class [60]. In contrast, the modeling approach that includes all species reduces variation by considering a wide range of tree species characteristics, such that similar results are seen in studies using traits to predict growth [61]. Previous studies have found that using traits at the species level enables growth rates to be predicted more accurately than at the individual level, and this could be due to the fact that the species level is more representative of a tree’s specific growth potential. More importantly, despite the fact that the validation dataset is in the same climate region as the training dataset and also has similar stand types, the H–D relationship varies considerably with the stand characteristics and competition [23,58,62]. The effects of these factors may be much higher than the effects at the species level. This will be further illustrated later in this paper.

The idea behind and approach for modeling single tree species have been validated in plantation forests and natural forests with clearly dominant species (target management species) [10,63,64]. However, in natural mixed forests with complex tree species, many studies only treat tree species as dummy variables or random effects in their model [46,60] and achieve satisfactory fitting results. This implies that tree species should not be considered to be the main factor affecting the H–D model in natural mixed forests. In addition, the tree species in different natural forest communities are so diverse and different that it is impractical to consider all species for modeling individually, which would be costly in terms of human and economic costs [9].

In conclusion, modeling ideas for plantations and pure natural forests should not be directly applied to mixed forests models [65]. Our results demonstrate a basic principle in natural mixed forest modeling: the focus should be on the abiotic effects and biotic interactions for allometric growth relationships rather than on species tags. This idea deserves to be further validated and extended in future studies of natural mixed forests.

4.2. Functional Traits Are Worthy of Attention in the H–D Model

Our results showed that the functional traits of plants could influence the H–D ratio (Figure 4). Typically, the tree height is considered to be a characteristic of a plant’s ability to obtain light and compete, while the diameter as the support structure represents its ability to survive and defend. Therefore, the H–D ratio could represent the trade-off between individual growth and survival [66]. Similarly, functional traits were demonstrated to characterize the fitness of species under different environmental changes. In this study, we used seven functional traits representing the main components of different economic spectra [40,67,68]. SLA, LA, and LN could indicate the characteristics of resource acquisition and allocation by plants, and increasing values of these traits reflect the high nutrient acquisition and photosynthetic capacity of the species. Therefore, the trade-off in allometric growth was biased towards height growth. In addition, the relationship between the WD and H–D ratios was positive. This finding broadly supports those of other studies in the field linking wood density to slow radial growth [41]. It is thus clear that a direct linear prediction of the allometric growth relationship between the tree height and diameter using traits is not appropriate (the R² is low) [69]. However, due to the significant relationship between the H–D ratio and traits (p < 0.05), traits could still be used as extended variables in generalized nonlinear H–D models to enhance the model’s superiority.

The most important finding of this study is that functional traits can represent species characteristics applied in the generalized H–D model with a higher accuracy than that of the basic model, both for individual traits and when combined with other explanatory variables (Table 5). Among all individual traits, wood density had the highest predictive performance, with a 5.15% reduction in error in the validation dataset compared to the base model (MSE: 10.0347 vs. 10.5796). Wood density, an important growth and survival trade-off for trees, expresses the efficiency of transporting water and nutrients [70]. Moreover, it represents a tree’s ability to resist external disturbances [67]. These characteristics are closely related to the trade-off between plant height and radial growth.

The functional traits of plants have long been used to predict plant growth and survival [71,72,73]. In recent years, there has been considerable controversy over whether traits are suitable for evaluating ecosystem functioning [74,75,76]. Functional traits are excellent tools for revealing plant responses to the environment as well as species interactions [77], including allometric growth [78]. Ecological theories and methods can offer a broad perspective on forest management. However, only biodiversity seems to have received significant attention in studies related to the H–D model [55,60,78], and no studies have yet tested functional traits in tree growth and harvest models.

As discussed above, the differences between species for H and DBH allometric growth relationships have been considered, but few studies have quantified species-level effects (see Hulshof et al., 2015 [13]). It is generally accepted that species effects on H–D models originate from their genetic and biological characteristics [79,80], and this influence can only be classified by replacing the base model form. Our results confirm that functional traits as the phenotypic attributes of plants can be used to quantitatively estimate allometric growth relationships and growth relationships across species.

It is worth noting that, as inferred from previous studies, functional traits cannot be used arbitrarily due to the presence of intraspecific variation [81,82], as in our study, the Hmax led to inaccurate predictions (with errors larger than the base model) during validation. This may be attributed to the fact that the Hmax of species is more subject to intraspecific variation due to environmental influences in different stands [14], which are more suitable to be investigated separately in different stands and used as an indicator to assess stand quality [10,55,63]. In contrast, other traits (e.g., leaf economic spectrum traits) are more stable at small scales of environmental variation [83]. To prevent similar accuracy problems in practical applications, we suggest using a composite trait indicator for species transformed via a principal component analysis (PC), as it represents the characteristics of species on each functional axis well [68]. We used such composite trait values in our multilevel expanded variables’ combination and mixed-effects modeling.

4.3. Multi-Level Expanded Variables’ Application and Combination

We considered expanded variables at the stand, species, and individual levels separately, and the combination of variables at different levels formed a generalized mixed-effects H–D model that reduced the MSE by up to 7.9% (Table 5). This hierarchical concept has only been applied in a few studies [60]. Adding expanded modeling information has become necessary to build more accurate H–D models [44,84]. However, most studies have only focused on the magnitude of the effects of a single type of expanded variable, without considering the hierarchy processing as well as combinations across the hierarchy. To build a more stable H–D model, we optimized two aspects: one is the filtering out of expanded variables that are representative at each level, such as competition at the stand level (BA) and the individual level (BAL), which has been used in generalized H–D models in several regions [23,85]; and the second is topography as an environmental variable which better explains allometric growth differences in tree height and diameter among stands at small scales [55]. We also considered all reparameterizations, and although some studies reparametrized for only one specific parameter (e.g., the asymptote parameter), our results suggest that more parameters deserve to be considered in the reparameterization process.

In this study, stand quality was the most effective expanded variable, and a generalized H–D model that included only the stand quality could improve the accuracy of the base model by 6.7% (Table 5). Stand quality can indicate the production potential of a forest, and it is widely used in growth and harvest models as a basic characteristic representing a specific stand [86]. However, some studies have argued that stand quality does not explain the variation pattern of tree height and diameter [87,88]. We believe that this may be caused by the inconsistency of indicator arithmetic. Moreover, since stand quality is a relatively abstract indicator, it needs to be expressed using specific stand survey factors [89,90] by referring to previous related studies, in which the dominant tree height was used as the stand quality and achieved significant accuracy gains [60,91]. Similarly, we also obtained good results using the maximum tree height of the dominant species per sample site to indicate the stand quality.

Covering multiple levels of information will maximize model accuracy, but it may be difficult to obtain all of the expanded variables during actual field surveys and applications. When this information is lacking, we recommend using at least the traits of the species and the stand quality of the stand to improve the prediction accuracy of tree height in natural mixed forests.

4.4. Review and Outlook

We conducted a series of tests and have presented arguments for the H–D model for natural mixed forests, and our results suggest that functional traits should be considered in the study of the relationship between tree height and diameter at the species level. Despite our promising results, questions remain. For example, the direct relationship between the functional traits and allometric growth was not close, which may be due to the fact that the selected traits did not cover the overall developmental strategy characteristics of the trees, or it may be that the functional traits used in this study were the average values of major tree species in forests in Northeast China. A large number of studies have also been based on the mean values of traits for data processing and mechanism studies [32,33]. Further research should be undertaken to investigate the intraspecific trait variation in different stands so that traits in different stands cover both stand-level and species-level information. This approach will obviously increase the workload of researchers investigating plant traits and would thus be more suitable for the conservation and management of specific vegetation areas. Forest growth and harvest models based on average traits can be directly extended to a broader spatial and geographical area. In addition to this, due to study scale limitations, it is not possible to further determine the suitability of our model for carbon stock estimations at larger scales, especially in the context of uncertain global climate change. Therefore, the application of functional traits should be extended to larger temporal and spatial scales in future research on allometric growth modeling [92].

Some modeling principles, such as the partitioning of datasets, should be considered in future research [7]. To test the actual performance of the developed models, some data methods (e.g., cross-validation) can greatly reduce errors [44]. However, the training trees and validation trees of the model should not come from the same location or dataset and should be extended using more places to demonstrate that a modeling approach or idea is generalizable [10,46,55,85].

With the development and progress of science and statistical technology, data Bayesian methods, machine learning, and deep learning [93] can be useful under various conditions. However, a greater emphasis on biological and ecological mechanisms is also required to meet the challenges of estimating the effects of climate change and human disturbance [94].

5. Conclusions

This study presents a new series of height–diameter models for a typical natural mixed-species forest in Northeastern China. We provide evidence that tree species’ functional traits can significantly influence tree H–D ratios and are suitable as expanded variables for a generalized H–D model. Moreover, expanded variables at different levels have different effects on the relationship between height and diameter and can be combined into a generalized model covering information from the individual to the stand level. Our mixed-effects combined model with a sample plot as a random effect can be applied to new stands and can improve the prediction accuracy of aboveground carbon stocks. Our results show that functional traits have a significant effect on the H–D ratio and can be applied to the H–D model and extended to similar allometric growth relationship studies, in which traits are used as a species-level stable factor to overcome the difficulty of modeling various tree species in natural mixed forests. In addition, the generalized mixed-effects model combining multilevel variables can improve the aboveground carbon stock estimation of the zone up to 82.42%.