Modeling Height–Diameter Relationships for Mixed-Species Plantations of Fraxinus mandshurica Rupr. and Larix olgensis Henry in Northeastern China

Abstract

The mixture of tree species has gradually become the focus of forest research, especially native species mixing. Mixed-species plantations of Manchurian ash (Fraxinus mandshurica Rupr.) and Changbai larch (Larix olgensis Henry) have successfully been cultivated in Northeast China. Height–diameter (H–D) models were found to be effective in designing the silvicultural planning for mixed-species plantations. Thus, this study aimed to develop a new system of H–D models for juvenile ash and larch mixed-species plantations, based on competition information and tree and stand attributes. The leave-one-out cross-validation was utilized for model validation. The result showed that the H–D relationship was affected not only by the tree attributes (i.e., tree size and competition information) but also by stand characteristics, such as site quality and species proportion of basal area. The best model explained more than 80% and 85% variation of the tree height of ash and larch, respectively. Moreover, model validation also confirmed the high accuracy of the newly developed model’s predictions. We also found that, in terms of total tree height, ash in middle rows were higher than those in side rows, while larch in the middle rows were higher in the early growth period but then became lower than those in the side rows, as the diameter increased. The newly established H–D models would be useful for forestry inventory practice and have the potential to aid decisions in mixed-species plantations of ash and larch.

1. Introduction

One-third of the Earth’s total surface area is covered by forests [1], which act as a major sink for carbon and has an essential role in maintaining the global climate system and carbon cycle [2]. Natural forests occupy 93% of the global forest area, with a declining annual rate of 6.5 million hectares per year (2010–2015), while the other 7% were occupied by plantations with an average annual growth of 3.1 million hectares per year (1990–2015) [1,3]. Plantations are gradually expanding in many countries, including China. According to the ninth forest resource survey report (2014–2018) [4], the total forest area of China has reached 220 million hectares (22.96%), where as much as 79.54 million hectares were identified as plantations. Thus, strong theoretical knowledge and comprehensive forest management planning is necessary to manage such a large plantation area.

To date, numerous studies have been found regarding both mixed-species plantations [5,6,7,8,9] and monoculture [10,11,12] along with their pros and cons. Kelty [13] and Piotto [14] revealed that trees in large-scale monocultures could provide higher yields per hectare with certain desirable wood quality and growth rates. Additionally, monocultures of fast-growing and low-density wood species used for the fuel, pulp, and paper industry took up a large proportion of plantations [15]. However, Liu et al. [7] criticized monoculture plantations due to their negative environmental impacts, including the loss of soil productivity and fertility, disruption of hydrological cycles and biodiversity, etc. As opposed to monocultures, mixed-species stands were found to be more resilient [16] and productive [17], while still providing higher biodiversity [18] and an enormous potential to supply ecological and socio-economical services [5,19,20,21]. However, it needs considerably more effort in finding the best combination of species. Species mixing could increase maximum stand density and biomass yield [20,22], especially for native species [6,12], since they have a better adaptation ability to local conditions and are considered to be more commercially and environmentally attractive. Therefore, one of the most important steps for successful mixed-species plantations is to properly choose the combination of two or more tree species, which might generate a positive mixture effect [23].

Manchurian ash and Changbai larch are the two most crucial fast-growing tree species in Northeastern China [24]. These species are in different shade-tolerance levels, which create a perfect fit to be cultivated as mixed-species plantation [25,26]. Previous literature revealed that mixed planting of ash and larch have a higher survival rate of sapling [27], more productive than their respective monocultures [28], and improve the soil fertility and microbial activity [29]. Wang et al. [30] also found that the net photosynthetic rate (NPR) of ash was higher in mixed-species than in monoculture plantations, while at the same time, it did not cause any reduction on the NPR of larch. Unfortunately, further application and promotion of mixed-species plantation management are limited by a lack of comprehensive studies and effective equations.

In forest resource inventory, total tree height (H) and diameter at breast height (D) are the two most basic factors, which hold crucial implications for biomass and carbon estimation, along with forest growth and yield models [31]. Generally, tree height is measured by hypsometer or LIDAR data, which are labor-intensive, time-consuming, and costly. In this context, D of all trees and H of the sub-sample trees are mostly measured in a sample plot for forest inventory [32]. Thus, height–diameter (H–D) models that can predict H developed on the allometric relationship between H and D are needed. H–D model is necessary and meaningful in terms of volume and biomass estimation, stand structural analysis, growth and yield system, and carbon model. It can help forest managers better understand the nature of various relationships that characterize, differentiate, and influence forest development [33]. To date, numerous H–D models were developed across different regions [34], stands, and species [35,36,37,38]. Simple H–D models that only utilize D as a sole predictor are unable to adapt to diverse growth environments. Hence, it had relatively low precision and was not preferable among forest practitioners [39].

An effective way to improve model suitability is by developing a generalized H–D model. Generalized H–D models have shown excellent performance in describing H–D relationships, either in even- or uneven-aged stands [40,41,42,43]. It used stand variables and tree social status and competition as covariates to describe site and environmental conditions [32]. These covariates might account for the diversity of tree height with a given diameter, which makes the H–D models more applicable. In mixed-species plantations, it is essential to consider the mixing effect by some special indices, for characterizing stands structures such as stand density, species proportion, species competition, etc. [44].

Developing H–D models for mixed-species plantations of ash and larch is more challenging for the relationship of species coexistence. The aims of this study were—(1) to develop and compare H–D models based on different predictors and methods; (2) to evaluate the fitting and prediction performance of different approaches for H–D models; (3) to analyze the effect of tree- and stand-level variables on the tree height used to develop the H–D models.

2. Materials and Methods

2.1. Study Sites

This study was conducted in six forest farms (the Maoer Mountain Forest Farm, Yimianpo Forest Farm, Xiaojiu Forest Farm, Shangzhi Forest Farm, Heilonggong Forest Farm, and the Yuanbao Forest Farm) in the Shangzhi city of the Heilongjiang province, Northeast China (from 127°17′ E to 129°12′ E and from 44°29′ N to 45°34′ N, Figure 1). The Shangzhi city has a total area of 8910 km² with varying landscapes, from high mountains in the east to undulating low hills in the west, and the altitude ranged from 116 to 1600 m above sea level. Dominated by dark brown soil, the region had a temperate continental monsoon climate characterized by less rain and low temperature in spring, warm and rainy summer, long and cold winter, and autumns that brought frequent early frosts and a steep decline in temperature. It had a low mean annual temperature of 2.5 °C, with the highest- and lowest-monthly temperatures being 21.6 °C in July and −20.2 °C in February, respectively. The Shangzhi City had an average annual rainfall of 652.2 mm, with single rain events ranging from 5 to 577 mm. The total area of the forest land was approximately 5100 km², out of which nearly 1300 km² were plantations. The dominant tree species for the plantations were Larix olgensis Henry, Pinus koraiensis Siebold & Zucc., Pinus sylvestris var. mongholica Litv., and Fraxinus mandshurica Rupr.

In the present study, the mixed-species plantations of Manchurian ash and Changbai larch in the Shangzhi city were relatively young; their stand age ranged from 10 to 25 years. They were established at an initial density of 3300 cutting·ha⁻¹, with reasonably narrow spacing (1.5 m $\times$ 2.0 m) and several variations of row-wise mingling patterns (Figure 2). The planting designs were distinguished as seven categories:

(A)

1-row larch: 1-row ash (1:1)

(B)

2-rows larch: 2-rows ash (2:2)

(C)

3-rows larch: 2-rows ash (3:2)

(D)

3-rows larch: 3-rows ash (3:3)

(E)

5-rows larch: 3-rows ash (5:3)

(F)

5-rows larch: 5-rows ash (5:5)

(G)

6-rows larch: 4-rows ash (6:4)

Based on those various planting designs, both larch and ash were divided into two categories—(1) trees in the middle of the rows (denoted as TMR); and (2) trees on the sides of the rows (denoted as TSR). TMR represents the situation in which the neighboring rows were trees of the same species, while TSR refers to the condition in which the adjacent rows were trees of different species.

2.2. Data

In July and August of the three consecutive years (2017, 2018, and 2019), 69 sample plots from 23 plantations (3 sample plots per plantation) were established, with a permanent length of 50 m and the width of at least three repeats, which varied according to the planting designs. The so-called “repeat” was determined as a basic unit for the width of each sample plot, and one repeat consisted of some rows in accordance with the summation of each planting design’s ratio. For example, for 1:1 planting design, one repeat equaled to two rows containing an equivalent proportion of both species. Meanwhile, for 3:2 planting design, one repeat was equivalent to five rows comprising three rows of larch and two rows of ash. A similar mechanism was also performed for other planting designs. Further, several basic information for all individual trees located in the sample plots was recorded, such as tree species, status (dead or alive), and type of rows (TMR or TSR). Then, D and H were measured using diameter tapes and Vertex IV Ultrasonic Hypsometer made by Haglöf Sweden, respectively.

After the basic measurement was conducted, the following stand variables were measured—(1) regardless of the tree species, such as the number of trees per hectare (N, trees·ha⁻¹), quadratic mean diameter (Dq, cm), the dominant height of 100 D-based largest trees per hectare (Hd, m), basal area (BA, m²·ha⁻¹), and the dominant diameter (Dd, cm) of 100 D-based largest trees per hectare; (2) by taking into account the tree species, such as the number of trees per hectare (N_L for larch and N_F for ash), quadratic mean diameter (Dq_L for larch and Dq_F for ash), basal area (BA_L for larch and BA_F for ash), dominant height (Hd_L for larch and Hd_F for ash), and dominant diameter (Dd_L for larch and Dd_F for ash). In order to describe the specific structure of the mixed-species stand, we also calculated other simple indices at stand level [44]—(1) the proportion of each species relative to total trees (RN_L = N_L/N for larch and RN_F = N_F/N for ash); and (2) the basal area proportion of each species relative to the total basal area (RBA_L = BA_L/BA for larch and RBA_F = BA_F/BA for ash). In addition, the information on the neighboring trees represented the competitive state of the selected tree, which had an important effect on the H–D relationship. Hence, we calculated the distance-independent indices—the ratio of each individual D to Dd (RD) [45] and the basal area of the trees larger than the subject tree (BAL, m²·ha⁻¹). Additionally, the corresponding intra- and inter-specific indices were also calculated: RD_intra, RD_inter, BAL_intra, and BAL_inter. Summary statistics of the tree and stand variables are presented in

Table 1. Summary statistics of the mixed-species stands along with the stand- and tree-variables.

Variable		Sample Sizes	Mean	Min.	Max.	S.D.
Stand level	Age (years)	69	18	10	25	4
	N (trees·ha−1)	69	1843	1060	2583	374
	NL (trees·ha−1)	69	765	240	1693	318
	NF (trees·ha−1)	69	1079	353	1717	312
	Dq (cm)	69	11.2	7.4	15.2	1.8
	DqL (cm)	69	12.1	7.3	15.8	2.0
	DqF (cm)	69	10.6	6.5	16.1	1.8
	Hd (m)	69	15.5	9.7	19.4	2.4
	HdL (m)	69	14.9	7.7	19.4	2.5
	HdF (m)	69	15.1	9.2	19.1	2.4
	Dd (cm)	69	16.81	9.63	21.51	2.70
	DdL (cm)	69	16.29	8.55	21.14	2.71
	DdF (cm)	69	14.68	8.55	21.13	2.56
	BA (m2·ha−1)	69	17.84	7.59	24.92	3.63
	BAL (m2·ha−1)	69	8.50	1.61	16.83	3.32
	BAF (m2·ha−1)	69	9.34	2.98	15.80	2.82
	RNL	69	0.41	0.17	0.73	0.14
	RNF	69	0.59	0.27	0.83	0.14
	RBAL	69	0.47	0.18	0.79	0.14
	RBAF	69	0.53	0.21	0.82	0.14
	RNLF	69	0.82	0.21	2.70	0.55
	RBALF	69	1.05	0.22	3.79	0.67
Tree level	Larch in side rows
	D (cm)	5190	11.4	5.0	25.6	3.4
	H (m)	5190	12.0	4.9	20.6	3.0
	RD	5190	0.69	0.25	1.48	0.17
	RDintra	5190	0.71	0.26	1.49	0.17
	RDinter	5190	0.78	0.29	1.73	0.20
	BAL (m2·ha−1)	5190	9.74	0.00	24.83	5.66
	BALintra (m2·ha−1)	5190	5.69	0.00	16.82	3.38
	BALinter (m2·ha−1)	5190	4.05	0.00	15.55	3.21
	Larch in middle rows
	D (cm)	4479	12.3	5.0	25.6	3.3
	H (m)	4479	13.2	5.1	20.0	2.3
	RD	4479	0.68	0.26	1.55	0.17
	RDintra	4479	0.70	0.26	1.63	0.18
	RDinter	4479	0.78	0.29	1.72	0.20
	BAL (m2·ha−1)	4479	10.93	0.00	24.74	6.21
	BALintra (m2·ha−1)	4479	6.71	0.00	16.80	3.79
	BALinter (m2·ha−1)	4479	4.22	0.00	14.56	3.15
	Ash in side rows
	D (cm)	9460	10.5	5.0	25.3	3.0
	H (m)	9460	12.8	5.4	20.0	2.5
	RD	9460	0.63	0.25	1.60	0.15
	RDintra	9460	0.72	0.25	1.90	0.17
	RDinter	9460	0.65	0.25	1.71	0.17
	BAL (m2·ha−1)	9460	11.76	0.00	24.91	5.37
	BALintra (m2·ha−1)	9460	5.73	0.00	15.79	3.41
	BALinter (m2·ha −1)	9460	6.03	0.00	16.83	3.29
	Ash in middle rows
	D (cm)	3699	10.1	5.0	22.0	2.5
	H (m)	3699	13.3	5.7	19.9	1.9
	RD	3699	0.56	0.25	1.08	0.14
	RDintra	3699	0.66	0.28	1.25	0.15
	RDinter	3699	0.58	0.26	1.20	0.14
	BAL (m2·ha−1)	3699	13.94	0.37	24.87	4.88
	BALintra (m2·ha−1)	3699	6.68	0.00	14.67	2.86
	BALinter (m2·ha−1)	3699	7.26	0.23	16.83	3.18

S.D.: Standard Deviation. N = number of trees per hectare; N_L = number of trees per hectare for larch; N_F = number of trees per hectare for ash; Dq = quadratic mean diameter; Dq_L = quadratic mean diameter for larch; Dq_F = quadratic mean diameter for ash; Hd = dominant height; Hd_L = dominant height for larch; Hd_F = dominant height for ash; Dd = dominant diameter; Dd_L = dominant diameter for larch; Dd_F = dominant diameter for ash; BA = basal area; BA_L = basal area for larch; BA_F = basal area for ash; RN_L = relative N_L for larch; RN_F = relative N_F for ash; RBA_L = relative basal area for larch; RBA_F = relative to basal area for ash.

.

2.3. Methods

2.3.1. The Comparison of Total Tree Height (H) Between TMR and TSR

Previous studies revealed that in ash and larch mixed-species plantations, H of TMR and TSR were significantly different under the same specific stand age [27,28]. Prior to the development of H–D models, we utilized the probability density curves of H in different rows, to compare the differences of H between TMR and TSR. Furthermore, the differences were examined statistically by using analysis of variance (ANOVA) or nonparametric tests, at a significance level of 0.05. The choice of using either ANOVA or nonparametric tests, depended on the normality test of each sample. If the normality test was successfully passed, ANOVA would be used; otherwise, nonparametric tests would be more appropriate.

2.3.2. Basic Model of H–D Relationship

In this study, an H–D allometry relationship was visually screened to appropriately choose 12 equations as candidate functions (Figure 3 and

Table 2. Height–diameter (H–D) equations selected for evaluation.

Functions Number and Forms	Reference
[1] $H={\beta }_1+{\beta }_2\mathit{log}D$	Curtis, 1967
[2] $H=1.3+{\beta }_1\left(1-e^{-{\beta }_2{D}^{{\beta }_3}}\right)$	Weibull, 1951
[3] $H=1.3+{\beta }_1{(1-{\beta }_2{e}^{-{\beta }_{3}D})}^{{\beta }_4}$	Richards, 1959
[4] $H=1.3+{\beta }_1{\left(\frac{D}{1+D}\right)}^{{\beta }_2}$	Curtis, 1967
[5] $H=1.3+{\beta }_1{D}^{{\beta }_2}$	Schreuder, 1979
[6] $H=1.3+\frac{{\beta }_{1}D}{D+{\beta }_2}$	Bates, 1980
[7] $H=1.3+e^{({\beta }_1+\frac{{\beta }_2}{D+1})}$	Wykoff, 1982
[8] $H=1.3+{\beta }_1{e}^{\frac{{\beta }_2}{D}}$	Schumacher, 1939
[9] $H=1.3+{\beta }_1\left(1-e^{-{\beta }_{2}D}\right)$	Farr, 1989
[10] $H=1.3+{\beta }_1{e}^{(\frac{-{\beta }_2}{D+{\beta }_3})}$	Ratkowsky, 1990
[11] $H=1.3+{\beta }_1{D}^{{\beta }_2{D}^{-{\beta }_3}}$	Huang, 1992
[12] $H=1.3+{(\frac{D}{{\beta }_1+{\beta }_{2}D})}^2$	Näslund, 1936

) [42,44,45,46,47,48,49,50,51,52,53,54,55,56]. Furthermore, the 12 selected models were fitted to the entire dataset by utilizing the base package in R software 3.5.3 [57], and the best fitting performance model was named as the basic H–D model (HD_B) for ash and larch.

2.3.3. Generalized H–D Model

The H–D relationship was strongly influenced by the stand characteristics and neighboring trees’ information [32]. Thus, a generalized H–D (GHD) model was necessary. Once the best H–D model was selected, we evaluated the potential of several stand variables and competition indices contributions to the H–D relationship. First, the selected basic model was fitted to the data for each plot to obtain the corresponding parameter estimation values. Furthermore, the relationships between model coefficients and stand variables or competition indices were scrutinized by scatter plot and correlation analysis [58]. The stand and competition variables needed to be significant, and the fitting statistics of the GHD model needed to be improved.

Dummy variable models are frequently used in forestry modeling, to describe the effect of species-, region-, and even origin-specific differences [40,59,60]. In this study, a dummy variable ($S_1$) was included in the GHD model to describe the row-specific effect ($S_1=0$ for TMR and $S_1=1$ for TSR), which was called a GHD_R model. An F-test was then used to test whether separate models were required for the different type of rows, from a statistical point of view [59]. Parameters of the GHD model were expressed by a linear function $S_1$, which was defined as follows:

$${\beta }_i={\beta }_j+{\beta }_k{S}_1$$

(1)

where ${\beta }_i$ is the parameter in the GHD models; β_j and β_k are the estimated parameters in the GHD_R model.

2.3.4. Mixed-Effects Model

The mixed-effects model is a model that includes both fixed- and random-effect variables. The application of a mixed-effect model has increased in the field of forest modeling [40,41]. In this study, a nonlinear mixed-effects (NLME) modeling framework was used for the hierarchical structure of H–D data, by setting the sample plots as random effects in HD_B, GHD, and GHD_R. Numerous studies have applied the mixed-effects models to describe H–D relationships and have improved the model fitting and prediction accuracy [32,43]. At a plot-level, the NLME model of the jth tree height in the ith sample plot was modeled as [61]:

$$\{\begin{array}{c}{H}_{ij}=f\left({\varphi }_{ij},{\nu }_{ij}\right)+{\epsilon }_{ij}; \\ i=1,2,3, \dots ,m; j=1,2,3, \dots , n_i\end{array}$$

(2)

where $H_{ij}$ represents the height of the jth tree in the ith plot; $m$ and $n_i$ are the number of plots and observations in the ith plot, respectively; $f$ is a nonlinear function of a plot-specific parameter vector ${\varphi }_{ij}$ and a covariate vector ${\nu }_{ij}$; and ${\epsilon }_{ij}$ is a within-plot error term, which subjects to a multivariate normal distribution with a mean value vector of 0 and variance–covariance matrix of $R$. ${\varphi }_{ij}$ is given as:

$${\varphi }_{ij}=A_{ij}\beta +B_{ij}{b}_i; b_i~N\left(0,G\right)$$

(3)

where $\beta$ is the fixed-effect parameter vector; $b_i$ is the random-effect parameter vector of the ith plot, which was assumed to have a multivariate normal distribution with a mean value vector of 0 and variance–covariance matrix of $G$, $A_{ij}$, and $B_{ij}$ are the incidence matrices of the appropriate dimensions, consisting of 0 or 1. In the weighted regression analysis, a power function was used to correct the heteroscedasticity of variance exhibited in the H–D models, and to find the most effective variance functions [58]. All regressions were calculated in the R software 3.5.3., by using both base and nlme packages [57,62].

2.3.5. Model Assessment and Validation

In the present study, six types of models were evaluated—(1) HD_B model; (2) GHD model; (3) GHD_R model; (4) mixed—effect HD_B model (HD_BM); (5) mixed-effect GHD model (GHD_M); and (6) mixed-effect GHD_R model (GHD_RM). These models were validated using the leave-one-out cross-validation [63], in which the H–D models were constructed with all-but-one sampled plot data (fitting data included m-1 sample plots) and the fitted model was used to predict the excluded sampled plot. Statistical indices were used to assess and compare model performance, and the best model was selected, based on the following criteria:

Statistical significance of the estimated parameters;

Statistics of model assessment, such as the adjusted coefficient of determination (R_a²), root mean square error (RMSE), and Akaike information criterion (AIC) [64].

Performance of model validation, such as mean absolute error (MAE) and mean absolute error percent (MAE%) [65].

The marginal adjusted coefficient of determination (R_am²) only involved fixed effect, while the conditional adjusted coefficient of determination (R_ac²) concerned both fixed- and random-effects, which was further used to evaluate the mixed-effects models [40]. However, the mixed-effects model validation needed prior information to calculate the random-effects parameters. The parameters of random effects were calculated using the best linear unbiased predictions (BLUPs) method [61]. A vector of random effects parameters of sampled plot k was calculated with Equation (4).

$$b_k=G{Z}_k^T{(Z_{k}G{Z}_k^T+R_k)}^{-1}{\epsilon }_k$$

(4)

where $b_k$ is a vector of random effects parameters of sampled plot k calculated by all sample trees; $G$ is the variance–covariance matrix estimated in the modeling process; $R_k$ is the corresponding variance–covariance matrix of within-group errors; $Z_k$ is the matrix of the partial derivatives of the nonlinear function, with respect to its random parameters; and ${\epsilon }_k$ is error terms of the sampled plot k using the fixed effects parameters of the mixed-effects models.

2.3.6. Comparison of Sample Designs

Although validation of the mixed-effects models used all sample tree height observations to calculate the random effect parameters, different numbers of available measured trees were then used to compare the prediction performance. In calculating the random effect parameters, researches have pointed out that the increasing number of sampled trees per plot would decrease the bias between the observations and the predicted values [64]. However, a large number of sample trees per plot would lead to a higher inventory cost, which only gives a slight improvement in the model’s accuracy. Thus, four trees per sample plot seemed to be adequate for the calculation of random effects [41]. Simulations were performed using subsamples of tree height measurements, where predefined sample sizes of 2, 3, 4, 5, 6, 7, 8, 10, 12, 14, 16, 18, 20, 24, and 26 trees were subsampled without replacement and the calculations were repeated continuously for 100 times. Within each subsample, we estimated the random effect parameter values and obtained the final prediction bias. These simulations were also performed in R software 3.5.3. using the nlme package [62].

3. Results

3.1. The Variation of H in Different Rows

The probability density plots of the two species were used to analyze the distribution of H across different statuses (TMR and TSR, Figure 4). Preliminary visual inspection showed that the two species’ heights were obviously different across both TMR and TSR, in which, the data centered on different values. As the non-normality (p < 0.01) of H in different rows of two species showed in Figure 4, a nonparametric test (Mann-Whitney test [66]) was applied to compare whether the value of the average height of TMR and TSR were significantly different between larch and ash. The result of the two-way Mann-Whitney test supported that the average height of TMR and TSR were significantly different (p < 0.01). However, limited evidence confirmed the different height of TMR and TSR by the data collected so far; hence, it was also necessary to use a modeling method for further analysis and verification.

3.2. Basic H–D Model Results

By fitting the 12 candidate functions listed in Table 2, R_a² ranged from 0.5780 to 0.6011 for ash and from 0.6295 to 0.6513 for larch (

Table 3. Fitting statistics of candidate functions.

Function NO.	Larch			Ash
	Ra2	RMSE	AIC	Ra2	RMSE	AIC
[1]	0.6490	1.6274	36551.17	0.5949	1.4980	48000.48
[2]	0.6505	1.6238	36509.48	0.5997	1.4891	47845.45
[3]	0.6511	1.6225	36495.03	0.6009	1.4869	47807.14
[4]	0.6509	1.6230	36499.19	0.5996	1.4893	47847.50
[5]	0.6295	1.6720	37069.57	0.5780	1.5289	48538.71
[6]	0.6406	1.6467	36776.98	0.5917	1.5039	48104.73
[7]	0.6502	1.6245	36517.13	0.5985	1.4913	47883.26
[8]	0.6513	1.6222	36489.04	0.6004	1.4878	47820.95
[9]	0.6428	1.6416	36717.57	0.5951	1.4977	47995.33
[10]	0.6513	1.6220	36488.64	0.6011	1.4864	47797.80
[11]	0.6510	1.6227	36496.60	0.6003	1.4880	47825.59
[12]	0.6479	1.6299	36580.49	0.5971	1.4939	47928.45

The function in bold is the best basic model. R_a² = adjusted coefficient of determination; RMSE = root mean square error; AIC = Akaike information criterion.

), which showed a relatively slight fluctuation in a small range. The other two fitting statistics (RMSE and AIC) also displayed a similar pattern. Table 3 indicated that the number 10th function from Ratkowsky was found to be better than others for both larch and ash, with a higher value of R_a² and a lower value of RMSE and AIC.

3.3. Generalized H–D Model

Correlation analysis of each parameter in the HD_B model indicated that the dominant height considering species and species proportion by basal area had a relatively greater impact on the parameters. Different competition indices had different effects on species; therefore, the final forms of the GHD model were:

Ash:

$$H=1.3+\left({\beta }_0+{\beta }_{1}H{d}_F+{\beta }_{2}RB{A}_F\right)e^{\frac{-\left({\beta }_3+{\beta }_{4}BA{L}_{intra}+{\beta }_{5}BA{L}_{inter}\right)}{D+{\beta }_6}}$$

(5)

Larch:

$$H=1.3+\left({\beta }_0+{\beta }_{1}H{d}_L+{\beta }_{2}RB{A}_L+{\beta }_{3}RD\right)e^{\frac{-{\beta }_4}{D+{\beta }_5}}$$

(6)

There was a substantial improvement in model fits when the covariates were added, (

Table 4. Parameters estimates, variance, and the fit statistics for ash.

Terms	HDB	GHD	GHDR	HDBM	GHDM	GHDRM
Fixed Parameters
${\beta }_0$	19.3542 (0.2054)	0.9403 (0.2014)	1.0989 (0.2045)	15.957 (0.3153)	0.0192 (0.2353)	0.0187 (0.2343)
${\beta }_1$	4.4967 (0.1812)	0.8947 (0.0081)	0.2319 (0.0277)	2.7639 (0.0998)	0.9842 (0.0105)	0.2901 (0.0299)
${\beta }_2$	−0.9106 (0.1621)	−0.6741 (0.1585)	0.8878 (0.0081)	−1.484 (0.1696)	−1.2728 (0.1861)	0.9773 (0.0103)
${\beta }_3$		0.7846 (0.0978)	−0.6955 (0.1597)		0.9683 (0.1858)	−1.2811 (0.185)
${\beta }_4$		0.0552 (0.0035)	0.8572 (0.1025)		0.0386 (0.0074)	0.8929 (0.1847)
${\beta }_5$		0.0431 (0.0033)	0.0545 (0.0036)		0.0437 (0.0114)	0.0390 (0.0074)
${\beta }_6$		−2.5907 (0.1104)	0.0443 (0.0034)		−2.5601 (0.1757)	0.051 (0.0112)
${\beta }_7$			−2.4928 (0.1149)			−2.5679 (0.1745)
Variance and Covariance
${\sigma }_{u_{i1}}{}^2$				6.0896	0.0022	0.0021
${\sigma }_{u_{i1}{u}_{i2}}$				−1.8709	0.0187	0.0170
${\sigma }_{u_{i2}}{}^2$				0.7876	0.6604	0.6326
${\sigma }^2$	1.9840	0.6526	0.6733	0.3357	0.3180	0.3214
Power of Variance Function
φ	0.2268	0.4106	0.3841	0.4412	0.4620	0.4565
Fitting Statistics
Ram2				0.5367	0.7739	0.7755
Rac2	0.6007	0.7799	0.7811	0.8073	0.8073	0.8086
RMSE	1.4866	1.1037	1.1007	1.0327	1.0329	1.0293
AIC	47781.92	39948.66	39878.79	38201.22	38205.99	38115.64

HD_B = H–D model for ash and larch; GHD = generalized H–D model; GHD_R = GHD model to describe the row−specific effect; HD_BM = mixed—effect HD_B model; GHD_M = mixed-effect GHD GHD_RM = mixed-effect GHD_R model

and

Table 5. Parameters estimates, variance and fitting statistics for larch.

Terms	HDB	GHD	GHDR	HDBM	GHDM	GHDRM
Fixed Parameters
${\beta }_0$	22.0687 (0.32)	−3.2305 (0.3365)	−2.8436 (0.3659)	22.1981 (0.3766)	0.8813 (0.9818)	1.3128 (1.0437)
${\beta }_1$	7.3134 (0.2993)	0.9926 (0.0132)	−0.6354 (0.1554)	10.7985 (0.6782)	1.1193 (0.0231)	−0.6404 (0.1943)
${\beta }_2$	−0.3217 (0.2003)	−0.8614 (0.1186)	1.0035 (0.0146)	3.7726 (0.4841)	−2.7711 (0.3547)	1.1434 (0.0254)
${\beta }_3$		4.5559 (0.2281)	−0.9276 (0.1225)		2.1637 (0.4704)	−2.8744 (0.3693)
${\beta }_4$		2.7125 (0.288)	4.5772 (0.2324)		6.1743 (0.7434)	1.9964 (0.4935)
${\beta }_5$		−0.7514 (0.3519)	3.2783 (0.3478)		1.7318 (0.4936)	6.7962 (0.8233)
${\beta }_6$			−0.5494 (0.1211)			−0.5872 (0.1477)
Variance and Covariance
${\sigma }_{u_{i1}}{}^2$				12.5038	0.6341	0.6223
${\sigma }_{u_{i1}{u}_{i2}}$				6.6646
${\sigma }_{u_{i2}}{}^2$				4.9502
${\sigma }^2$	0.8771	0.8660	0.8956	0.6397	0.6729	0.6708
Power of Variance Function
φ	0.4953	0.1894	0.1615	0.1922	0.1740	0.1748
Fitting Statistics
Ram2				0.5849	0.8361	0.8365
Rac2	0.6513	0.8399	0.8404	0.8590	0.8568	0.8572
RMSE	1.6220	1.0992	1.0975	1.0314	1.0393	1.0382
AIC	36488.65	29029.72	29002.23	27810.19	27957.6652	27937.7944

). The R_a² increased from 0.6011 to 0.7800 for ash and from 0.6513 to 0.8418 for larch. The corresponding RMSE and AIC also showed a significant reduction. Stand variables and competition indices explained more variations on the basis of HD_B. Then, a dummy variable $S_1$ was added into the different parameters, and the fitting performance was evaluated by a smaller AIC value and the statistical significance of parameters. The results of the F-test between the GHD and GHD_R models turned out to be statistically significant for both ash and larch (p < 0.01). The final forms of the GHD_R models for both species were as follows:

Ash:

$$H=1.3+\left({\beta }_0+{\beta }_1{S}_1+{\beta }_{2}H{d}_F+{\beta }_{3}RB{A}_F\right)e^{\frac{-\left({\beta }_4+{\beta }_{5}BA{L}_{intra}+{\beta }_{6}BA{L}_{inter}\right)}{D+{\beta }_7}}$$

(7)

Larch:

$$H=1.3+\left({\beta }_0+{\beta }_1{S}_1+{\beta }_{2}H{d}_L+{\beta }_{3}RB{A}_L+{\beta }_{4}RD\right)e^{\frac{-\left({\beta }_5+{\beta }_6{S}_1\right)}{D+{\beta }_7}}$$

(8)

Based on the GHD_R models’ formula, the effects of $S_1$ (TMR and TSR) on the H–D relationship for larch and ash were simulated at the average values of other predictors with three random sampled plots (Figure 5). The GHD_R model results showed that there was a significant difference in H between TMR and TSR for both species. The H of TMR was found to be higher than those of TSR for both ash and larch. However, it is worth noting that for larch, the H of TMR was higher in the early growth period, but then the difference decreased with time, before finally being lower than the H of TSR.

3.4. Mixed-Effects H–D Model

The plot-level random effect was introduced in the HD_B, GHD, and GHD_R models; and was found to be significantly different among the three models, in which a variation of the H–D models was observed among the different plots. However, employing a plot-level random effect on both the GHD and GHD_R models for the two species showed a slight improvement. The HD_BM model described nearly 20% more variations than the HD_B model in the H–D relationship, and there was an apparent difference between the fit performance evaluated by only considering the fixed-effect parameters and the mixed-effect parameters (R_am² and R_ac² were 0.5367 and 0.8073, respectively, for ash, and 0.5849 and 0.8590 for larch). The R_am² values of GHD_M and GHD_RM models were slightly smaller compared to the corresponding R_ac² values. The final form of the HD_BM, GHD_M, and GHD_RM models were:

HD_BM models:

Ash:

$$H=1.3+\left({\beta }_0+u_{i1}\right)e^{\frac{-\left({\beta }_1\right)}{D+{\beta }_2+u_{i2}}}$$

(9)

Larch:

$$H=1.3+\left({\beta }_0\right)e^{\frac{-\left({\beta }_1+u_{i1}\right)}{D+{\beta }_2+u_{i2}}}$$

(10)

GHD_M models:

Ash:

$$H=1.3+\left({\beta }_0+{\beta }_{1}H{d}_F+{\beta }_{2}RB{A}_F\right)e^{\frac{-\left({\beta }_3+{\beta }_{4}BA{L}_{intra}+\left({\beta }_5+u_{i1}\right)BA{L}_{inter}\right)}{D+{\beta }_6+u_{i2}}}$$

(11)

Larch:

$$H=1.3+\left({\beta }_0+{\beta }_{1}H{d}_L+{\beta }_{2}RB{A}_L+{\beta }_{3}RD\right)e^{\frac{-{\beta }_4}{D+{\beta }_5+u_{i1}}}$$

(12)

GHD_RM models:

Ash:

$$H=1.3+\left({\beta }_0+{\beta }_1{S}_1+{\beta }_{2}H{d}_F+{\beta }_{3}RB{A}_F\right)e^{\frac{\left({\beta }_4+{\beta }_{5}BA{L}_{intra}+\left({\beta }_6+u_{i1}\right)BA{L}_{inter}\right)}{D+{\beta }_7+u_{i2}}}$$

(13)

Larch:

$$H=1.3+\left({\beta }_0+{\beta }_1{S}_1+{\beta }_{2}H{d}_L+{\beta }_{3}RB{A}_L+{\beta }_{4}RD\right)e^{\frac{-\left({\beta }_5+{\beta }_6{S}_1\right)}{D+{\beta }_7+u_{i1}}}$$

(14)

where $u_{i1}$ and $u_{i2}$ are the random parameters, while the other symbols are defined above. The parameter estimates, variances, and fit statistics for all H–D models are listed in Table 4 and Table 5. The scatterplots between the standardized residuals and the estimated height by utilizing weighted regression are shown in Figure 6. Within the sampled plot, heteroscedasticity in the residuals was reduced by applying a weighted regression in the form of a power function, at the base of the predicted height. The values of the power variance functions are also listed in Table 4 and Table 5.

The effects of Hd, RBA, and BAL considering tree species along with Rd on the H–D relationship were simulated with H–D curves by species-specific GHD_RM model of ash (Figure 7) and larch (Figure 8). The variables of interest were roughly divided into four equal intervals, where other predictors were at the mean value of our data. The simulation showed that all variables, except Hd_F, had a negative effect on H for both species, which was similar to others [40]. As shown in Figure 7 and Figure 8, Hd gave a more considerable contribution to the H–D model, among all variables. The differences of H between RBA_F and RBA_L were larger as the D increased, and the species competition significantly affected the initial H–D relationship for ash and increased the influence of the competition to larch.

3.5. Model Validation

The GHD_RM model was found to have the best performance among the six models, followed by GHD_M, HD_BM, GHD_R, GHD, and HD_B. The GHD_RM model gained the lowest value of MAE (0.8095 m for ash and 0.8111 m for larch) and MAE% (6.4927% for ash and 6.7883% for larch) (

Table 6. Indices of model validation for larch and ash.

Model	Ash		Larch
	MAE (m)	MAE% (%)	MAE (m)	MAE% (%)
HDB	1.1967	9.8683	1.3046	11.5190
HDBM	0.8299	6.6421	0.8789	7.6137
GHD	0.8701	7.0340	0.8686	7.3402
GHDM	0.8156	6.5440	0.8117	6.7934
GHDR	0.8682	7.0166	0.8684	7.3366
GHDRM	0.8095	6.4927	0.8111	6.7883

MAE = mean absolute error.

). The results showed that the fit statistics and the performance of model validation of all models did not have a significant increase (Table 4 and Table 5). For larch, apart from the HD_B model, the value of MAE ranged from 0.8095 to 0.8299, and MAE% decreased from 6.6421% to 6.4927%, and a similar regular pattern also appeared for ash. Overall, the prediction performances were within a reasonable range and improved as the model’s complexity increased.

The residuals of the six models were plotted across four D classes (Figure 9). The average values of MAE and MAE% were relatively stable in all D classes and six models, except the HD_B models for both small size (5 < D < 10 cm) ash and larch, which indicated a stable prediction accuracy within the range of D. This could improve the prediction accuracy by applying both the join of covariates and NLME modeling approach for larch and ash. Compared to the ordinary least square, the NLME method showed a significant improvement from the perspective of MAE and MAE%.

3.6. Result of the Sampled Designs Comparison

The prediction performance of the GHD_RM model for each species with different numbers of prior measured trees are shown in Figure 10. The value of MAE and MAE% were decreased with increasing numbers of sampled trees, which implicated that the prediction accuracy increased when more prior height information was used to calculate the random effect parameters. There was no significant difference between six, seven, or even more sampled trees, based on the significant test of MAE and MAE% in the different numbers of previously measured trees.

4. Discussion

4.1. Basic H–D Model and Generalized Model

In order to develop the H–D model for mixed-species plantations of ash and larch, 12 alternative models based on D as the only predictor were evaluated in this study. The function number 10 in Table 3 and Table 4, which were previously used to predict tree height in British Columbia [67], gave the best performance to be fitted in our data. In this function, ${\beta }_0$ was considered to be the asymptotic height, and ${\beta }_1$ represented steepness, known by its derivative function that would be extended to the formula of stand variables and competition information, to develop a generalized model.

H_d mostly represented site quality and was frequently used in the H–D models for both monocultures and mixed-species plantations [32,40]. Species proportion by basal area is an important attribute in mixed-species plantations [44], which had an effect on stand structure and tree growth. Beside stand variables, competition information expressed by some indices, including distance-independent (i.e., BAL) [32,68], is necessary for the H–D model. The distance-dependent competition indices might show a better performance in describing the competition status. However, it is time-consuming and labor-intensive, since each neighboring tree’s coordinates need to be measured precisely. Due to the limitation of both fieldwork duration and economic resources, we chose the indices of BAL and RD to explain the competition variations in the H–D models. When these covariates were added into the basic H–D models of ash and larch, R_ac² increased by 34.39% and 31.89%, while RMSE reduced by 30.94% and 36.41%, respectively. Hence, the extended basic H–D model was preferable, since it gave a relatively high prediction accuracy. The results of the model validation by utilizing the leave-one-out cross-validation also showed that the generalized H–D models were better than the basic models, in terms of prediction. Thereby, H of the mixed-species plantations of ash and larch would be accurately estimated with the newly developed models, which were found to be more flexible and effective.

GHD_R used a dummy variable to analyze the differences between TMR and TSR for both species, which proved to be significantly different by the F-test result between GHD_R and GHD model (p < 0.01). Wang et al. [26] found that the mixture of ash and larch improved intraspecific competition of ash, which might be the reason why H of ash in middle rows was higher than those in the side rows, in three different random sample plots, for the same value of D (Figure 5). In contrast to ash, H of larch was found to be higher in TMR than TSR in small D; then the H of larch was found to be higher in TSR than TMR as the D increased, which was similar to others [27,28]. Larch is classified as shade intolerant-species. In mixed-plantations of larch and ash, larch in TSR was often shaded by ash, primarily the middle and lower part of their crown. As a result, larch in TMR got more lights and grew faster, mainly in the early phase after planting [69]. In the present study, we found that there was a large number of dead larches in mixed-species plantations, especially for those that grew in the side rows. One possible reason behind this phenomenon was the effect of interspecific relationship, since the remaining larches in TSR were found to grow healthier with broader space.

4.2. H–D Model with NLME

The NLME model considers the potential differences between plot-levels and is widely applied in H–D modeling [40,41,64]. The performances of the GHD, GHD_R, and HD_BM models were quite similar, which indicated that predicting a random effect can be thought of as an empirical surrogate for stand-level variables that affect the H–D relationship [68]. Generally, the convergence of mixed-effects models become harder to achieve as the number of random parameters increase. It was almost impossible to achieve the convergence when the models contained more than two random parameters, in which further improvement on the algorithm was required [70]. Actually, the mixed-effects model by only considering the fixed parameters could be used to predict height [40], but the prediction accuracy would be significantly lower and was not recommended to be applied. In our study, the R_am² of HD_BM, GHD_M, and GHD_RM were lower than the same models estimated by nonlinear least squares. Therefore, the predictions of the NLME models containing the random effects were better and had relatively high precision.

As shown in Table 6, the prediction ability of the mixed-effects model was better by using the random effect parameters calculated by all measurements of H. However, it was very difficult or even unrealistic to measure the H of all trees in forest resource inventory, which implied that choosing several trees for H measurement would be more favorable. Many studies have analyzed the impact of different numbers of available measured trees on the prediction of mixed-effects model [40,41,43]. The prediction performance of mixed-effects model validation with a priori samples of different sizes was compared, and the result showed that the model’s performance increased as the prior information increased. Based on the significant test of MAE and MAE% in Figure 10, six trees would be a reasonable choice and close to the number of dominant trees per sample plot, as slightly different from other studies [41,43].

4.3. Influence on Species Mixing in H–D Relationship

One of the advantages of mixed-species over monoculture plantings is that it yields a higher productivity [7], and the mixed-species plantations of ash and larch in Northeast China is thought to be a more productive plantation than their corresponding monocultures [28,30]. Zhang et al. [28] concluded that H of both species in mixed-species plantation were higher than in monoculture plantations. The final H–D model for larch in this study differed from the model formula developed by Zang et al. [32] using data from monoculture plantations of larch across a broader geographic area; both included stand variables and competition indices that improved fitting performance and increased the precision of predictions. The predictions of H utilizing H–D model developed by Zang et al. [32] were significantly lower compared to both the observed and predicted H by using the model developed in this study (Figure 11). Therefore, preliminary analysis showed that H of larch in mixed-species plantation of ash and larch were higher than that in monoculture plantations of larch. However, a further comparison of the H–D model using more detailed data between ash in mixed-species and monoculture plantations is essential, since an H–D model in monoculture plantations is still lacking.

Generally, the H–D models developed in this study show strong practicability and operability in predicting height by relying on diameter and dominant tree height data, which could be used to calculate other covariates. Therefore, the developed models in this research were useful for predicting tree height and in future work, it will be an essential tool for analyzing the differences in total biomass, stand volume, wood quality, etc. of mixed-species plantations and monoculture plantations.

5. Conclusions

Tree height has important implications on stem volume, biomass, carbon storage, and others, with a growing interest in mixed-species plantations. The main objective of the developed H–D models in this study was to predict the missing height data in mixed-species plantations, which might reduce field workload and improve the prediction accuracy. Compared to previous research work on H–D models, we considered the row-specific effect along with intra- and inter-specific competition effects to describe the H–D relationship in mixed-species plantation of ash and larch. The GHD model included D as a predictor, and two covariate predictors (stand variables and competition information), which showed a significant improvement as compared to the basic models. According to the actual condition, we could use the different H–D models to predict tree height. The three mixed-effects models HD_BM, GHD_M, and GHD_RM performed well, but the use of GHD_RM was preferable. However, the newly developed H–D models was not only important to predict the tree height but was also essential to comprehend the vertical structure of mixed-species plantation of ash and larch in Northeast China.