Modeling growth and yield of loblolly pine stands under intensive management

The objective of this work was to develop and validate a prognosis system for volume yield and basal area of intensively managed loblolly pine (Pinus taeda) stands, using stand and diameter class models compatible in basal area estimates. The data used in the study were obtained from plantations located in northern Uruguay. For model validation without data loss, a three‐phase validation scheme was applied: first, the equations were fitted without the validation database; then, model validation was carried out; and, finally, the database was regrouped to recalibrate the parameter values. After the validation and final parameterization of the models, a simulation of the first commercial thinning was carried out. The developed prognosis system was precise and accurate in estimating basal area production per hectare or per diameter classes. There was compatibility in basal area estimates between diameter class and whole stand models, with a mean difference of ‐0.01 m2 ha‐1. The validation scheme applied is logic and consistent, since information on the accuracy and precision of the models is obtained without the loss of any information in the estimation of the models’ parameters.


Introduction
Natural forest formations in Uruguay are few and consist mostly of treeless pampas, known as the Campos grasslands of South America, concentrated primarily along river beds (Six et al., 2013).According to White & Pou (1980), due to the degradation of native forests, by 1980 about 130 thousand hectares of exotic forests had been planted in Uruguay, which was suffering from a lack of wood supply.The authors estimated that during this period Uruguay spent approximately US$16 million more importing wood products than it earned by exporting meat and derived products, the country's main export items.This scenario changed in the last three decades, when the country underwent a land-use transformation, from free-ranging cattle husbandry towards industrial fast-growing plantations, especially in the central and northern parts of the country (Vihervaara et al., 2012).This change was brought about by the implementation of forestry incentive laws (law 15939 in 1987 and law 16906 in 1998), which encouraged large-scale forestry plantations through the use of subsidies, tax relief, and targeted loans to investors (Mendell et al., 2007;Redo et al., 2012).Following these laws, there was a significant increase in the planted forest area of Uruguay, reaching 991 thousand hectares by 2013 (Uruguay, 2014).
Of the forest plantations in Uruguay, the majority is destined for pulp mills and domestic supply, and the main planted species are Eucalyptus grandis, E. globulus, and E. dunnii.The other major planted species in Uruguay are Pinus elliottii, P. pinaster, and P. taeda (Uruguay, 2014).Olmos & Siry (2009) estimated that 70% of Eucalyptus plantations in the country are managed for pulping and the rest for sawtimber production, whereas 100% of Pinus plantations are destined for sawtimber and plywood production.Many of the plantations managed for solid wood production are located in northern Uruguay, in the departments of Rivera, Paysandú, and Tacuarembó (Mendell et al., 2007;Six et al., 2014).In northern Uruguay, of the planted species, 52% are Pinus and 48% are Eucalyptus, differing from the species distribution in the rest of the country, which is of about 8% Pinus and 92% Eucalyptus (Uruguay, 2014).
The type of management applied in the plantations for solid wood production can be considered as intensive.Intensive forest management is defined as the production of high-quality timber in a short time frame, accomplished through the improvement of the stand's growth environment.This improvement is done mainly through soil preparation, fertilization, weeding, and thinning practices, among others.Intensive forest management has been proven to increase the quality and quantity of timber production in southern United States (Fox et al., 2007;Jokela et al., 2010), and is also a profitable investment option in other parts of the world.Rubilar et al. (2008) reported production gains of 20 to 50% and reductions in final rotation age of 2 to 5 years, for intensively managed Pinus and Eucalyptus plantations in Chile and Argentina.Bussoni & Cabris (2011) compared a 24-year rotation for mixed sawlog diameters to a 12-year rotation for small-log production of Pinus taeda stands in northeast Uruguay, and found that the land expectation value was about 50% higher for the longer rotation than for the shorter one.These authors also concluded that the investment in forest plantations without the adequate early silvicultural treatments, such as pruning and pre-commercial thinning, which characterize intensive management, result in a high risk of financial loss.
Given the long rotations of 22 years for Pinus and 16 years for Eucalyptus, associated with solid-wood production management regimes in Uruguay (Cubbage et al., 2014), growth models that predict future wood production are needed.The choice of a prognosis system to portray forest growth and yield depends on the level of desired detail, as well as on the management practices to be evaluated.Stand-level models are more suitable for forests destined for cellulose and energy wood production.However, when more detailed answers are required, as for forests for sawlog production, diameter class models are more indicated (Burkhart et al., 1981).Some of the most sought-after tools in forest management are mathematical models capable of estimating the impacts of different management alternatives, such as thinning regimes, applied to forest stands (Scolforo, 1993).
Examples of prognosis systems that depict Pinus taeda thinned-stand growth can be found in Eisfeld et al. (2005) and Costas et al. (2006Costas et al. ( , 2007) ) for Argentinian and for Brazilian plantations, respectively.For Uruguay, some studies on model sets for P. taeda have been published, which include volume and taper equations (Rachid et al., 2014), diametric frequency distribution functions (Hirigoyen & Rachid, 2014), and several mensurational relationships, such as height to diameter at breast height (DBH), crown diameter to DBH, crown volume to DBH, bark thickness to DBH, and bark thickness to tree height (Leites et al., 2013).
The objective of this work was to develop and validate a prognosis system for volume yield and basal area of intensively managed loblolly pine (Pinus taeda) stands, using stand and diameter class models compatible in basal area estimates.

Materials and Methods
The data used in the present study was obtained from approximately 60 thousand hectares of forest stands located in northern Uruguay, in the departments of Rivera, Tacuarembó, and Paysandú, situated between 60 and 230 m above sea level (31°22'S, 55°44'W).The data was obtained from plots of the Continuous Forest Inventory, which are circular in shape and vary in size from 300 and 500 m 2 .A total of 1,360 permanent plots were used, of which 163 were reversed for validation, comprising 7,164 measurement pairs.The characterization of the model parameterization and validation data was done (Table 1).The mean height of the 100 largest trees in DBH per hectare was used to represent the dominant height of the sample plots.
The majority of the plots (about 90%) used in the present study were less than 8 years old.This is an obstacle for the development of the prognosis system, since very few plots had been subjected to the first commercial thinning operation.The silvicultural operations planned for the stands are presented in Table 2.
Different models were used to estimate the following characteristics (Table 3): site index classification by the Chapman-Richards model in anamorphic form; stand volume by the Sullivan & Clutter (1972) model; stand tree survival; mean diameter variance by the adaptation of the Chapman-Richards model; and diameter class height model.
The last model needed to finish the prognosis system is a taper function, which allows estimates of the different product sizes that each diameter class can provide.Since a scaling database was not available in the present study, a taper model was not fitted.The models presented in Table 3 were chosen after exhaustive testing of several models.Additional information about the different tested models can be obtained in Ferraz Filho (2009).
A two parameter Weibull probability density function was used to estimate tree density in the different diameter classes.Due to its ability to simulate different distributions, such as the exponential and normal distributions, the Weibull distribution is widely used in many forestry studies, as those carried out by: Cao ( 2004), Merganic & Sterba (2006), Qin et al. (2007), Miguel et al. (2010), Petráš et al. (2010), and Soares et al. (2010).The probability density function of the Weibull distribution is obtained by the equation, in which: b>0 and c>0; b and c represent scale and form parameters, respectively; and x is the variable of interest, in this case DBH.
For the compatibilization between stand models (Table 3) and the diameter class model, the method of moments was applied to estimate the parameters of the Weibull function.The method of moments uses the first two non-central moments of the diameter distribution -mean diameter and quadratic diameter, respectively -to estimate the parameters of the Weibull function.
The shape parameter c is firstly recovered using the following equation: in which CV DBH is the coefficient of variation of DBH.Using parameter c, the scale parameter b can be recovered by the equation: DBH=b. .1+1 c ,

Γ ( ) in
which Γ is the gamma function.
As shown by these two equations, only the DBH and its variance are needed when the method of moments is used to recover the shape and scale parameters of the Weibull distribution.The key to archive the Table 1.Descriptive statistics of the dendrometric characteristics of the permanent plots of loblolly pine (Pinus taeda) evaluated (1) .compatibility between modeling resolutions is to obtain the DBH from the stand models.To achieve this, the future stand basal area should be estimated using the modified stand volume equation (Table 3): The future stand quadratic mean diameter can be estimated using the results from the basal area and survival projection models, using the equation below: Dg 2 = G/0.0000785398N.
Since the variables mean diameter, quadratic mean diameter, and diameter variance follow the same logic as a Pythagorean equation, the mean diameter at breast height can be estimated using the following equation: DBH = (Dg 2 -Sd 2 ) 0,5 .
With the stand's future DBH and diameter variance values, it is possible to estimate the shape and scale parameters of the Weibull distribution.
The validation of the prognosis system was made by separating an entire sub-region from the parameterization database.This guarantied spatial independence between the parameterization and validation databases (Vanclay & Skovsgaard, 1997).The equations were used to assess the precision, expressed as the mean absolute residual.
Validation was distinguished between young and mature input data.The young data was comprised of the first measurement available after the pre-commercial thinning, whereas the mature data was the penultimate measurement available.This was done to evaluate the behavior of the prognosis system when different projection lengths are used.Both the young and mature  input data were projected to the last measurement date, and the observed values were compared with the estimated ones.
After validation, the two databases were united and the models were reparametrized.This allows for any model inconsistencies to be detected in the validation process and for the resulting models to be fitted using all available data.The site index model was the only one fitted initially with the complete data set.
In order to evaluate the adherence of the estimated diameter distribution to the observed data, the Kolmogorov-Smirnov test was used in the 163 validation plots.
After the parameterization of the final models of the prognosis system, a thinning simulation was carried out (Figure 1).The simulation was done considering a thinning operation at 11 years of age, and the remaining forest data was projected to 17 years of age, i.e., the age of the second thinning operation (Table 2).
Since the used database did not have information on thinned plots, a thinning simulator was used (Arcebi Junior et al., 2002).This simulator estimates the proportion of trees to be removed in each diameter class, using the following equation, in which: P is the proportion of trees to be removed from the diameter class (DC) i; and β 1 and β 2 are the estimated regression parameters (β 1 = -0.49935and β 2 = 1.815105).A reference age of 10 years was used for site classification due to the young nature of the database.The sample plots were stratified into four different productivity classes with 5-m amplitude (Figure 2).The thinning simulator was fitted by Arcebi Junior et al. (2002) with data from the second thinning operation for P. taeda stands planted in the state of Paraná, Brazil.This simulator describes selective thinning from below in stands with less than 600 trees per hectare, and was chosen since it is the one that most resembles the thinning operations that are to be applied in Uruguay.To use this simulator, it is necessary to calculate the proportion of trees to be removed in the first diameter class of the stand, which should be done for each following diameter class until the desired number of trees to be removed is reached.

Results and Discussion
The dominant height values that distinguish the different productivities at the reference age are: 5 to 10 m (site index, SI=7.5); 10 to 15 m (SI=12.5);15 to 20 m (SI=17.5);and 20 to 25 m (SI=22.5).These values classify these plantations as being very productive.In southern United States, the dominant height range of planted loblolly pine varies between 4.7 and 10 m at 10 years of age (Diéguez-Aranda et al., 2006), which are values similar to the lowest productivity class found in the present study.The four site classes obtained in the present study are coherent with values reported for areas near the study location, in the states of Paraná and Santa Catarina, in southern Brazil.Scolforo & Machado (1988) recorded mean dominant height at the age of 10 years ranging from 7.7 to 22.2 m for P. taeda stands located in these same areas.
The initial parameters and goodness of fit statistics for the models that compose the prognosis system are shown in Table 4.All the parameters were statistically significant at a 99% confidence level, except for the β 2 parameter of the diameter class height model, which was significant at a 90% confidence level.
The accuracy and precision of the prognosis system for young and mature input plot values are presented in Table 5.The mean age for the young input data was 4.4 years, whereas for the mature plots it was 5.8 years.The mean projection age, the last measurement of the permanent plots, was 6.9 years.The accuracy of the models was satisfactory, with all models presenting mean residuals tending to zero and below 10% deviation, with the exception of the variance of DBH model.Precision was also considered satisfactory and presented values below 20% deviation, with the exception of the variance of DBH model.These values are in sync with the ones presented by Soares et al. (1995), who observed accuracy and precision estimates below 11% when modeling P. pinaster growth in Portugal.
The young input data had a tendency to result in less precise projection estimates, when compared to the mature input data.This behavior was expected since input data that is closer to the projection age is generally more correlated to the stand's future data (Coble et al., 2012).For instance, errors up to 6.17 m 3 ha -1 were found when projecting volume from young plots under 4.4 years of age.This error was reduced in -0.18 m 3 ha -1 when the input data was obtained from older plots with 5.8 years.
The less precise projection values of the young input data had an effect on the estimation of diameter distribution.Of the 163 validation plots, 11 did not present adherence to the observed values when using the mature input data, according to the Kolmogorov-Smirnov test, at 1% probability.For the young input data, the number of plots that did not present adherence rose to a total of 41.Hirigoyen & Rachid (2014), while evaluating Uruguayan P. taeda, also applied the Weibull probability density function using the method of moments and reported that 5.3% of the plots did not present adherence to the observed Table 4. Initial parameters and goodness of fit for the models of the prognosis system for loblolly pine (Pinus taeda) (1) .data.This value was very similar to the one of 6.7% found for the mature plots in the present study.The ability of the system to predict diameter distribution was largely affected by mortality estimates.For example, in the mature projections, 9 out of the 11 plots that did not present adherence to the observed data had a high mortality rate, above 360 trees per hectare in some plots.
The final model parameterization was carried out using the combined parameterization and validation database.The final parameters and goodness of prediction statistics of the models used in prognosis system are shown in Table 6.All the parameters were statistically significant at a 99% confidence level, except for the β 2 parameter of the stand volume model, which was significant at a 95% confidence level.
The parameterizations of the initial and final models were considered adequate, according to the statistical values of goodness of fit (Tables 4 and 6).The variance of DBH model presented the highest residual standard error (Syx=47.12%).Considering that the variance of DBH is a difficult characteristic to predict, the Chapman & Richards projection type model was considered satisfactory for this end.Eisfeld et al. (2005) also found that projection models outperform prediction models, when working with diameter variance models.The clutter volume model, or some modification of it, has been successfully used to model growth in (1) Y, young data; and M, mature data. (2)DBH, diameter at breast height.
Although slight, most of the statistical values of model quality presented a change for the worst from the initial to the final model parameterization.This occurred because more data was added to the parameterization database, increasing the distribution range of all the tree characteristics and, therefore, resulting in a greater variability.None of the fitted models presented any residual tendencies, evaluated by residual graphs (not presented), which might affect the quality of the predicted data.
As stated by Zhang et al. (2010), an important aspect of a prognosis system is the compatibility between estimates from different levels of resolution, in this case, between stand estimates and diameter class estimates.To check this, the difference in basal area estimates was tabulated using the validation database.The mean difference between basal area estimates using the stand model (Table 6) and the Weibull data was of -0.01 m 2 ha -1 , with values ranging from -0.05 to 0.16 m 2 ha -1 , confirming the compatibility between the different modeling approaches.
The thinning simulation was conducted using the parameters of the models related to the final parameterization (Table 6).The first step taken was to select four plots, with mean age of 5.7 years, that represented the mean characteristics of each site class.The thinning operation was conducted until a final number of 400 trees per hectare was obtained.The diameter distribution before and after the first commercial thinning operation at 11 years of age was determined for the two intermediate sites (Figure 3), which were chosen for being the most representative of the study area, comprising 94% total plots.
As can be expected in selection thinning from below, the greatest tree removal occurred in the smaller diameter classes, whereas the larger classes remained practically unaltered.There was an evolution of the basal area with and without thinning, observed through the four selected plots for the thinning simulation (Figure 4).This prognosis was made starting from 60 to 204 months of age, which corresponds to the age of the second commercial thinning.
The basal area values for the age of the second commercial thinning operation, for thinned and unthinned areas, are of similar magnitude (Figure 4), with the greatest difference occurring in the higher productivity sites.By 18 years of age, site 22.5 m presented a difference of about 10 m 2 ha -1 between thinned and unthinned areas, whereas site 7.5 m presented a difference of about 1.0 m 2 ha -1 .This can be explained by the fact that, in order to obtain a remaining stand of 400 trees per hectare, a greater amount of basal area must be removed from the higher productivity stands due to the larger size of the trees.
Figure 3. Loblolly pine (Pinus taeda) tree diameter distributions before and after the first commercial thinning simulation for sites 12.5 m (A) and 17.5 m (B).Pesq. agropec. bras., Brasília, v.50, n.8, p.707-717, ago. 2015 DOI: 10.1590/S0100-204X2015000800009 Even though the basal area values for the thinned areas are close to those for the unthinned ones, the former most likely will not pass the latter, confirming the results obtained in several thinning studies carried out by Hasenauer et al. (1997), Simard et al. (2004), Río et al. (2008), and Skovsgaard (2009).This statement can only be confirmed after the study plots undergo the first commercial thinning operation, allowing the modeling of post-thinning growth.
Although thinning reduces the basal area of the stand, individual tree growth is benefitted.The mean DBH value for site 22.5 m at 18 years of age in unthinned stands was of 40 cm, whereas in thinned stands it was of 43 cm.This result is in alignment with Zeide (2001), who found that thinning increases the growth of individual trees, despite reducing their number and the volume growth of the entire stand.

Conclusions
1.The prognosis system developed is precise and accurate in estimating loblolly pine (Pinus taeda) basal area production per hectare or per diameter classes.
2. There is compatibility in basal area estimates between diameter classes and whole stand models, with a mean difference of -0.01 m 2 ha -1 .
3. The validation scheme applied is logic and consistent, since information on the accuracy and precision of the models is obtained without the loss of any information in the estimation of the models' parameters.
used to assess the accuracy of the models, expressed as the mean residual; whereas y y n and y

Figure 1 .
Figure 1.Scheme illustrating the different stages of the thinning simulation.S, site index; N, number of trees per hectare; S 2 d, mean diameter variance; Hd, mean dominant height; G, basal area; and t, age.

Figure 2 .
Figure 2. Upper and lower limits of four site index classes at different ages for loblolly pine (Pinus taeda).The dashed line represents the reference age of 10 years.

Figure 4 .
Figure 4. Basal area projection for unthinned (A) and thinned (B) loblolly pine (Pinus taeda) stands in four different site indexes.

Table 2 .
Silvicultural operations planned for the evaluated stands of loblolly pine (Pinus taeda).

Table 3 .
Models used in the prognosis system for loblolly pine (Pinus taeda).

Table 5 .
Accuracy and precision of the prognosis system for loblolly pine (Pinus taeda).