Prediction Comparison of Stand Parameters and Two Ecosystem Services through New Growth and Yield Model System for Mixed Nothofagus Forests in Southern Chile

: Forest managers need tools to predict the behavior of forests not only for the main stand parameters, such as basal area and volume, but also for ecosystem services such as timber volume and carbon sequestration. Useful tools to predict these parameters are growth and yield model systems with several possible options for modeling, such as the whole stand-level model, with or without diameter distribution generation, individual tree-level model, and compatibility models. However, those tools are scarce or developed mainly for forest plantations that are mostly located in the northern hemisphere. Thus, this study focuses on analyzing predictions of several growth and yield models built for native mixed Nothofagus forests from southern Chile, using the simulator Nothopack. A dataset of 19 permanent plots with three measurements were used for comparing the different models. Individual tree-level simulation presented the best goodness-of-ﬁt statistics for stand parameters and ecosystem services. For example, the basal area gave an R 2emp of 0.97 and 0.87 at 6 and 12 years of projection. However, the stand-level simulations with a generation of diameter distribution and both compatibility models showed satisfactory performance, both in accuracy and bias control. The simulator Nothopack, which has the capability of obtaining detailed outputs, is a useful tool to support management plans for these forest ecosystems.


Introduction
Growth and yield (G&Y) models help to represent the dynamics of natural and artificial forest stands and they include growth, mortality, and other modules to describe changes in stand structure [1]. Forest G&Y models can be classified into stand-level (low resolution) or individual-level (high resolution) models [2]. Stand-level models are those in which the modeling units are aggregated parameters such as basal area, stocking, and site productivity [3]. In contrast, individual-level models keep track and describe each tree in the stand. Both levels have advantages and disadvantages: stand-level models present well-behaved predictions in the long-term for stand parameters; however, they are inadequate to predict tree variables (such as diameter distributions or individual competition). In contrast, individual-level models are better at predicting the variables of trees and hence key information concerning timber assortments, while often lacking precision when summarizing to stand-level parameters [4].
To exploit the advantages of both model types and to improve predictions, mathematical methods have been developed to link stand-and individual-level models into a compatible system [5][6][7][8]. For instance, predicted individual tree basal areas can be adjusted so that their sum equals a predicted stand-level total basal area [4]; thus, making these predictions

Data Description
Data for this study originated from 19 permanent plots (PP) established in secondgrowth RORACO forests in southern Chile, dominated mainly by N. obliqua and located between the 38 • 30 and 42 • 40 S latitude ( Figure 1). The research project D97I1065 from the Universidad Austral de Chile established the PP network in 2000, with areas of 250 or 500 m 2 (<4800 trees ha −1 ) depending on the stand density. These plots have the following two remeasurements: all plots in 2006, and a subset of six plots in 2012.
Trees above 5 cm in diameter at breast height (DBH) and taller than 2 m were measured for DBH (cm). Total height (H, m) was sampled sequentially for a subset of 15 trees per plot. Nothofagus species were identified (SP, with N. alpina (1); N. obliqua (2); and N. dombeyi (3)), and the rest were recorded as other species (4). Sociological status (SS), according to vertical stratification was also obtained (dominant (1); codominant (2); intermediate (3); and suppressed (4)). Increment cores, at breast height, were extracted from the same subset of trees for total height to determine breast height age (years) at the first measurement.

Simulations
Growth simulations were developed to compare the behavior and uncertainties in stand-parameters and the provision of ecosystem services, such as timber volume and carbon sequestration. For this study, the available simulator Nothopack [19], built for RO-RACO second growth forest, was used. The simulation types considered were stand-level without generation of a diameter distribution (S1), stand-level incorporating a generated diameter distribution (S2), individual-level (T), yield compatibility (YC), and growth compatibility (GC). These two compatibility simulations integrated results from S1 and T (further details of these methods are presented below). For all simulations, growth estimates and predictors were updated annually until projection age. Further details of the models associated with these simulations are presented below. The independent dataset used to assess these simulations corresponded to the remeasurements of the 19 plots from the PP network. The first projected period corresponded to six years between 2000 and 2006 and from 2006 to 2012. A second projected period for analyzing corresponded to 12 years from 2000 to 2012.

Stand-Level without Diameter Distribution Simulation
The Nothopack's stand-level simulation has a growth module and a mortality module. The first performs the estimation of the basal area annual growth (BA), specifically, the sum of basal areas of N. alpina, N. obliqua, and N. dombeyi (BAN), and the basal areas of other companion species (BAC). The second module incorporates a modification from Reineke's self-thinning models to estimate future stand density. Further details are presented in Palmas et al. [11].
The compatible projection models to predict futures values (BAN 1 and BAC 1 ) based on current conditions (BAN 0 and BAC 0 ) are as follows: where theβ coefficients are the same parameter estimates from Equations (1) and (2). PBAN was kept constant over the projection time (PBAN 1 = PBAN 0 ), and PNHAN 1 was modified over time, incorporating the difference between year 0 and 1, using the following equation for each year: where theβ coefficients also are shown in the Appendix A and other terms were explained previously.
The stand-level mortality model is an annual projection of NHA based on the expression from Reineke [20] and defined by the following two parameters:α andk. The original model was presented by Palmas et al. [11] and parameter estimates are shown in the Appendix A, and its expression is as follows: where NHA 0 is the current number of trees per hectare, ∆t is the number of years between growth periods, QD 0 is the current quadratic diameter, and QD 0 max is the current maximum quadratic diameter defined by the following: whereα varies depending on if the stand is dominated by N. alpina (α 1 ), N. obliqua (α 2 ), N. dombeyi (α 3 ), or a mix of them (α 4 ), andβ is a constant for all stands [18]. Total stand volume (VHA, m 3 ha −1 ) was estimated using the following general volume equation proposed by Gezan and Ortega [21] with BA, HD, PNHAN, and PBAN as predictors: The parameter estimates are presented in the Appendix A.

Stand-Level with Diameter Distribution Simulation
This simulation used the same models for basal area growth, and mortality as in Section 2.2.1, but with the generation of the following three-parameter Weibull-based diametric distribution model: where f (x) is the probability evaluated at x, A is the position parameter (set at 5 cm), and B and C are scale and shape parameters, respectively. These were obtained using the parameter recovery method as described by Gezan and Ortega [21], and they differ among species according to the following equations. N. alpina (1): C =β 0 +β 1 B +β 2 QD 1 (11) N. obliqua (2): B =β 0 +β 1 QD 2 +β 2 PNHAN +β 3 ln(RS) C =β 0 +β 1 B +β 2 QD 3 +β 3 QD Companion species (4): where sub-index for stand parameters correspond to values associated with the specific species, relative spacing was calculated using the expression RS = [(10,000/NHA)0.5]/HD * 100, and parameter estimates are shown in the Appendix A.
The above models allow the estimation of BA and NHA for each diameter class and species. Volume was later calculated per diameter class using available parametric height model (H, m) (Equation (18)) and a taper equation (Equation (19)) reported by Gezan et al. [22]. The values of BA, NHA, and volume for DBH classes were aggregated to stand-level.

Tree-Level Simulation
The individual tree simulation for Nothopack has an annual growth tree DBH module, and it uses the same stand-level mortality module as in previous simulations. Volume was calculated by each tree using the taper Equation (19) and then aggregated to stand-level. Here, in each year of the time projection, the DBH was updated by adding the following estimated annual diameter increment (AIDBH, mm year −1 ) based on the model [12]: ln(AIDBH) =β 0 +β 1 ln(BALn + 10) +β 2 SDI +β 3 ln(DBH) +β 4 ln(AGE) +β 5 SS (20) where BALn is the basal area of larger trees for Nothofagus (m 2 ha −1 ) and other terms were described before. All coefficient values can be found in the Appendix A.

Compatibility Simulations
Two compatibility methods were considered to calibrate the projections of the individuallevel models, referred to as proportional yield compatibility (YC) and proportional growth compatibility (GC; also known as disaggregation) [4,23]. Both first calibrate individual mortality probabilities based on the stand-level mortality projection (NHA 1 , from Equation (6)), and then calibrate each tree's DBH to equate the summation of their basal area to the stand-level BA simulation (BA 1 , from Equation (3)). The expressions for the YC method are as follows: where FT 1i and FT 1i are the calibrated and predicted expansion factors (the number of trees that each sample tree represents) for ith tree at time 1, respectively; NHA 1 is the trees per hectare at time 1; DBH 2 1i and DBH 2 1i are the calibrated and predicted squared diameters for the ith tree at time 1, respectively. BA 1 is the predicted basal area at time 1 from Equation (3); and K = π/40, 000, is a constant.
The expressions for the GC method are as follows: where DBH 2 0i is the squared diameter of the ith tree at time 0 and all other terms were previously defined. Note that Equation (23) requires a search algorithm to estimate a power value m that makes the sum of the predicted expansion factors equal to NHA 1 .

Evaluation of Ecosystem Services
The following two ecosystem services were considered based on the growth and yield simulations: provisioning and regulating services. The first is associated with the provision of fiber wood, specifically saw-log (SL), pulpwood (PW), and residual lumber (RES). The second is associated with carbon sequestration, which, in this study, aimed at total tree living biomass. Using taper models, it was possible to estimate different stem products for S2 and T simulations; however, for S1 simulation, this was not possible.
At present, the typical market log products in Chile are [24] (1) SL with a length of 3.7 m and a minimum log diameter of 24 cm, (2) PW with a length of 2.44 m and a minimum diameter of 10 cm, and (3) the remaining individual stem volume corresponding to RES. In the case of carbon sequestration, the carbon stock calculation (C, tonne C), above-ground tree dry matter biomass (AGB, tonne d.m.) was estimated based using specific models for each species previously reported by Milla et al. [25], which are as follows: where model terms were already described, and parameter estimates can be found in the Appendix A.
A root-shoot ratio (R, tonne d.m.) of 0.24 was considered the same for all species using tier 1 for temperate broadleaf forests [26]. Carbon fraction (CF, tonne d.m. −1 ) was specific to species with values of 0.4445, 0.4248, 0.4370, and 0.4282 for N. alpina, N. obliqua, N. dombeyi, and companion species, respectively [27]. The formula considered to calculate carbon stock was as follows: where i = 1, 2, 3, and 4 corresponding to N. alpina, N. obliqua, N. dombeyi, and companion species, respectively. Finally, carbon sequestration over the simulation period was calculated by subtracting the carbon stock at the initial time of simulation from the carbon stock at the end of the projection.
The following goodness-of-fit measures were used to evaluate all the following simulations: empirical coefficient of correlation (R 2 emp ), relative root mean square error (RMSE%), relative bias (BIAS%), and Theil's inequality coefficient (U2), calculated from the following expressions [28,29]: where y i andŷ i are the ith observed and projected value, respectively; y is the mean response of the observed value, and n is the number of observations. For graphical outputs comparing scenarios, relative residuals were used, which were defined as the difference between observed and predicted values divided by the mean observed value and represented as a percentage. All database manipulations and statistical analyses were performed using R v. 3.6.2 software [30]. The R package Nothopack v. 1.0.0 used for all simulations was developed by Gezan et al. [19] and can be accessed at https://github.com/sgezan/Nothopack on 10 September 2020.

Number of Trees Goodness-of-Fit
The simulation models studied were stand-level without diameter distribution (S1), stand-level incorporating a generation of a diameter distribution (S2), tree-level (T), yield compatibility (YC), and growth compatibility (GC). These were assessed based on different goodness-of-fit measures for their predictions of stand parameters such as the BA, QD, and VHA (Table 2 and Figure 2). As expected, all the simulations presented similar results for the NHA, given that all the simulations used the same mortality model. The relative residuals for the NHA show similar errors for the simulation models at 6 and 12 years of projection ( Figure 2A). The R 2 emp at six years presented a value of 0.97, and with 12 years, it suffered a minimal drop to 0.95. The RMSE% was 8.47 and 11.76% at 6 and 12 years of projection, respectively; but the simulation model, T, was slightly worse with 8.61 and 12.27%, respectively. The BIAS% was low for both projection years, observing a higher bias for simulation T but still with values close to zero. The Theil's inequality coefficient (U2) shows values close to zero, indicating a high performance in predictions of the simulation models where they have similar behavior (Table 2). Table 2. Goodness-of-fit statistics for each simulation (S1, S2, T, YC, and GC, see text for definitions) for number of trees (NHA, trees ha −1 ), basal area (BA, m 2 ha −1 ), stand volume (VHA, m 3 ha −1 ), saw-log volume (VOLP_SL, m 3 ha −1 ), pulpwood volume (VOLP_PW, m 3 ha −1 ), residual volume (VOLP_RES, m 3 ha −1 ), and carbon stock (C, tonnes C ha −1 ). Statistics are presented at 6 and 12 years of projection in second-growth RORACO forests in southern Chile.

Basal Area Goodness-of-Fit
In contrast to what is expected from the theory, where stand models should have better accuracies for BA than individual-tree models for long periods, the present T simulation resulted in better goodness-of-fit statistics. The R 2 emp values of 0.97 and 0.87 were found at 6 and 12 years of projection for the T simulation model, while other simulations gave lower values of 0.93 and 0.70 for the same periods. The RMSE% reached errors less than 6% at six years of projection and errors close to 11% at 12 years, showing high performance in the prediction of growth in the basal area. The BIAS% was low, with values close to 3% at six years of projection and 8% at 12 years. For the BA, the T simulation had a smaller bias and different behavior, presenting an under-estimation bias, while the other simulations presented over-prediction biases ( Table 2). The relative residuals exhibited similar errors for all the simulations, where the T simulation resulted in lower bias and an opposite behavior to the other simulations ( Figure 2B).
The QD is obtained using the NHA and BA; thus, the analysis is similar to these precedent stand parameters. However, it is possible to observe a small error at six years of projection in all the simulations, which increased at 12 years. Almost no bias was found at six years, changing to some over-estimation at 12 years, where, in all cases, the T simulation was less biased ( Figure 2D).

Volume Goodness-of-Fit
The observed volume growth reached in the PP was 8.81 m 3 ha −1 year −1 with a 95% CI of [5.78, 11.83] at six years and 11.21 m 3 ha −1 year −1 with a 95% CI of [7.25, 15.17] at 12 years. The VHA presented greater differences between the simulation models. As expected, S1 showed deficient performance, indicating that a general volume equation (Equation (8)) is not accurate enough to predict the behavior of this parameter in the validation sample available (Table 2, Figure 2C). In contrast to BA, volume is probably the best parameter for the tree-level simulation model where following each tree's development leads to better stand volume estimations. This is reflected in the T simulation that showed lower errors and bias for both time projections, 6 and 12 years, reaching R 2 emp of 0.94 and 0.83, RMSE% of 6.85 and 11.56%, and BIAS% of 0.11 and 5.05%, respectively.
Despite the added uncertainty in the generation of a diameter distribution and in the calculation of volume through these estimated diameter classes, S2 showed good predictions for both projection periods. At six years, S2 performed worse than the compatibility simulations, YC and GC; however, at 12 years of projection, S2 was more stable than the YC and GC simulations, and even presented lower bias than the T simulation model with values of 1.51 and 5.05% for 6 and 12 years, respectively. The compatibility simulations presented over-estimation bias against S2 and T that showed an under-estimation of VHA predictions (Table 2). Nevertheless, once observing the relative residuals, these simulations presented an error close to 0% at 12 years of projection, but the T simulations still presented the lowest prediction errors ( Figure 2C). These differences in VHA become more critical when the tree volume is subdivided into different stem products.

Products Volume Goodness-of-Fit
The following three objective products were defined to evaluate ecosystem services: saw-log (SL), pulpwood (PW), and residual lumber (RES). The simulation models' performances were contrastingly different depending on the type of product and the projection time. SL, a product associated with higher income, was better predicted for by the T simulation with an excellent R 2 emp of 0.98 and 0.99 at 6 and 12 years, respectively. The RMSE% was 11.77% at six years of projection, decreasing with 12 years to 7.24%. At the same time, the T simulation presented almost no bias at six years, marginally increasing to 2.61% at year 12. The U2 statistics were all close to zero, indicating a good prediction for saw-log volume ( Table 2).
On the contrary, the S2 simulation model had deficient performance for predicting SL; here, both the RMSE% and BIAS% were more deficient at six years (35.33 and 29.50%, respectively) than at 12 years of projection (19.63 and 11.96%, respectively). The compatibility simulation models presented better predictions than S2 at six years, but their behaviors were worst at 12 years of projection, indicating a high-volume over-estimation. The Weibull distribution diameter model used in building the model for the S2 simulations had difficulties estimating the proportion of large trees, and therefore, under-estimates the log-volume for larger trees (those subject to a minimum diameter of 24 cm) ( Figure 3A).
In the case of volume for the pulpwood product (PW), better simulation results were found for the compatibility models, where the GC was slightly poorer than the YC. However, the performance of the T simulation was similar to both compatibility models, having good statistics for R 2 emp (>0.92 at six and >0.79 at 12 years), RMSE% (<8.5% at six and <17.8% at 12 years), BIAS% (<1% at six and <8% at 12 years), and U2 (<0.08 at six and <0.15 at 12 years). As occurred previously with SL, neither of the S2 simulations presented reasonable PW estimations. For example, its R 2 emp value at six years was only 0.43; however, as with SL, the predictions improved at 12 years of projection with an R 2 emp of 0.68 ( Figure 3B, Table 2).
The volume for the residual lumber (RES) presented a more consistent behavior, where the T simulation showed a slightly overall better performance. For this product, all the simulations showed a decreasing BIAS% but an increased RMSE% at 12 years of projection ( Table 2). The relative residuals showed similar trends as the goodness-of-fit statistics but with lower bias for the compatibility models at six years ( Figure 3C).

Carbon Stocks Goodness-of-Fit
In general, the predictions of carbon stocks presented good performance, where the T simulations showed better results with an R 2 emp of 0.97 at six years and 0.89 at 12 years of projection. However, the compatibility simulations, YC and GC, showed similar statistics (0.96 and 0.86), for the same periods, respectively. Given that all the above-ground biomass equations (Equations (25)-(28)) used are strongly associated with the DBH and H, the S2 simulation had some difficulties to predict C, but still had reasonable accuracy with R 2 emp values of 0.74 and 0.80 at 6 and 12 years, respectively ( Table 2).
The statistic RSME% presented lower values at six years for the T, YC, and GC simulation models, but these doubled at 12 years of projection. In contrast, the errors of the S2 simulation stayed at~17% for both periods, but these were still higher than other simulations. S2 and T presented an under-estimation for C, while the YC and GC presented an overestimation in both periods of projection ( Figure 3D). However, the estimations were almost unbiased with values lower than 3 and 7% for both years of prediction, except for the prediction of C at six years using simulation S2 (13%). The carbon sequestration observed for the growth period of six years was 19.82 tonnes C ha −1 with a 95% CI of [14.98, 24.67]. In the case of 12 years of projection, the observed growth was 44.11 tonnes C ha −1 with a 95% CI of [20.19, 68.02]. The CI interval for carbon sequestration over 12 years is an indication of the high variability found in the remeasured plots; thus, simulations have to be able to deal with these diverse natural conditions. The T simulation presented the best performance at six years, increasing its under-estimation as sequestration values got larger. This is probably due to issues with under-estimation arising from the mortality model, or otherwise, due to the absence of an ingrowth component in the simulations (Figure 4A). The compatibility models showed an over-estimation at smaller carbon sequestration levels (<25 tonnes C ha −1 ) and showed better behavior with larger levels. The S2 simulation presented an under-estimation for the majority of the plots.
Both simulations at 6 and 12 years of projection had similar behavior, but with more errors in the estimation at 12 years, as expected ( Figure 4B).

Discussion
The model evaluations performed in this study, based on the data used for validation, suggest that stand level models, S1 and S2, present limitations for prediction, especially with volume (VHA). The use of the same mortality model for all the simulations allowed the isolation of its effects in the evaluation of the different simulation approaches. However, some problems associated with the prediction of other stand variables, such as the BA, VHA, QD, volume of products, and C stocks, could originate in the limited NHA predictions. Two important observations of this study are that these simulation models were fitted for Nothofagus species mainly, and they do not have an ingrowth module (hence, ingrowth is assumed to be zero). The model evaluations performed found problems in the NHA prediction once separated by species (N. alpina, N. obliqua, N. dombeyi, and companion species) ( Figure 5). The companion species had a poor performance compared to the Nothofagus species, accumulating errors and bias over the prediction periods. In addition, it was detected in the inventory data that in 2012, a few trees belonging to the companion species were observed; however, this recruitment was not considered in any of the simulated models (Figure 2A). Nonetheless, given these limitations, the NHA presented, in general, an excellent performance with small errors and an almost null bias given that RORACO forests, by definition, are mainly composed of Nothofagus species.
In most G&Y models, it is expected that the basal area (BA), over long periods, will be better predicted by stand-level models [4,6,31] than by tree-level models; however, in the present study, S1 and S2 did not show good accuracy. On the contrary, the treelevel simulation presented better BA estimations, not only with lower errors but also with less bias. This improvement is likely due to the incorporation of individual treecompetition predictors, such as BAL or RS, that deal better with the future response of RORACO forests [16]. This performance in the T simulation also suggests that the individual DBH growth for RORACO stands is highly correlated with its current size (DBH and H), competition status, and site characteristics. Similar conclusions have been reported in other mixed-species forest G&Y models [32,33]. Additionally, it is possible that the accuracy of stand-level BA growth models is affected by the higher tree density (NHA) of companion species, adding not only more uncertainty in their BA growth but, as mentioned before, affecting mortality patterns [11]. In contrast to the T simulations, the poor BA estimation of the stand-level simulation models affected the compatibility simulation models as well, given that these are fixed to the stand-level BA estimation. It is expected that the compatibility models will result in improvements in the BA estimation at the tree-level, as they will match the BA estimated using stand-level models. However, in the present study, the T simulation had a better performance than the compatibility models, YC and GC. This result contrasts with improvements in individual volume predictions when using the proportional growth compatibility method, as reported for Douglas-Fir stands [34]. Nonetheless, the YC and GC methods showed similar performance in predictions for stand parameters, with the GC slightly poorer than the YC. Other contrasting results have been reported for birch [6] and loblolly pine plantations [23], where the predicted goodness-of-fit statistics for the BA were better using the GC than the YC method.
The high resolution of the T simulation (tree-level), as expected, led to a better prediction of stand volume (VHA). On the contrary, the whole-stand simulation, S1, presented a low capacity to accurately predict the VHA. Similar behavior has been reported in other species [7,35]. However, the inclusion of a generated diameter distribution, as conducted in the S2 simulation, increased the ability to accurately predict the volume per hectare, a result that agrees with other studies [7].
The volume disaggregation required to increase the accuracy in the economic evaluation of ecosystem services, including saw-log timber and pulpwood, together with carbon sequestration estimations, showed different behavior depending on the simulation model, time of projection, and specific output analyzed. The T simulation was more successful in the prediction of the saw-log timber volume (SL), a positive result associated with its higher resolution and capacity to model thicker stem diameters properly. For saw-log volume, the stability in the DBH growth model translated into relatively constant errors with a low bias at 12 years. All the other simulation models presented poorer estimations with high increases in the RMSE%. However, for the small diameters associated with pulpwood volume (PW), the compatibility models, YC and GC, presented a slightly better performance than the T simulation. The product estimation errors for S2 simulation are higher than for other simulations, but also, for those products with small diameters, such as PW and RES. Here, the prediction tends to have decreased errors as the projection period increased. Qin and Cao (2006) found the same behavior, where the stand-level simulations showed a tendency to produce more accurate estimates for long-term growth.
The ability to predict carbon stocks was directly associated with the prediction of VHA; thus, the T simulation had better accuracy for both projection periods. The compatibility simulations also had excellent predictions with a bias lower than 7% at 12 years of projection, indicating a good performance from these models. S2, despite presenting poorer goodness-of-fit statistics, is still capable of predicting C with a reasonable R 2 emp value of 0.80 at 12 years of projection. However, S2 showed an underestimation at six years of projection, even simulating a bigger emission than the sequestration of C ( Figure 4A). This underestimation was induced by the low capacity in volume prediction, a high rate of mortality modeled, and the lack of ingrowth simulation.
The observed annual carbon sequestrations in this study were, on average, between 3.3 and 3.7 tonnes C ha −1 year −1 . Similar values were reported by Parada et al. [15] with 2.5 tonnes ha −1 year −1 for the carbon sequestration of above-ground tree dry matter biomass in a forest with the presence of the emergent N. dombeyi in Parque Nacional Puyehue, Chile. The values found in the present study indicate that the second-growth RORACO forest type, specifically the stands dominated by N. obliqua, are sequestering substantial amounts of carbon. However, due to the high variability found in the estimation of this variable, it is recommended, for studies focused on estimation of carbon sequestration for these forests, to use inventories with increased statistical power to reduce uncertainties.
As with any G&Y simulations, the compatibility adjustments, as well as the unadjusted predictions, suffer from an accumulation of errors and reduction in the goodness-of-fit statistics on their predictions as the simulation period increases. The compatibility methods are not only limited in their accuracy by the length of the simulation, but they are also strongly influenced by the accuracy of the stand-level simulations. Nevertheless, all the simulations evaluated had U2 values close to zero, indicating a better performance than the null hypothesis; that is, these models outperform the mean value estimation.
The sample of permanent plots used to validate the simulation models presented here incorporate the variability and different stand conditions of RORACO forests; however, these data are primarily dominated by N. obliqua, limiting inference to this forest type. Thus, there is a need for additional permanent plot data to validate further conditions and improve these models, particularly on the other two dominant species, N. alpina and N. dombeyi.
The G&Y system, packed in the R library Nothopack, to our knowledge represents the first effort to incorporate stand-level, individual-level, and compatibility models for the total distribution area of RORACO forests in Chile. The lack of time series or remeasured permanent plots limits this model, in aspects such as the incorporation of an ingrowth module or the inclusion of specific models for companion species. However, this is a dynamic process and further long-term data and evaluations should improve these aspects in the future, leading to an improved G&T system. In any case, this simulator is, without a doubt, a baseline for further improvements using and developing new models built with more robust information.
The implications of the results from Nothopack are useful for the sustainable management of second-growth RORACO forests in southern Chile, particularly given the importance and need of having some estimation of its value and a general idea of its management (thinning) that also includes some stimulus from the government. The principal users of Nothopack are small owners of these Nothofagus forests, forester consultants, and government agencies, as they require this type of tool to monitor the management of these resources, both for fiber products and ecosystems services.

Conclusions
The fitted individual-tree simulations show reasonably accurate predictions and a low bias for basal area, volume, specific products volume, and carbon stocks. The whole standlevel with the generation of a diameter distribution and both compatibility simulations displayed reasonable estimates, even with projections of 12 years. The present evaluations with independent remeasured plots confirm that the simulator Nothopack is a useful tool to support management plans for these forest ecosystems, where growth estimates might be limited to stands under 90 years old, or to merely update forest inventories over a short period (e.g., less than 10 years). Finally, this study provides the first evaluation and implementation of compatibility methods to link stand-and individual-level models in second-growth stands for RORACO forests in southern Chile.    (6)) 0.003595746 (α) 1.014 (k) ----QD 0max (Equation (7)) 11.6167(α 1 ) 11.3770(α 2 ) 11.7630(α 3 ) 11.6167(α 4 ) −1.4112 β -VHA (Equation (8)