Dagujia forest farm and Changlinggang forest farm, characterized by different climate characteristics, were the locations to conduct our field investigation. Dagujia forest farm (124°47′–125°12′E, 42°22′–42°16′N) is from Liaoning province of Northeast China, and Changlinggang forest farm (110°00′–110°04′E, 30°47′–30°50′N) is from Hubei province of central China. The present study was conducted at these two forest farms because the two sites are nearly the northernmost and southernmost areas suitable for the normal growth of Larix kaempferi plantation. The elevation of Dagujia forest farm ranges from 200 to 600 meters. The annual average temperature is 6.3°C, the extreme maximum temperature is 34°C, and the extreme minimum temperature is −34°C. The average annual precipitation is 600 mm, mainly concentrated from July to September. Frost-free days in 1 year is about 137 days, and the type of soil is mainly Alfisol. Changlinggang forest farm is characterized by a north subtropical monsoon mountain temperate climate. The range of altitude is 1,600–1900 m. The annual average temperature is 13.5°C, the extreme maximum temperature is 29°C, and the extreme minimum temperature is −17.2°C. The average frost-free days per year is about 250 days. The annual precipitation is 1,600 mm, mainly concentrates from July to September. The type of soil is mountainous brown soil. Larix gmelinii, Pinus massoniana, and Cunninghamia lanceolata are the main tree species.

An investigation at Dagujia forest farm was conducted in 2017 and 2018. A total of 17 forest stands with complete appearance and different ages were selected to develop permanent sample plots. In 2017, a total of 96 permanent sample plots, each with an area of 0.08 hectare (approximately 28.3 × 28.3 m²), were developed from 14 forest stands with ages ranging from 4 to 34 years. In 2018, a total of nine permanent sample plots with an area of 0.06 hectare for each plot (20 × 30 m²) were developed from the remaining three forest stands with age ranging from 12 to 20 years. Thus, 105 sample plots were developed with current densities from 388 to 2,850 trees per hectare, and with an initial planting density of 2,500 or 3,300 trees per hectare. The field investigation at Changlinggang forest farm was conducted from June to August in 2019 and 2022. In 2019, six permanent sample plots, each measuring 0.06 hectares (20 × 30 m²), were developed for each of the forest stands aged 5–42 years. In 2022, a supplementary field study was conducted at three additional forest stands aged 19–36 years. Six plots, each sized 0.06 hectare (20 × 30 m²), were developed from each forest stand. Therefore, a total of 78 permanent sample plots were developed at Changlinggang forest farm, with the current density ranging from 183 to 1,667 trees per hectare, and an initial planting density of forest stands of 1,667 or 2,500 trees per hectare. The detailed description of the initial planting density and historical forest management practices for each sample plot are shown in Table 1.

The measurements for the attributes of all trees from each plot were conducted and recorded in the field. Diameter at breast height (DBH, cm), tree height (HT, m), and height to the crown base (HCB, m) were measured. Crown width (CW, m) from four directions was measured with a measuring tape. Height-diameter ratio (HDR) was calculated as the ratio of total tree height to diameter at the breast height. Dominant tree height (HT_dom, m) was defined as the average tree height for the eight trees with the largest tree height from the 0.08-hectare plots and six trees from the 0.06-hectare plots. The sum of the basal areas of trees larger than the subject tree (BAL, m²·ha⁻¹), the number of trees with DBH larger than the subject tree (LDN, trees·ha⁻¹) and the stand basal area per hectare (BA, m²·ha⁻¹) defined as the ratio of the sum of basal area for each tree in the plot-to-plot area were calculated. Trees from each plot were screened for the preliminary inspection, and those with unreasonable measurements or other tree species in the sample plots were excluded. Only the plots with complete records of thinning times and thinning intensity were used in the height to crown base model development. Finally, a total of 87 plots from Dagujia forest farm and 78 plots from Changlinggang forest farm were used to fit the height to crown base model. Table 2 shows detailed information about the stand characteristics of each sample plot. A total of 13,822 trees from 165 sample plots were used to develop the height to crown base models for Larix kaempferi plantations.

Among all the predictors, height to crown base is mostly related to tree height. Therefore, height to crown base is first modeled by using total tree height as the sole variable. Logistic function and exponential function (Wykoff et al., 1982; Soares and Tome, 2001) were used as base models, and the three models finally used and compared were listed in Table 3. All three base equations were fitted by ordinary least squares estimation. Ra², AIC, and RMSE were used to compare models (Equations 1–3).

Tree variable, competition variable, site quality, climate variable, initial planting density, and thinning intensity were added to the base model. All the variables and combinations of different variables have been tried to incorporate into the model. The fitting accuracy of each model and significance of parameter test were combined to select the best model. In addition, the biological sense of the parameter estimates was also considered. In detail, the tree level variables included HT, DBH, CW, and HDR. The competition variables considered were BA and BAL, site factor of HT_dom, logarithmic transformation of altitude (lnE), slope percentage multiplied by cosine of slope direction (SLC), slope percentage multiplied by sine of slope direction (SLS), and LDN. The climate variables were generated following Ma et al. (2024). Initial planting density was added using a dummy variable approach. To further analyze the effects of thinning on the change of height to crown base of individual trees, we further developed an index by considering the thinning intensity and the thinning intervals, and added it to the height to crown base model. Jia et al. (2021) developed an index for forest stands, conducted only once at the time of thinning, using a ratio between the thinning intensity and stand age when the thinning practice was implemented. In our study, we developed a comprehensive index for conducting two and three times of thinning by taking the thinning time effect into consideration. For the forest stands conducted thinning practice one, two, and three times, a comprehensive thinning effect could be shown as Equations 4–6.

Where TE_i were the thinning indices for the forest stands conducting thinning for one (i = 1), two (i = 2), and three times (i = 3), respectively; TI₁, TI_2, and TI₃ were the thinning intensities for the first time, second time, and third time, respectively; Age is forest stand age for the current time, Age₁, Age₂ and Age₃ were the forest stand ages for the first thinning, second thinning and third thinning.

To avoid the problem of over-parameterization results from the increasing number of variables included in the model, only one variable from each type was selected and added to the final model. The multicollinearity of the model was examined using the variance inflation factor (VIF). VIF is larger than 10, indicating a problem of multicollinearity between variables. Considering the hierarchical structure and the possible correlation between the trees nested in the plots, a nonlinear mixed effects model was developed by introducing random effects at the plot level (Yan et al., 2024; Yang and Huang, 2018). A residual plot was used to detect the efficiency of heteroscedasticity elimination by using the power equation and exponential equation. All analyses were conducted through the R software using nlme package (R Core Team, 2023). Hierarchical partitioning analysis was used to quantify relative importance of the independent variables on height to crown base performed by “hier.part” package of R software (Tong et al., 2025).

Based on estimates and biological sense of parameters, as well as goodness of fit for the three candidate models in Table 3, the logistic equation was regarded as the best base model with Ra² of 0.76 and RMSE of 2.48 m. HDR from tree level variables and stand basal area per hectare from competition variables were finally incorporated into height to crown base model. DD_0_sp of climate and logarithmic transformation of altitude of site quality variables were included in the final model. Subsequently, the thinning index was successfully incorporated into the height to crown base model. Finally, initial planting density and study area were added into the model using a dummy variable approach. The final model is shown as Equation 7, and the parameter estimates and goodness of fit statistics are shown in Table 4. The final model had Ra² of 0.9297, AIC of 46,789, and RMSE of 1.3137 m. All parameter estimates were significant at the significance level of 0.05, and they had biological significance. It indicated that competition, climate, site quality, and thinning practice significantly affected height to crown base. All parameters for the dummy variables were significantly different at the significance level of 0.05, indicating that initial planting density had a significant impact on the height to crown base of trees from Liaoning and Hubei provinces.

Where HCB is height to crown base; HT is total tree height; BA is stand basal area per hectare; HDR is height-diameter ratio; DD_0_sp is spring degree-days below 0°C; TE is the comprehensive index related to thinning practice. Y₁ = 1, Y₂ = 0, Y₃ = 0 indicated the initial planting density of 1,667 tree·ha⁻¹; Y₂ = 1, Y₁ = 0, Y₃ = 0 indicated the initial planting density of 2,500 tree·ha⁻¹; Y₃ = 1, Y₁ = 0, Y₂ = 0 indicated the initial planting density of 3,300 tree·ha⁻¹, respectively; X₁ = 1 and X₂ = 0 represented Liaoning province, X₁ = 0 and X₂ = 1 represented Hubei province.

A nonlinear mixed effects HCB model at the plot level for Larix kaempferi was developed. The combinations of more than two random effect parameters lead to non-convergence or insignificant estimates. Finally, the parameter of m₁, regarded as a random effect, had the largest Ra², and the smallest AIC and RMSE (Table 4). Compared to the base model, Ra² of the final mixed model was 0.9316, AIC was 46,642, and RMSE was 1.2954 m, which were better than the base model. Based on the likelihood ratio test (LRT), the performance of nonlinear mixed effects model was significantly improved compared to the nonlinear model (LRT = 151.42, p < 0.0001). By comparing the residual graphs of the base model and nonlinear mixed effects model, it was shown that heteroscedasticity was improved and the range of residuals was reduced (Figure 1).

Based on the parameter estimates for the fixed effect in the mixed model, we further analyzed the influence of each variable on height to crown base (Figure 2). The results revealed that height to crown base of individual trees increased with the increasing of total tree height (Figure 2A), competition level (Figure 2B), HDR (Figure 2C), and the thinning intensity (Figure 2F). Height to crown base decreased with the increasing of DD_0_sp (Figure 2D) and logarithmic transformation of altitude (Figure 2E). The results from Figure 2 also indicate that a large initial planting density will increase the height to crown base in both Liaoning province and Hubei province. Increased thinning intensity decreased the competition level, and the effect of interaction of thinning and competition on the height to crown base was also studied, as shown in Figure 3. It indicated that the height of crown base rose with the increase in thinning intensity for forest stands with all the kinds of initial planting densities. The intense competition level also increased the height to crown base. It firmly showed that the effect of competition on the height to crown base was reduced with an increase in thinning intensity, especially for trees from Hubei province.

Relative contribution for each variable (Figure 4) and different groups of the data divided by different study areas (Figure 5, A1,A2), initial planting densities (Figure 5, B1–B3), and different thinning intensities (Figure 5, C1–C4) were studied. It indicated that the relative contributions of each group of variables to height to crown base were as follows: tree size (41.7%) > competition (27.3%) > thinning (24.5%) > interaction of competition and thinning (3.2%) > climate (2.4%) > site quality (0.9%). Total tree height, stand basal area per hectare and thinning intensity were the most important variables in determining height to crown base with tree height (36.9%) > stand basal area (24.6%) > thinning intensity (20.71%) for Liaoning province, and total tree height (33.8%) > thinning intensity (22.28%) > stand basal area (26.76%) for Hubei province. For the trees from different initial planting densities, the most important variables for the trees from the stand with initial planting density of 1,667 trees per hectare were thinning intensity (31.53%) > total tree height (28.97%) > stand basal area (27.57%), total tree height (23.28%) > logarithmic transformation of altitude (20.76%) > thinning intensity (19.05) for trees from the stand with initial planting density of 2,500 trees per hectare, and total tree height (39.67%) > stand basal area (24.30%) > thinning intensity (21.06%) for trees from the stand with initial planting density of 3,300 trees per hectare. We divided the thinning intensity into four groups; the most important variables for most of the groups were tree height, stand basal area, and thinning intensity.