Development of Variable-Density Yield Models with Site Index Estimation for Korean Pines and Japanese Larch

: The objectives of this study were to develop site index and variable-density yield models for Korean red pine ( Pinus densiﬂora Siebold & Zucc.), Korean white pine ( P. koraiensis S. & Z.), and Japanese larch ( Larix kaempferi (Lamb.) Carri è re) in Korea. The data were collected between 2012 and 2021 from repeatedly measured empirical plots in each target stand in the North Central region of Korea: Gangwon and North Gyeongsang provinces. To develop the site index for each species, a dominant height growth model by species was developed using the Chapman–Richards function. The site index was computed with a base age of 40 years and used as an independent variable to predict the stand volume. To develop the variable-density yield models, three stand density variables, the stand age, and the site index were applied. The stand density variables used were the stand basal area, the number of trees per hectare, and the relative density. All the models were successfully developed with signiﬁcant parameters and reasonable ﬁt statistics, and the residuals analyzed presented unbiased scatter plots. Yield models with the stand basal area, or the number of trees, can be used to predict the stand volume. The yield model with relative density was ﬂexible to apply across the stand age because the input of the absolute stand density was not required. Model simulation and comparisons with other studies also supported the applicability of the models developed in this study. The models were found to be highly applicable for predicting and simulating these targeted stands, particularly in Korea. variable for yield model development. This study was performed to develop variable-density yield models using Demonstration of the stand volume and simulation comparison with previous studies The predicted stand volume for this study was calculated using Equation (9) with yield model type 3 based on the ﬁxed-effect in Table 4. To compare the simulation, the variables in previous studies were set to a moderate level, because the model by Yoon et al. [49] was based on the stand age, and the models by KFRI [37] and NIFoS [38] were based on the stand age and the site index. The site index was applied as 14, 16, and 18 from KFRI [37] and 14, 16, and 18 from NIFoS [38], respectively, the Korean red pine, the Korean white pine, and the Japanese larch. The site index for this study was to moderate level: 16, 18, and 22, respectively, each species. this study was


Introduction
The growth and yield model at the whole stand level has the advantage that information on growth and yield can be obtained from the variables used to predict the stand growth [1,2]. The yield table is the table summarizing the expected yield by age and site index [3,4]. It expresses the predicted yield by age in even-aged forests and is one of the oldest methodologies used for estimating the yield [1,5]. The yield table can be divided into three groups, namely the normal yield table, the empirical yield table, and the variabledensity yield table [6]. The normal yield table provides the yield amount expected by the age and site index for ideally normal stands or fully stocked stands [7]. Data for normal tables is usually collected by the forest inventory. When establishing plots, it is important to ensure that they are installed in forest stands where the stand density is evaluated as being fully stocked or in normal stands at various ages [2,5].
In addition, the empirical yield table is termed in contrast to the normal yield table and is applied to stands with an average stand density [8]. To develop an empirical yield table, data should be collected from stands with an average density. The advantage of the empirical yield table over the normal yield table is that it provides accurate information on the representative yield of the stand under silvicultural treatment [9][10][11]. Yield models at Data on the Korean red pine, the Korean white pine, and Japanese larch stands were collected from the region of Gangwon and the North Gyeongsang provinces in Korea. The total number of sites was 49 for Korean red pine, 52 for Korean white pine, and 45 for Japanese larch. Three plots were established together at each site with a plot size of 0.04-0.09 ha, in line with the thinning intensity. The plots were originally set up to monitor the growth and yield under three thinning treatments, the control, light thinning, and heavy thinning. After the plots were established along with the first forest inventory, the plots were re-measured at 3 year intervals. Additional details on the methods used in this study can be found in Lee and Choi [43,44]. The total inventory period for this study was from 2012 to 2021, but the number of remeasurements were not the same among all the sites; 2-4 times depending on the established year and plot maintenance. For reasons related to natural disasters such as forest disease, forest fire, and storm damage, some sites were not retained during the whole measurement period. Accordingly, the total collected data were used independently following the objectives of the model development on site index and the prediction of the variable-density yield. Since tree heights in some sites were measured only in one plot by site during the first measurement, the total number of available samples was different between the model types.
For the site index models, repeatedly measured empirical plot data were used instead of samples from felled trees. The models were more stable and accurate based on the increased number of samples rather than relying on the few standard tree samples collected and the stem analysis [32,54,55]. In this study, we defined the dominant height as the average height of the top 20% of tree heights in a plot, as reported by Lee and Choi [43]. The stand age was applied based on the number of rings plus 2 on the wood disc, at approximately 0.2 m above ground, which was collected for the stem analysis of the felled tree during the first inventory. The number of years with remeasurement were then added to the stand age to calculate the corresponding age for each measurement year. The number of samples with summary statistics used for the site index model is shown in Table 1. All the age ranges, from the young to the old stands were considered in developing the site index model. For the variable-density yield models, the independent variables analyzed were the stand age (AGE), the stand basal area (BA), the number of trees per ha (TPH), the relative density (RD), and the site index (SI) with the stand volume (VOL) as a dependent variable. The TPH was calculated based on the number of trees per plot times the expansion factor of the plot size by hectare. The BA was calculated based on the sum of the tree basal area (ba) in a plot times the expansion factor of the plot size by hectare. The RD was calculated using the coefficients of the stand density index derived from previous literature [43]. The detailed methodological concept and approach have been explained in detail in Lee [55] and Lee and Choi [43]. The specific calculation procedure was as per Equations (1)- (5). (1) where D (cm) is the tree diameter at breast height (1.2 m above ground), or dbh; ba is the tree basal area (m 2 ); D q is the quadratic mean diameter (cm); ba is the average tree basal area in a plot (m 2 ); π is the mathematical constant approximately equal to 3.14159; SDI max is the maximum stand density index (trees/ha at 25 cm of D q ); SDI is the stand density index (trees/ha); exp is the mathematical constant approximately equal to 2.71828; ln is the natural logarithm. a 0 and a 1 are the coefficients to calculate the SDI provided by Lee and Choi [43] and RD is the relative density. For the Korean red pine, the Korean white pine, and the Japanese larch, the coefficient a 0 is 30,241.86, 54,235.40, and 83,349.42 and a 1 is −1.0466, −1.2141, and −1.3727, respectively [43]. All other variables are as defined previously.
The SI was calculated using the parameters for each species based on the site index model developed in this study. The SI curves with the Chapman-Richards model were applied with a base age of 40 years, representing the dominant height at the base age in m. To calculate the stand volume, the individual tree volume in a plot was calculated by using the tree diameter at breast height and the tree height based on stem volume models as described in Equation (6) with the parameters provided in Lee et al. [56]. In cases where tree height was not observed for the first measurement, height-dbh allometry via the Chapman-Richards model was applied to estimate the tree height by species [54,57].
where v is the tree volume (m 3 ) and H is the tree height (m). For the Korean red pine, the Korean white pine, and the Japanese larch, the coefficient c 0 is 0.0228, 0.01097, and 0.05489 and c 1 is 0.00003506, 0.00003772, and 0.00003366, respectively [56]. The stand volume was calculated from the sum of the individual tree volumes in a plot on a hectare basis with the unit of m 3 ha −1 . The plot data for the variable-density yield model were all inventoried data, similar to the data for the site index model. However, some of the plots with a particularly young stand age and/or a high stand density were excluded. We removed the empirical stands with a stand age of <15 and/or a relative density of >1.3 among all the plots. This was predominantly to prevent any unexpected bias in residuals from the model fitting and to provide accurate parameter estimates. The number of samples with summary statistics used for the variable-density yield model is shown in Table 2.

Modeling Approach and Statistical Analysis
To develop the site index, the dominant height growth model over age was developed based on the Chapman-Richards function. This type of model has been widely applied in many countries due to the flexibility of the model curve as well as its accuracy [31][32][33][34]38,39,54,55,[57][58][59][60]. In our study, each plot-level observation at the site was used as a sample for the site index model.
Many studies have also applied the function to explain the effects of additional variables on dominant height growth, such as soil and climate factors, as well as the effect of the initial planting density, especially in plantations [33,61]. The parameters of dominant height growth were fitted onto the base Chapman-Richards function, which is comprised of three parameters as shown in Equation (7). To develop the model more conservatively and readily comparable with those in the previous studies, the anamorphic site index curves were applied as shown in Equation (8). The detailed derivation of the algebraic difference equation from Equation (7) to Equation (8) is explained in Lee et al. and Lee [32,54].
where H d is the dominant height (m); e, also known as Euler s number, is the same with exp as previously defined; a, b, and c are parameters, which refer respectively to the asymptote, the growth rate, and the shape (or inflection) of the original Chapman-Richards function; A 0 is the base age at 40 years. The other variables are the same as defined earlier. After developing the model and estimating the site index for each species, the variabledensity yield models were developed with three types of stand density variables, namely BA, TPH, and RD, in addition to AGE and SI. The transformation of the stand volume in the logarithmic scale was analyzed to control the heteroscedasticity. For the predictor of stand density, BA, TPH, and RD were examined with several types of transformation, including the reciprocal value, the square root value, and the logarithmic scale. With the best predictors selected, each yield model in this study was designated as type 1 with BA, where VOL is the stand volume (m 3 ha -1 ); b 0 -b 3 are the parameters estimated in this study; Density is BA (m 2 ha -1 ) for type 1, TPH (trees ha -1 ) for type 2, or RD for type 3 by yield model. Other variables are the same as defined earlier.
In the developed models, a mixed-effect modelling approach was applied by considering the nested plot design [44]. The site, or experiment location, was used as a random effect [44]. For the dominant height model of the Chapman-Richards function, the random effect was added as a constant within the term of parameter a rather than other parameters due to the convergence problem [33] as shown in Equation (10). For the variable-density yield model, the random effect was added as a constant term [44] as shown in Equation (11).
Owing to the nested plot design, the random-effect parameter (u) and the random error term (ε) are incorporated, whereas all other variables and parameters remain the same as previously defined. All the parameters were fitted using the nlme package in R statistical software [62].
To evaluate the performance of the model, conditional R-squared (R 2 ) was used, which describes the proportion of variance explained by both the fixed and random effects. The root mean square error (RMSE), the mean absolute error (MAE), and the Akaike information criterion (AIC) were examined as fit statistics. Residual plots were checked over the fitted values and all the other predictors. In the modelling stage, both the logarithmic and the arithmetic scale in the yield model were examined, and some of the residual plots selected were illustrated for unbiased prediction of our models in the main manuscript.
All the other residuals are provided in the Supplementary Files. When the predicted values were calculated from the logarithmic scale to the arithmetic scale, we corrected the residuals by adding the variance of the error term [13,63]. The developed site index models were displayed with the observations and the site index curves, and the suitability was checked. Based on our research findings on the major site index range by species, the yield model simulations were illustrated and compared by species, as an example, to examine the effect of the site index and the stand density as the stand age increases. The variable-density yield model outputs were compared with the results from other studies to assess their reliability.

Parameters Estimates for Site Index Models
The models of dominant height with the Chapman-Richards function fitted well for all three species in Equation (10) and all the fixed-effect parameters were statistically significant (p < 0.001). The standard deviation of the random effect by species was relatively similar, ranging from 2.2389 to 2.9209. The standard deviation of the residuals was largest for the Japanese larch, given that the mean, maximum, and range for the dominant height in the samples were the largest among the species.
Parameter a was the highest with 30.6312 for the Japanese larch and the lowest with 21.4977 for the Korean red pine (Table 3). It reflected the asymptote of the three species well with the order of the Japanese larch followed by the Korean white pine and the Korean red pine. The shape parameter c was the highest with 1.9108 for the Korean red pine and the lowest with 1.3682 for the Japanese larch. It indicated that the inflection point for the dominant height for the Japanese larch was the earliest among the three species, and the growth performance came first for the young stage, given that parameter c was the smallest. Note: a, b, and c are fixed-effect parameters via Equation (10) with the standard error in parentheses. All the fixed-effect parameters were highly significant (p < 0.001). std(u) is the standard deviation of the random effect in the experiment. std(ε) is the standard deviation of the residual in model performance. R 2 is the coefficient of determination. RMSE is the root mean square error. MAE is the mean absolute error. AIC is the Akaike information criterion.
In the overall fit statistics, the model for the Korean white pine performed the best, given that R-squared was the highest at 0.9651 and the RMSE and the MAE were the lowest at 0.8485 and 0.6464, respectively. Compared to the standard deviation of the dominant height in the sample data, the RMSE and the MAE were smaller for all the species, and the models were found to be valid for use in prediction. The residual plots were all unbiased for the fitted values in all three species (Figure 1). The residuals were also shown to be randomly dispersed, with an unbiased pattern over the predictor of stand age ( Figure S2 of the Supplementary Material).

Parameters Estimates for the Variable-Density Yield Models
In the variable-density yield model, the model fitting performed best with the Korean white pine, followed by the Japanese larch and the Korean red pine (Table 4). All the fixed effect parameters were highly significant (p < 0.001) except for the intercept of model type 2 for the Japanese larch (p = 0.068). The standard deviation of the random effect was different among model types, but it was the smallest in the models for Japanese larch. When the model types were compared within each species, model type 1 with BA was the best for all species, with R 2 greater than 0.99. The model with TPH was estimated with the lowest fit statistics for all three species compared with the other model types.
The parameter signs for all the model types were the same among species and logical in accordance with the growth characteristics. The reciprocal AGE variable presented the negative parameter sign, which indicated a positive estimate of growth. The parameters of SI and the logarithmic form for all the stand density variables, such as the BA, TPH, and RD, were analyzed with positive signs. This was performed logically with the stand development of the stand volume ( Table 4) Table 3 using Equation (10).

Parameters Estimates for the Variable-Density Yield Models
In the variable-density yield model, the model fitting performed best with t white pine, followed by the Japanese larch and the Korean red pine (Table 4). A effect parameters were highly significant (p < 0.001) except for the intercept of m  Table 3 using Equation (10). Table 4. Parameter estimates and fit statistics for the variable-density yield models to predict the stand volume at the whole stand level by species.

Species
Model Type

Independent Variables
Fixed Effects Random Effect Residual Fit Statistics  (11) with the standard error in parentheses. All the fixed-effect parameters are highly significant (p < 0.001) except for b0 of type 2 for the Japanese larch (p = 0.068). std(u) is the standard deviation of the random effect in the experiment. std(ε) is the standard deviation of the residual in the model performance. R 2 is the coefficient of determination. RMSE is the root mean square error. MAE is the mean absolute error. AIC is the Akaike information criterion.  Table 4 using Equation (11).

Predictions and Simulation
To illustrate the site index, the fixed-effect parameters b and c by species from Table  3 were applied to Equation (8). The mean, with the range from minimum to maximum of the site index by species, was 17.0 (8.7-26.6) for the Korean red pine, 18 (Figure 3). The proportion of observations within these ranges was 87%, 97%, and 93%, respectively, of the total observations by species. Considering the site index range for most of the observations, the Korean white pine was shown to have the most compact range, followed by the Japanese larch and the Korean red pine. Overall, the site index curves in this study were relatively stable and were in line with the general site index curves in previous studies from the Korea Forest Research Institute [37] and the National Institute of Forest Science [38].  Table 4 using Equation (11).

Predictions and Simulation
To illustrate the site index, the fixed-effect parameters b and c by species from Table 3 were applied to Equation (8). The mean, with the range from minimum to maximum of the site index by species, was 17.0 (8.7-26.6) for the Korean red pine, 18 (Figure 3). The proportion of observations within these ranges was 87%, 97%, and 93%, respectively, of the total observations by species. Considering the site index range for most of the observations, the Korean white pine was shown to have the most compact range, followed by the Japanese larch and the Korean red pine. Overall, the site index curves in this study were relatively stable and were in line with the general site index curves in previous studies from the Korea Forest Research Institute [37]   After examining the standard range for the site index by species, the w yield simulations were demonstrated over the stand age. One example showed of stand volume by site index with the fixed relative density (Figure 4a), and After examining the standard range for the site index by species, the whole-stand yield simulations were demonstrated over the stand age. One example showed the range of stand volume by site index with the fixed relative density (Figure 4a), and the other example showed the range by relative density with the fixed site index (Figure 4b). The range for the site index and the relative density was determined to express the wide range of variables, based on the current study for the site index and Lee and Choi [44] for the relative density. The simulation was presented up to age 80, encompassing a high proportion of the age range in the modeling dataset.   Table 4. In plot (a), the simulation used a varying site index with a moderate level of fixed relative density with a value of 0.6. The representative range for the site index found in this study was applied by species. In plot (b), the simulation used a varying relative density with a moderate level for the fixed site index by species. SI: site index. RD: relative density.
In the results presented in Figure 4a, the stand volume increased the fastest in the Japanese larch stand, but it later yielded results similar to those of the Korean white pine stand. Also, the stand volume was the lowest in the Korean red pine stand. In the result shown in Figure 4b, the stand volume was generally highest in the Japanese larch stand at an early stage, but the Korean white pine stand attained a similar stand volume after a stand age of 60 years was reached. The lowest stand volume was recorded in the Korean red pine stand. In Figure 4a, the minimum and maximum stand volume ranged between 210-332 m 3 ha −1 for the Korean red pine, 316-438 m 3 ha −1 for the Korean white pine, and 311-458 m 3 ha −1 for the Japanese larch at 50 years of age. Similarly, in Figure 4b, the stand volume ranged between 247-394 m 3 ha −1 , 301-499 m 3 ha −1 , and 320-548 m 3 ha −1 by species. This comparison indicates a greater difference in the stand volume derived from the stand density.
The variable-density yield models in this study were demonstrated in line with similar prior studies in Korea, where the stand yield was predicted by the stand age and the The simulations were performed using Equation (9) with model type 3 based on the fixed-effect parameters in Table 4. In plot (a), the simulation used a varying site index with a moderate level of fixed relative density with a value of 0.6. The representative range for the site index found in this study was applied by species. In plot (b), the simulation used a varying relative density with a moderate level for the fixed site index by species. SI: site index. RD: relative density.
In the results presented in Figure 4a, the stand volume increased the fastest in the Japanese larch stand, but it later yielded results similar to those of the Korean white pine stand. Also, the stand volume was the lowest in the Korean red pine stand. In the result shown in Figure 4b, the stand volume was generally highest in the Japanese larch stand at an early stage, but the Korean white pine stand attained a similar stand volume after a stand age of 60 years was reached. The lowest stand volume was recorded in the Korean red pine stand. In Figure 4a, the minimum and maximum stand volume ranged between 210-332 m 3 ha −1 for the Korean red pine, 316-438 m 3 ha −1 for the Korean white pine, and 311-458 m 3 ha −1 for the Japanese larch at 50 years of age. Similarly, in Figure 4b, the stand volume ranged between 247-394 m 3 ha −1 , 301-499 m 3 ha −1 , and 320-548 m 3 ha −1 by species. This comparison indicates a greater difference in the stand volume derived from the stand density.
The variable-density yield models in this study were demonstrated in line with similar prior studies in Korea, where the stand yield was predicted by the stand age and the site index. Using a type of empirical yield model, the previous studies did not encompass the range of variables for the stand density. This was because the samples were taken from the empirical stand where silvicultural treatments were generally executed with the stand density index being relatively low [37,38,49]. Therefore, the comparison in this simulation was carried out by applying a similar stand density to the previous studies. After checking the stand conditions, the mid or mid-upper SI with 0.5 of RD was applied to demonstrate SI 16 for the Korean red pine, SI 18 for the Korean white pine, and SI 22 for the Japanese larch. As a result, the model demonstration for this study was similar to the other yield estimations from previous studies ( Figure 5). The demonstration in this study for the Korean red pine and the Japanese larch was plotted within the predicted range of the preceding studies up to the mature stage of the stand development, while the yield prediction for the Korean white pine was higher in this study than in prior studies. site index. Using a type of empirical yield model, the previous studies did not encompass the range of variables for the stand density. This was because the samples were taken from the empirical stand where silvicultural treatments were generally executed with the stand density index being relatively low [37,38,49]. Therefore, the comparison in this simulation was carried out by applying a similar stand density to the previous studies. After checking the stand conditions, the mid or mid-upper SI with 0.5 of RD was applied to demonstrate SI 16 for the Korean red pine, SI 18 for the Korean white pine, and SI 22 for the Japanese larch. As a result, the model demonstration for this study was similar to the other yield estimations from previous studies ( Figure 5). The demonstration in this study for the Korean red pine and the Japanese larch was plotted within the predicted range of the preceding studies up to the mature stage of the stand development, while the yield prediction for the Korean white pine was higher in this study than in prior studies. Figure 5. Demonstration of the stand volume and simulation comparison with previous studies [37,38,49]. The predicted stand volume for this study was calculated using Equation (9) Demonstration of the stand volume and simulation comparison with previous studies [37,38,49]. The predicted stand volume for this study was calculated using Equation (9) with yield model type 3 based on the fixed-effect parameters in Table 4. To compare the simulation, the variables in previous studies were set to a moderate level, because the model by Yoon et al. [49] was based on the stand age, and the models by KFRI [37] and

Assessment of the Site Index Models
In many countries, the site index has been developed with the appropriate base age of the target species considering not only the dominant height growth characteristics and the modeling data, but also the final felling age [31][32][33][34]39,54,55,[57][58][59][60]. In the 1960s and 1970s, the base age of 20 years was applied by reflecting the stand condition to estimate the current growth and yield in Korea. Since then, the Korea Forest Research Institute has updated the site index with a base age of 30 years when developing the yield table in the 2000s [35][36][37][38][39]. The base age applied in the current study was 10 years higher than in previous studies by the Korea Forest Research Institute and National Institute of Forest Science, but several studies have reported a higher base age than 30 years. In terms of the Korean red pine, for example, the base age was applied as being 40 years by Lee [54] and 60 years by Lee et al. [59]. In terms of the Korean white pine, Lee et al. [32,55] reported site index models with a base age of 40 years in the region of Gangwon and the North Gyeongsang provinces. In terms of the Japanese larch, the base age for the site index models in the region of the Jeolla provinces were applied with a 40 year base age by Jeon et al. [58] and with a 50 year base age by Kim et al. [31].
This was considered suitable taking into account the final felling age currently being implemented in Korea. The final felling age is 60, 70, and 50 for the Korean red pine, the Korean white pine, and the Japanese larch, respectively, in national forests and 40, 50, and 30, respectively, by species, in private forests [64]. The average age by species was between 40 and 50 according to our modeling data (Table 1). Our data also presented continuing vigorous dominant height growth at age 30, and it remained relatively stable after the age of at least 40-50 years ( Figure 3). Therefore, considering these perspectives, our results support the concept of a base age of 40 years in this study.
The site index models developed in this study were featured with data based on observations from the repeatedly measured plot inventory but not from the felled trees for the stem analysis. Given that this type of dataset is known as being more descriptive and suitable for modeling, the site index curves that were developed were evaluated to provide an improved estimate by species [33,42]. The proposed site indices for the Korean red pine, the Korean white pine, and the Japanese larch ranged between 12-18, 12-16, and 16-24 at the Korea Forest Research Institute [37], 12-20, 10-18, and 16-26 at the Korea Forest Research Institute [65], and 8-14, 12-16, and 16-22 at the National Institute of Forest Science [38], respectively, by species with a base age of 30 years. The site indices used in this study were predominantly proposed for 14-20, 16-22, and 22-30 with a base age of 40 years, as explained in our results. The range difference in the major site indices at 2-6 m by species were considered acceptable by taking into account the difference in the base age and, accordingly, the amount of dominant height growth between prior studies and the current study.
When compared with previous studies, there were some similarities and differences in the range of site indices by the species. For the Korean red pine, for example, several studies have proposed a site index range in line with the region and the base age: 6-14 and 8-12 with a base age of 30, respectively, by Pyo [60] and Pyo et al. [66] in the region of the Chungcheong provinces, 6-18 with a base age of 30 by Park and Lee [34] in the region of the Jeolla provinces, 8-18 with a base age of 30 by Yoon et al. [39] and 6-20 with a base age of 60 by Lee et al. [59] in the region of the Gyeongsang provinces. Lee et al. [32] studied the site index with a base age of 40 and proposed a site index range of 14-22 and 18-30, respectively, for the Korean white pine and the Japanese larch in the region of Gangwon and North Gyeongsang provinces. A site index for the Japanese larch was suggested for the region: 12-22 with a base age of 40 by Jeon et al. [58] and 14-22 with a base age of 50 by Kim et al. [31] in the region of the Jeolla province. The site index was analyzed differently by various studies, which supported a relatively higher site index range of this study than that of the previous studies. Considering these results, the application of the site index curves developed in this study is more suitable for Gangwon and the North Gyeongsang provinces where the field data were collected.
The majority of site indices in our observations indicated that the Korean white pine stands grew within a narrow, compact range of the site index compared with the other two species (Figure 3). Based on this finding, in the case where a stand is managed at a suitable site by species, it can be inferred that the stand volume of the Korean white pine stands may be less volatile and less strongly influenced by the site index than the other two species. However, the minority of site indices represented were within a certain range at an extremely low level, especially in the case of the Japanese larch (Table 1, Figure 3), and those could not be regarded as being within the standard site index, even in other studies [31,34]. In this context, it would not be practical and economical in a point of view from trees suitable for a site. Therefore, plots with a poor site index were not suitable for use in plantations, and they would be considered an inappropriate site in terms of recommendations for silviculture. Consequently, the overall site index models were analyzed with significant parameters, an appropriate base age, and site index curves for each species.

Practicability and Applicability of the Variable-Density Yield Models
The parameter estimates generally presented logical behaviors as the predicted stand volumes were higher with increasing AGE, SI, and stand density (Table 4). In contrast to the normal or empirical yield tables, which contain only two variables, namely the stand age and site index, in areas with over-(normal yield table) or fully-stocked levels of density (empirical yield table), these variable-density yield models make it simple to estimate the stand volume in a variety of stand density levels. In a designated stand with a fixed site index, a yield model type 1 with BA or type 2 with TPH can be used to predict the stand volume at a specific time by setting the BA or the TPH to the values required. However, as the stand development is progressing, the BA increases, and the TPH decreases due to the carrying capacity in the stand [43,55], making it difficult to simulate the stand volume. In this case, by applying the yield model type 3 with the RD, the prediction of the stand volume can be easily calculated so that one can use it for their simulation without considering the absolute stand density such as the BA or the TPH, as demonstrated in Figures 4 and 5.
In this context, a stand can be simulated over stand age with a specific RD by varying the site index range across the species (Figure 4). In the simulation example, it can be inferred that stand density affects the stand volume, more than the site index, according to the situation (Figure 4). This example highlights the necessity of variable-density yield models, given that the stand volume can vary depending on the level of the stand density. To easily illustrate our models with other studies, 0.5 of the RD was applied in the simulation for each species under a pre-determined site index ( Figure 5), because other studies predominantly used data from empirical stands as part of research that has been performed alongside silvicultural treatments such as thinning [37,38]. In this simulation comparison, the estimated stand volumes for this study were similar to the other studies. Considering that our estimation resulted in relatively central values, in comparison with the predictions from other studies, our data also represented the empirical forest condition effectively for each species in Korea. In the example, the estimated stand volume for the Korean white pine was the largest in the simulation result, with the same RD by species. This is supported by previous studies [43,55], which have reported that the Korean white pine can retain a higher maximum stand density, compared with the Korean red pine and the Japanese larch.
There were several trials on these studies to develop yield models targeting these species in Korea [5,[46][47][48][49]. However, models from previous studies were developed based on stand conditions, meaning that sufficient data for the mature stands were not collected at that time. The models developed in these studies may be targeted for use in the context of immature stands or different development of stands from the present conditions. Therefore, renewed, updated models are vital for better prediction. Empirical yield tables have also been provided in Korea [37,38]. However, these types of stand yield tables only provide the fixed value of the stand volume in line with the stand age and the site index, regardless of the level of stand density. These were developed based on empirical stands, using a base age of 30 years for the site index for all species, where the relative density was prone to being solely moderate, which makes it difficult to represent a variety of stand densities. The models were developed and updated using the data from intensive experiments and from the National Forest Inventory, which may include other species in the plots used, as well as in plantations with the target species. Some of the impure stands, including unevenaged or mixed forests, could be considered to be included depending on the definition of pure stands. This is because the species in a stand is defined as a dominant species when its basal area is greater than 50% of the total in a stand [67,68]. In these contexts, yield models of pure stands should be examined in more detail in Korea, especially for the main commercial species.
To resolve the limitations and disadvantages in previous studies and to update the models, our models were evaluated to provide better predictions with flexibility and practicability. The yield models developed helped to estimate the stand volume with stand density variables according to the absolute stand density (BA and TPH) or the relative stand density (RD). Given that prior research proposed the base age by data characteristics from their studies, the suggested base age of 40 years in this study can be an alternative solution considering the current condition of stand development in Korea based on examination of the site index in this study [31,32,39,54,59]. The data used in this study represented the pure stands better than the data from the national forest inventory. It should be more accurate and predictable for the monoculture of the Korean red pine, the Korean white pine, and the Japanese larch in Korea. Therefore, the variable-density yield models were highly applicable in predicting the stand volume as whole-stand models for the Korean red pine, the Korean white pine, and Japanese larch stands in Korea.
Overall, the maximum and minimum range of the stand volume estimations, were covered by or similar to the range of preceding studies ( Figure 5). Given that a wide range of stand densities were used in our study, and the stand density variable by model type was added to the models developed, it was evaluated to be more practical for yield prediction by stand density. Accordingly, this can be extended to an interactive application with biomass and forest carbon conversion [37,38,46,48]. It was not feasible to provide tree volumes by diameter class or the estimation according to the thinning treatment because the developed models in this study relied on stand-level data without thinning predictors. To address these types of detailed predictions, more specific model types should be further developed in future studies with tree-level data and/or repeated measure thinning experiments [69].
The models developed in this study were available for stand yield prediction from young to mature stands using three types of stand density variables (Table 4). However, it is important to apply caution in predictions for stands younger than 15 years because the age range was not included in the modeling due to an insufficient number of samples and a biased fit. The monoculture of the Korean red pine, the Korean white pine, and the Japanese larch was targeted, especially those grown in the region of Gangwon and the North Gyeongsang provinces, contrary to previous studies [37,38,49] which targeted the whole region of Korea. The different amounts of stand volume were substantiated in line with the specific regions and their data between our study and previous studies [37,38,49]. Therefore, the models developed in this study will be more accurate when applied within the stands of these targeted regions or provinces.

Conclusions
To provide the stand volume for the Korean red pine, the Korean white pine, and the Japanese larch, variable-density yield models with the necessary site index models were developed using linear and nonlinear mixed-effect models. The empirical plot data were successfully applied for both the yield and site index models. The site index models with the Chapman-Richards function for each species were suitably developed with a base age of 40 years. It was considered to be more feasible to use a greater base age than the current one of 30 years, by taking into account the final felling age, the dominant height growth characteristics, the age range in the dataset, and the more mature forests in Korea. The site index curves developed in this study were considered to effectively represent suitable trees at the site. The major range of the site index was 14-20 for the Korean red pine, 16-22 for the Korean white pine, and 22-30 for the Japanese larch. These site indices generally indicated superior growth characteristics in the study region compared to the entire country.
The variable-density yield models for each species fitted reasonably with the types of stand density variables, BA, TPH, and RD, with AGE and SI. All the models were unbiased as to the predicted values and all the independent variables. Given that the models developed encompassed a wide range of stand densities, they are applicable to any stocking level from normal stands to empirical stands regardless of the thinning intensity. The yield model type 3 with RD is flexible enough to simulate a development of stands by setting a relative density instead of the absolute density, such as BA and TPH. It is important to notice that the modeling data originated from stands older than or equal to age 15 and less than or equal to RD 1.3 in Gangwon and the North Gyeongsang provinces of Korea, suggesting that extrapolation should be carefully handled considering the temporal and spatial range of the data. Consequently, all the models were considered to be highly practicable for the stands of Korean red pine, Korean white pine, and Japanese larch, especially in Korea.