Stand-Level Components of a Growth and Yield Model for Nothofagus Mixed Forests from Southern Chile

: Reliable information on stand dynamics and development is needed to improve management decisions on mixed forests, and essential tools for this purpose are forest growth and yield (G&Y) models. In this study, stand-level G&Y models were built for cohorts within the natural mixed second-growth Nothofagus -dominated forests in Chile. All currently available (but limited) data, consisting of a series of stratiﬁed temporary and permanent plots established in the complete range of this forest type, were used to ﬁt and validate these models. Linear and nonlinear models were considered, where dominant stand age, number of trees, and the proportion of basal area of Nothofagus species resulted in signiﬁcant predictors to project future values of stand basal area for the di ﬀ erent cohorts (with R 2 > 0.51 for the validation datasets). Mortality was successfully modeled ( R 2 = 0.79), based on a small set of permanent plots, using the concept of self-thinning with a proposed model deﬁned by the idea that, as stands get closer to a maximum density, they experience higher levels of mortality. The evaluation of these models indicated that they adequately represent the current understanding of dynamics of basal area and mortality of Nothofagus and companion species in these forests. These are the ﬁrst models ﬁtted over a large geographical area that consider the dynamics of these mixed forests. It is suggested that the proposed models should constitute the main components of future implementations of G&Y model systems.


Introduction
Computational tools such as growth and yield (G&Y) models can be used by forest professionals to plan and implement management strategies at the local or regional scale. For instance, knowledge of growth dynamics, such as basal area growth and mortality, can be combined with forest inventories to determine timber production and examine the potential impacts of alternative management and harvesting regimes on the value and sustainability of the forest [1].
Growth and yield (G&Y) models for mixed forests were first developed during the 20th century, more than 150 years after such models were first implemented on commercial tree plantations [2]. Most models for mixed forest stands are limited and mainly focused on North American and European forests. In contrast to mono-specific plantations, mixed forests present additional complexity, Forests 2020, 11, 810 3 of 15 of currently available, but still limited, information from temporary and permanent plots established in the complete range of this forest type. The specific objectives included to fit and validate stand-level models for (1) basal area according to the cohorts of Nothofagus and companion species, (2) changes over time in the proportion of Nothofagus trees in a stand, and (3) tree mortality that considers the concept of self-thinning. These models should constitute the main component of future implementations of G&Y systems for this resource.

Data Description
The data for this study originated from three independent sets of sampling series, i.e., two temporary plot (TP1 and TP2) and one permanent plot (PP) series. All plots were established in second growth RORACO forests in Chile, and they are located between the 36 • and 42 • S latitude. The TP1 and TP2 series were established by the Universidad Austral de Chile between 1999 and 2000 (see [15] for original sampling methodology). The TP1 and TP2 series were sampled according to a stratification representative of the RORACO forest type based on the latest national forest inventory [11]. The TP1 data had a total of 50 plots with an area of 250 m 2 formed by a conglomerate of two subplots. Meanwhile, the TP2 data had a total of 120 rectangular plots with areas ranging between 250 and 500 m 2 . In contrast, the PP series consisted of three sites under silvicultural thinning measured between 1980 and 1999 with each site remeasured up to four times [25]. However, for the latter, only the subset of plots without treatment (i.e., controls), low thinning (less than 5% of removed basal area), and girdling treatments were considered in the present study. Hence, the PP series dataset consisted of 48 plots with 183 unique measurements, providing several growth and mortality periods. A map of the locations of plots of these three datasets are presented in Figure 1.
For all plots, trees above 5 cm of diameter at breast height (DBH, cm) were inventoried for DBH and total height (H, m). The Nothofagus species were identified and the rest was recorded as companion species. For all plots, the following stand-level variables were calculated: dominant age at breast height (AGE, years), dominant height (HD, m), site index (SI, m), total basal area (BA, m 2 ha −1 ), and total number of trees (NHA, trees ha −1 ). Quadratic diameter (DQ, cm) was measured and defined as the diameter of the tree with average basal area. Dominant age at breast height (AGE) is defined as the average age of 100 trees per hectare with the largest DBH. Dominant height (HD) is the average total height of the thickest 100 trees per hectare. Site index (SI) is the stand dominant height at 20 years (see [26] for further details). Additionally, for each of the cohorts, basal area for Nothofagus and companion species (BAN and BAC, respectively, m 2 ha −1 ), and number of trees of Nothofagus and companion species (NHAN and NHAC, respectively, trees ha −1 ) were calculated. Finally, the proportion of basal area and number of trees of Nothofagus (PBAN, PNHAN) and companion species (PBAC, PNHAC) were also obtained. To study the effect of geographic location in growth and productivity, all plots were assigned to their growth zone (ZONE) (further details on zoning are presented in [26]).
In order to focus on stands that were dominated by Nothofagus species, only those plots with PBAN > 0.6 were selected for this study. Summary statistics of these plots are presented in Table 1. Additionally, the dominant species (DOM-SP) of a given plot was defined as the Nothofagus species that had more than 70% of BA. In terms of composition, the TP1 and TP2 data are primarily of N. dombeyi but all dominant species are present; however, the PP data contains only plots dominated by N. alpina (Table 2).     To predict basal area for the two cohorts, BAN and BAC, this study fitted two independent models. All three Nothofagus species were combined into the same cohort based on reported low levels of differentiation between the species [26]. Given the limitation of the data, the TP1 and TP2 plots were used as training data, whereas only PP plots were used as validation data. Hence, validation is focused on the performance of N. alpina, which is a limitation. As more long-term data become available, further validations are possible. For BAN and BAC, linear models were fitted using a log transformation of these responses with different combinations of predictors, including AGE, HD, SI, NHA, NHAN, NHAC, PBAN, and PBAC. These predictors were considered in their original units or transformed using the functions of natural logarithm, inverse, square of the inverse, and square root of the inverse. ZONE for each stand was evaluated separately with no interactions with other predictors. To assist with model selection, a backward selection procedure was implemented based on a significance level set to α = 0.05, and models with variance inflation factors (VIFs) larger than four in any of their predictors were discarded. It is assumed that the BAN and BAC form the total stand basal area; hence, BA = BAN + BAC.
To estimate the future values of BAN and BAC given a starting condition, the prediction equations fitted above were derived into projection models by differentiating with respect to age following the methodology described by Clutter [27]. This study defined projection models as those that use current stand conditions to project parameter values into the future; hence, this methodology converted yield Forests 2020, 11, 810 6 of 15 models into compatible growth equations using the same model parameters. Simultaneous fitting of prediction and projection equations was not possible given that the TP1 and TP2 datasets contained only temporary plots, where projection models require basal area growth data measured over intervals. To evaluate these models, the PP series data, consisting of 217 growth intervals, were used, but note that this dataset, as indicated earlier, is dominated by N. alpina stands.

Proportion of Nothofagus Trees
To predict the proportion of trees corresponding to the Nothofagus cohort (PNHAN), a linear model was fitted based on the logit transformation of PNHAN using, as similar to the basal area fit, the TP1 and TP2 plots as training data and the permanent plot measurements from PP as validation data. The same predictors (and their transformations) used in the BA model were tested and a final model was selected using a backward selection procedure as indicated above. As stated before, all Nothofagus species were combined into a single cohort given their similarity in biological behavior.

Mortality
Given the limitations on long-term mortality data that adequately represent all dominant species, we followed a two-step procedure in this study. First, coefficients from Reineke's self-thinning models [28] were first obtained. These were then used to fit a proposed nonlinear model defined by the concept that, as stands get closer to a maximum density, they experience higher levels of mortality.
In a previously published study using the information from TP1 and TP2 [24], the following self-thinning expression was fitted [28]: where ln() is the natural logarithm and α and β are the parameters to estimate; NHA and DQ were described previously. The above model was obtained to stands dominated by N. alpina, N. obliqua, or N. dombeyi, separately, and parameters are reported by [24]. However, a single slope (β = −1.41) was found for all three species but with different intercepts (11.61, 11.37, and 11.76, respectively). The above model can be used, together with the current density value (NHA 0 ), to estimate the current maximum (or limiting) quadratic diameter (DQ 0max ), as follows: Here, DQ 0max is interpreted as the maximum stand quadratic diameter that is allowed at a density of NHA 0 . Hence, the model from Equation (2) implies that, as the current quadratic diameter (DQ 0 ) approaches DQ 0max , there is an increase in mortality and that stands dominated by the same species respond to the self-thinning rule evenly.
For the second step, the 217 growth intervals from the PP series were used to fit the following proposed nonlinear mortality model that uses both the current conditions (NHA 0 , DQ 0 and its DQ 0max ), as follows: where θ is the parameter to estimate and can be interpreted as a maximum mortality rate when the stand is at DQ 0max , and ∆t is the years between growth intervals.

Model Evaluation
Predictions and projections for the model of basal area of Nothofagus and companion species, proportion of Nothofagus trees, and mortality were evaluated by calculating the following goodness-of-fit measures: R 2 emp , RMSE%, and Bias%, which are detailed below. These measures were obtained for the training and validation datasets providing two assessments of the models.
where y i andŷ i are the ith observed and predicted (or projected) value, respectively;ȳ is the mean response observed value; and n is the number of observations. All linear and nonlinear models were fitted using the statistical software R version 3.3.2 [29]. To detect multicollinearity between predictors (see Table 1), variance inflation factors (VIFs) were checked. All goodness-of-fit statistics were evaluated using the back-transformed response variables to their original units. Because the models for BAN and BAC use the natural logarithm transformation, their back-transformed estimates were adjusted using the correction, i.e.,ŷ i * =ŷ i exp (σ 2 /2), where σ 2 is the mean square error. For graphical outputs, relative residuals were used, which were defined as the difference between observed and predicted values divided by the mean observed value. As an additional check, projection models were evaluated by using all 217 possible growth intervals within the PP data; here, the earlier plot measurement was used as the initial conditions to perform the simulations, and these were projected over the growth interval. Then, the later measurement was used to contrast observed against predicted values. For this dataset, the time between growth intervals ranged between 2 and 12 years.

Basal Area
For all plots considered in this study, the average total BA for Nothofagus and companion species corresponded to 38.48 and 3.41 m 2 ha −1 , respectively. BAN ranged from 12.66 to 89.57 m 2 ha −1 , and BAC from 0.00 to 26.40 m 2 ha −1 .
The final selected models for BA of Nothofagus and companion species were as follows: ln(BAN) =β 0 +β 1 ln(AGE) +β 2 ln(SI) +β 3 ln(NHA) +β 4 ln(PBAN) ln(BAC) =β 0 +β 1 ln(AGE) +β 2 ln(PNHAN) +β 3 ln(PBAN) The logarithmic transformation of the predictors returned the best results and had the additional advantage that it derives the projection models easily (see below). For the selected models, all final predictors showed low VIF values (<2.1), reflecting negligible levels of multicollinearity between them. The resulting fitted model for BAN had R 2 emp = 0.54, and the fitted model for BAC had a higher R 2 emp with a value of 0.85 (Table 3). The prediction of total basal area had an R 2 emp = 0.56. This moderate correlation reflects the wide range of conditions found in these forests. All models presented negligible bias (<1%). Both BAN and BAC models had good goodness-of-fit measures with the PP validation data (Table 3), where the BAN, BAC, and BA predictions returned slightly unfavorable higher Bias% values when compared to the training data, but these were all lower than 4%.
According to the estimated coefficients (Table 4), AGE was positively associated with both BAN and BAC (with slope coefficients of 1.21 and 0.09, respectively). Hence, as the stand gets older basal area increases, with a larger effect for the Nothofagus cohort. For BAN, the positive coefficients for SI (0.65) and NHA (0.52) indicate that better site quality and higher levels of stocking result in higher Nothofagus basal area. In the BAC model, PNHAN and PBAN have negative coefficients (−0.22 and −1.87, respectively) indicating that higher proportions of Nothofagus species abundance affect the amount of basal area of companion species. Table 3. Goodness-of-fit measures for models for the basal area of Nothofagus species (BAN, Equations (7) and (9)), basal area of companion species (BAC, Equations (8) and (10)), total basal area (BA), and the proportion of number of Nothofagus species trees (PNHAN, Equation (11)).   (7) and (9)), basal area of companion species (BAC, Equations (8) and (10)), and the proportion of number of Nothofagus species trees (PNHAN, Equation (11)). All model parameters were found to be significant (p < 0.001).  (Figures 2 and 3). However, this correspondence decreased with larger observed BAN and BAC values, and some under-prediction was found for observed BAN values above 75 m 2 ha −1 . Similar results were found for BA, as this mostly corresponds to Nothofagus species basal area ( Figure 2C). Normality and heterogeneity of residuals were also checked without important departures from these assumptions. Both basal area equations were used to derive their compatible projection equations. These models project future values (BAN 1 and BAC 1 ) based on the current stand conditions (BAN 0 and BAC 0 , respectively). These are as follows:

Model
where the β coefficients are the same parameters from the fitted Equations (8) and (9). For the evaluation of the projection equations using the validation dataset, all basal area models showed excellent goodness-of-fit measures (all with R 2 emp > 0.94). The relative residuals obtained over time for BAN, BAC, and BA projections ( Figure 4A-C) are centered around zero for shorter projections (i.e., low bias), whereas they tend to depart for increasing projection times (i.e., under-estimate). For projections under 6 years, the BAN model returned relative residuals lower than 10% and were centered about zero. After 6 years, the relative residuals reached higher values with a tendency to under-predict basal area. The BAC model had residuals centered around zero with no notable deviations even at 12 years of projections. However, there were some projections with residuals over 30%, which are not of relevant concern because of the low proportion of basal area from the companion cohort in the sampled plots. Predicted BAN and BAC values corresponded well with observed values in both training and validation data (Figures 2 and 3). However, this correspondence decreased with larger observed BAN and BAC values, and some under-prediction was found for observed BAN values above 75 m 2 ha −1 . Similar results were found for BA, as this mostly corresponds to Nothofagus species basal area ( Figure  2C). Normality and heterogeneity of residuals were also checked without important departures from these assumptions.
Both basal area equations were used to derive their compatible projection equations. These models project future values (BAN1 and BAC1) based on the current stand conditions (BAN0 and BAC0, respectively). These are as follows: where the β coefficients are the same parameters from the fitted Equations (8) and (9).   (3)). All panels are from the Temporary Plot series (TP1 and TP2).
For the plots considered in this study, the average proportion of Nothofagus species trees corresponded to 82%, where most of them presented values greater than 72%. The final selected model for PNHAN was the following:  For the training data, this model had reasonable goodness-of-fit measures with R 2 emp = 0.68 and Bias% = −1.50. Additionally, predicted PNHAN values tended to correspond with observed values, but high levels of uncertainty still existed ( Figure 2D). Moreover, for PP validation data, these measures were R 2 emp = 0.56 and Bias% = 3.46. The predictions for PNHAN tended to have less uncertainty with higher observed PNHAN, as observed in Figure 3D.
The estimated parameters of this model are shown in Table 4. The slope coefficient for PBAN (10.29) reflects the high association between this predictor and PNHAN (these predictors present a correlation of 0.89). For AGE, its coefficient (−0.01) indicates a reduction of PNHAN with increasing stand age, reflecting the pioneer behavior of Nothofagus species and the gradual establishment of companion species over time. These selected predictors all show low VIF values (<1.02).

Mortality
Among the remeasured plots from the PP data, the observed annual mortality rates had an average of 3.0% with a maximum of 14.2%. Their patterns were consistent over time for most plots evidenced from the parallel trajectories, as shown in Figure 5A.   (3)) using the PP data as validation.

Proportion of Nothofagus Trees
For the plots considered in this study, the average proportion of Nothofagus species trees corresponded to 82%, where most of them presented values greater than 72%. The final selected model for PNHAN was the following: For the training data, this model had reasonable goodness-of-fit measures with R 2 emp = 0.68 and Bias% = −1.50. Additionally, predicted PNHAN values tended to correspond with observed values, but high levels of uncertainty still existed ( Figure 2D). Moreover, for PP validation data, these measures were R 2 emp = 0.56 and Bias% = 3.46. The predictions for PNHAN tended to have less uncertainty with higher observed PNHAN, as observed in Figure 3D.
The estimated parameters of this model are shown in Table 4. The slope coefficient for PBAN (10.29) reflects the high association between this predictor and PNHAN (these predictors present a correlation of 0.89). For AGE, its coefficient (−0.01) indicates a reduction of PNHAN with increasing stand age, reflecting the pioneer behavior of Nothofagus species and the gradual establishment of companion species over time. These selected predictors all show low VIF values (<1.02).

Mortality
Among the remeasured plots from the PP data, the observed annual mortality rates had an average of 3.0% with a maximum of 14.2%. Their patterns were consistent over time for most plots evidenced from the parallel trajectories, as shown in Figure 5A. The fitted nonlinear model had a good fit with R 2 emp = 0.79, RMSE% = 18.46 and Bias% = −2.80. The single parameter estimate corresponded to θ = 0.003595746 (SE = 0.000213), indicating that, for future projections, the estimated number of trees will always be smaller than the current condition. The predicted mortality values had good correspondence with observed ones over the entire range of the data ( Figure 5B) and with relative residuals ranging from −30% to 30% ( Figure 4D). While these residuals may be considered as large model uncertainty, the fact that residuals are generally centered around zero, even after 12 years of projection, suggests a good overall accuracy. The fitted nonlinear model had a good fit with R 2 emp = 0.79, RMSE% = 18.46 and Bias% = −2.80. The single parameter estimate corresponded to θ = 0.003595746 (SE = 0.000213), indicating that, for future projections, the estimated number of trees will always be smaller than the current condition. The predicted mortality values had good correspondence with observed ones over the entire range of the data ( Figure 5B) and with relative residuals ranging from −30% to 30% ( Figure 4D). While these residuals may be considered as large model uncertainty, the fact that residuals are generally centered around zero, even after 12 years of projection, suggests a good overall accuracy.

Discussion
The fitted model for BAN presented here is biologically realistic as it accounts for stand age, productivity, and stocking using AGE, SI, and NHA, respectively. These are common predictors found in other reported Nothofagus species growth models [30]. In the present study, the effect of productive growth geographic zones (ZONE) on basal area growth was not significant, which is likely to be the result of the incorporation of site index (SI) in the model, an in situ stand parameter that describes that productivity. However, other studies have found that growth zone is a significant component [4,24,[30][31][32]. Additional future data should allow exploring these responses further together with the evaluation of effects of other environmental and biological factors that may affect Nothofagus species growth, such as light conditions, elevation, and nitrogen availability [22,33,34]. The fitted BAN models resulted in moderate correlations for the training and validation datasets (<0.54). This reflects the level of heterogeneity found on these forests, with many aspects that might influence the accuracy of this model, such as uncertainty on determining age, anthropic alterations, variable climate at plant establishment, etc.
The model for the proportion of Nothofagus trees reflects the pioneer behavior of the Nothofagus species [7], evidenced by the negative coefficient related to dominant age. It also shows the increasing establishment of the companion (shade-tolerant) species as the stand gets older, reflecting some gradual forest succession.
The self-thinning rule, which has been mostly applied in pure stands with some examples in mixed forest stands [24], worked successfully here to predict mortality. This natural self-thinning is occurring in stands as young as 25 years of dominant age, 15 years earlier than shown in stands dominated by N. obliqua [8]. Additional long-term measurements of permanent plots should allow the use of more widely tested stand-level mortality models, such as the ones presented by Thapa and Burkhart [35].
Based on the mortality and basal area projection models, it is possible to propose a flow that builds a simple simulation system of equations that can realistically project BAN values into the future for RORACO forests. To illustrate this system, simulated stands dominated only by Nothofagus species (i.e., PBAN equal to 1), with an initial BAN of 15 m 2 ha −1 , stand age of 15 years, SI of 10 m, and several initial NHAN values are presented in Figure 6. These simulated stands showed that basal area patterns grow asymptotically with larger BAN growth rates for forests with fewer initial trees ( Figure 5A). Furthermore, as expected, higher initial tree densities resulted in higher rates of mortality ( Figure 6B) following the patterns described by the proposed mortality model.  Figure 6. These simulated stands showed that basal area patterns grow asymptotically with larger BAN growth rates for forests with fewer initial trees ( Figure 5A). Furthermore, as expected, higher initial tree densities resulted in higher rates of mortality ( Figure 6B) following the patterns described by the proposed mortality model. In this study, the basal area models for the two cohorts (Nothofagus and companion species) were assumed to be independent. This assumption can be challenged by studies that suggest that additive effects allow higher yields of Nothofagus species when companion species are present [36,37]. However, the choice of independent growth models is supported by research that showed a lack of correlation between basal area of Nothofagus and companion species [38]. Furthermore, the sum of our proposed models accurately estimates total basal area (BA), but further measurements should allow for a better evaluation of these hypotheses.
The use of cohorts for G&Y models presents additional flexibility, and it is a reasonable In this study, the basal area models for the two cohorts (Nothofagus and companion species) were assumed to be independent. This assumption can be challenged by studies that suggest that additive effects allow higher yields of Nothofagus species when companion species are present [36,37]. However, the choice of independent growth models is supported by research that showed a lack of correlation between basal area of Nothofagus and companion species [38]. Furthermore, the sum of our proposed models accurately estimates total basal area (BA), but further measurements should allow for a better evaluation of these hypotheses.
The use of cohorts for G&Y models presents additional flexibility, and it is a reasonable compromise between whole-stand models and single-tree models [1,6], particularly for mixed forests, providing some additional granularity on the construction of these models. The similarity on growth behaviors of N. alpina and N. dombeyi [31] and the reliability of stand-level projections presented in this study justified the grouping of these species into a single cohort. In addition, there is evidence of low differentiation among these species, as reported in a multivariate analysis that used the same dataset from the current study [26] and from other authors indicating that stands dominated by either N. alpina or N. obliqua do not greatly differ in the total basal area or canopy height [30]. While it could be useful to have individual stand models for each of the Nothofagus species for further granularity, the present limitations of the dataset (which has not adequately represented each species) is likely to affect the accuracy of the final fitted models.
Unfortunately, the permanent plot (PP) series is, at present, our only currently available source of remeasured data to fit and validate the proposed Nothofagus species models, and its main limitations is that it is dominated by N. alpina. Hence, further improvements to our proposed models can only be achieved with additional permanent plots that better represent the other Nothofagus species, and that span over a wider geographical range. This lack of time series is particularly relevant for mortality that, in this case, was based on Reineke's expression and not on repeated measurements of the same plot.
Regardless of these limitations, these models represent critical components of growth and yield models for the natural mixed forests of the RORACO forest type in Chile, and importantly they are using information from stratified temporary and permanent plots established in its complete geographical range.

Conclusions
For this study, several stand-level models were built to improve the predictability of stand dynamics for natural mixed secondary forests of the RORACO forest type in Chile. Stand mortality was successfully modeled with a function defined by the concept that, as stands get closer to a maximum density, they experience higher levels of mortality. To our knowledge, these are the first broadly applicable growth and yield models for the RORACO forest type with dynamics of both companion species and Nothofagus cohorts. These models can be incorporated as part of a system of modules for a growth and yield simulator of this forest type. This system could also include local conditions and model operational thinning. The models reported here constitute simple and valuable tools to support management decisions for this resource in Chile.