Modeling Survival, Yield, Volume Partitioning and Their Response to Thinning for Longleaf Pine Plantations

Longleaf pine (Pinus palustris Mill.) is an important tree species of the southeast U.S. Currently there is no comprehensive stand-level growth and yield model for the species. The model system described here estimates site index (SI) if dominant height (H dom) and stand age are known (inversely, the model can project H dom at any given age if SI is known). The survival (N) equation was dependent on stand age and H dom , predicting greater mortality on stands with larger H dom. The function that predicts stand basal area (BA) for unthinned stands was dependent on N and H dom. For thinned stands BA was predicted with a competition index that was dependent on stand age. The function that best predicted stand stem volume (outside or inside bark) was dependent on BA and H dom. All functions performed well for a wide range of stand ages and productivity, with coefficients of determination ranging between 0.946 (BA) and 0.998 (N). We also developed equations to estimate merchantable volume yield consisting of different combinations of threshold diameter at breast height and top diameter for longleaf pine stands. The equations presented in this study performed similarly or slightly better than other reported models to estimate future N, H dom and BA. The system presented here provides important new tools for supporting future longleaf pine management and research.


Introduction
Longleaf pine (Pinus palustris Mill.) once dominated forests in the southeast U.S., occupying about 36 million ha prior to European settlement [1].By 1935, the area was reduced to about 8 million ha.Currently, there are only about 1.2 million ha of longleaf pine stands left [2], extending along the Gulf and Atlantic Coastal Plains from Virginia, south into central Florida, and north into the Piedmont and mountains of northern Alabama and Georgia [2].In recent years various organizations have begun promoting longleaf plantation establishment, directing most of their effort to private landowners with objectives that include production, but also aesthetics and wildlife habitat enhancement.
In order to improve stand management planning, researchers, managers and landowners need reliable information about stand dynamics and development.While a number of long-term experimental and permanent plot datasets exist for the species, these data, for the most part, have not been summarized in a comprehensive manner.As forest management decisions are based on information about current and future resource conditions, forest growth and yield modeling plays an important role by quantifying and summarizing relationships observed in field studies, and by providing stand projections under alternative management scenarios.Whole-stand-level growth and yield models predict future yields as a function of previous stand-level attributes such as age, stand density and site quality [3].
A number of models have been produced that predict elements of longleaf pine plantation stand dynamics [4][5][6][7][8].To our knowledge, however, no comprehensive stand-level growth and yield system has been produced for longleaf plantations that includes growth for thinned stands and merchantable volume estimations, and that can be applied to planted stands across a wide range of ages.
The objective of this study was to develop a stand-level growth and yield model system for thinned and unthinned longleaf pine plantations, using a long-term (>40 years) dataset measured and maintained by the U.S. Forest Service.A system of equations was developed to summarize the dynamics observed in the extensive, long-term dataset, and to provide a tool to predict and project stand growth and yield, including merchantable volume breakdown functions.

Data Description
The dataset used to estimate growth and yield parameters comes from 267 permanent plots measured and maintained by the U.S. Forest Service's Laboratory at Pineville, LA [9].The data were collected from regularly remeasured plots in a combination of seven studies exploring the effects of spacing and thinning on longleaf plantations distributed through the Western Gulf Coastal Plain from Santa Rosa County in Florida to Jasper County in Texas (Figure 1) and representing its current range in the Western Gulf Coastal Plain [9,10].Plantations were established on both old field and cutover sites.Soil texture for plots were primarily silt loams, very fine sandy loams, or fine sandy loams, characteristic of the U.S. Upper Coastal Plain.Most plots were burned regularly by prescribed burns or wild fires.Each plot was measured for ~40 years at ~five-year intervals, averaging eight measurements per plot.Plots were rectangular and ranged in size from 0.04 to 0.1 ha −1 [5].For each tree, stem diameter outside bark at 1.37 m height (dbh, cm) was measured to the nearest 0.25 cm and total tree height (H, m) was measured to the nearest 0.3 m on a subsample of 4 to 15 trees per plot.Using the model proposed by Quicke and Meldahl [11] that relates H to the inverse of dbh, individual plot by measurement time regressions were determined (P < 0.0001) to estimate H in all trees without H measurement (MAE% = 4.3%; RMSE% = 5.6%; Bias% = −0.4%;R 2 = 0.97).Mean dominant height (H dom , m) was determined for each plot at every measurement time as the mean of the top 25th percentile tree height.Site index (SI, m) was defined as H dom at a reference age of 50 years after planting.As several plots were not measured at exactly age 50 yrs, SI was assessed using H dom at index age plus/minus one year if necessary (i.e., 49 and 51 years).
From the complete dataset, 30 plots were randomly selected and removed to use for model evaluation and the rest (i.e., 237 plots) for model fitting.A total of 81 plots were thinned to constant basal area levels at five-year intervals; however, only pre-thinning measurements were considered on those plots.Summary statistics of individual trees and stand characteristics of both sub-datasets are shown in Table 1.For all trees with H measurement (28,083 observations) a form factor F = H/dbh (m cm −1 ) was calculated.In order to eliminate broken and malformed individuals, trees with F less than 0.54 m cm −1 and trees with F greater than 13.5 m•cm −1 were discarded for H dom determination (5.4% of H observations). Plots with less than four trees with H measured were also not considered for H dom fitting (7.9% of total plots).
The distribution of total observations (used for model fitting and model evaluation) by age, SI and surviving trees (N, ha −1 ) is shown in Table 2.Only 2% of plots had N > 2000 ha −1 , while 28% of data points had N < 500 ha −1 .About 29% of data had N between 1000 and 1500 ha −1 .Most stands (52%) had SI between 25 and 28 m, and only 4% of data had SI < 22 m.About 22% of stands had SI > 28 m.In term of age distribution, about 20% of data had age < 20 yrs. and 9% had age > 60 yrs.

Survival and Yield Models for Unthinned Stands
Data from all unthinned plots and from thinned plots prior to thinning were used to estimate survival, H dom , stand basal area and stand volume (inside and outside bark) parameters.
Following the guide curve method to produce an anamorphic model [12], a dominant height function was fitted based on the Chapman-Richards function using the following expression: where Age is the stand age (yrs.), a 0 , a 1 and a 2 are curve fit parameters and ε 1 is the error term, with ).For the selected site index age of 50 yrs., the model can be re-written as: This equation can be inverted to determine SI if stand age and H dom are known.This anamorphic model has the assumption that the shape of the height-age curve is independent of SI, and differences between any two curves are proportional to the ratio of their SI's [13].
A negative-exponential survival model that includes H dom was used to estimate survival using a modified version of the model proposed by Dieguez-Aranda et al. [14,15] and Zhao et al. [16]: where N j is the number of trees ha −1 at age j (yr.), N i is the number of trees ha −1 at age i (yr.) (i < j), H is the dominant height (m) at age i (yr.), b 1 to b 4 are curve fit parameters and ε 2 is the error term, 2 ).Several models proposed by Dieguez-Aranda et al. [14], Zhao et al. [16] and Burkhart and Tome [17] were also tested, but the model that we selected showed the best fit.After model testing, similar to Zhao et al. [16], the parameters b 1 and b 3 were set equals to 0 and 1, respectively, due to no improvement in predictive ability and convergence difficulties.The final model to estimate H dom was: The following generic equation, proposed by Borders [18], was initially used to predict basal area: where BA is the stand basal area (m 2 •ha −1 ), N is the survival (trees ha −1 ), c 1 to c 7 are curve fit parameters and ε 3 is the error term, with ε 3 ~ N(0, σ 3 2 ).After step-wise procedure and variance inflation factor (VIF) analysis, parameters non-significant and/or with high multicollinearity were discarded, resulting in the following final model to estimate BA: The equations reported by Gonzalez-Benecke et al. [19], which depend on the individual dbh and stand parameters N, H dom and SI, were used to estimate individual tree volume outside and inside bark for each living tree in the dataset.After aggregating all individual tree volumes within each plot, stand volume outside and inside bark was determined for each plot.This information at the plot level was used to fit a model for stand volume prediction, which was initially based in the following generic model proposed by Borders [18] and Pienaar [20]: where VOL is the stand stem volume outside or inside bark (m 3 •ha −1 ), d 1 to d 9 are curve fit parameters and ε 4 is the error term, with ε 4 ~ N(0, σ 4 2 ).Similar to BA, after the step-wise procedure and VIF analysis, parameters non-significant and/or with high multicollinearity were discarded, resulting in the following final model to estimate VOL:

.2. Basal Area Growth Model for Thinned Stands
As the effects of thinning on survival and H dom are small for southern pines (Westfall and Burkhart, 2001;Sharma et al. 2006) [21,22], we only modeled the response in BA growth after thinning.Several models, presented in Burkhart and Tome [17], were also tested in order to simulate BA growth after thinning, but the methodology reported by Pienaar [20,23] was selected.Following Pienaar [23], BA projection for thinned stands (BA t , m 2 •ha −1 ) was determined by using a competition index (CI) and the basal area of an unthinned counterpart stand (BA u , m 2 •ha −1 ), assuming that BA t can be expressed as a proportion of the basal area of an unthinned stand of the same age, H dom and number of surviving trees (i.e., an unthinned counterpart) that changes over time.The CI is the rate of competition decline, a measure of the relative degree of competition affecting tree size in the thinned compared to the unthinned stands, and it was determined as follows: From Equation 7, when BA t equals BA u with the same number of trees, the CI is zero.Similarly, when BA t is less than BA u (that is the general case in operational thinnings), the CI is larger than zero, but approaching zero as stand ages, as BA t will converge to BA u [20].As the permanent plots of thinned stands do not have an unthinned counterpart, projections of BA u over time were estimated using Equations 2, 4 and 5.The BA growth response after thinning was determined indirectly by projecting the time trend of CI, assuming an asymptotic trajectory towards a value of zero.Thus, reflecting that the thinned stand, which has the same age, H dom and number of trees as the unthinned counterpart, will approach, over time, the unthinned stand in terms of total BA [23].The projected CI after thinning was estimated using a modified version of the model proposed by Pienaar [23], including the effect of stand age on the rate of competition decline: where CI i and CI j are the competition index at age i and j (yr.) (i < j), respectively, and f 1 is the curve fit parameters and ε 5 is the error term, with ε 5 ~ N(0, σ 5 2 ).The exponential of the coefficient term , A , corresponds to the annual decline rate of the CI as the stand ages after thinning.After combining Equations 5 and 6, BA of a thinned stand is estimated using the projected CI as: where BA and BA are the projected BA (m 2 •ha −1 ) in the thinned and unthinned counterpart stands at age j (yr.), respectively.

Merchantable Volume
For each tree in the fitting dataset, merchantable stem volume (both outside and inside bark), from the stump to any top diameter, was estimated using the equations reported by Gonzalez-Benecke et al. [19] for a range of combinations of threshold dbh values (d, from 5.08 to 40.64 cm) and top diameter limit values (t, from 5.08 to 45.72 cm) that incremented at steps of 5.08 cm.Finally, for each plot, merchantable stem volume per hectare for each combination of d and t was calculated based on all living trees within each plot.
The merchantable volume yield breakdown at the stand level function was determined following Amateis et al. [24], where total volume yield outside or inside bark (i.e., VOL OB and VOL IB , m 3 •ha −1 ) was proportionally assigned to product classes defined by two variables: top stem diameter outside bark merchantability limit (t, cm) and a dbh threshold limit (d, cm): where VOL m is merchantable volume per hectare (m 3 •ha −1 ) for trees with dbh larger than d cm, to a top diameter of t cm outside bark, VOL is the total volume per hectare (m 3 •ha −1 ), Dq is the quadratic mean diameter (cm), N the survival (trees ha −1 ), and m 1 to m 5 are curve fit parameters.

Model Evaluation
The predictive ability of the fitted models for N (Equation 4), H dom (Equation 2), BA of unthinned stands (Equation 6), VOL of unthinned stands (Equation 8) and BA of thinned stands (Equation 11) was assessed with the 30 plot evaluation dataset.In the case of VOL m (Equation 12), two combinations of d and t were selected to evaluate this model.The selected values of d and t correspond to the threshold values used to estimate breakdown volume yield for other southern pine species [25][26][27], corresponding to chip-and-saw (d = 21.6 cm; t = 10.2 cm) and sawtimber (d = 29.2cm; t = 20.3cm) products.
Four measures of accuracy were used to evaluate the goodness-of-fit between observed and predicted (simulated) values for each variable originated from the dataset obtained in the model evaluation: (i) mean absolute error (MAE); (ii) root mean square error (RMSE); (iii) mean bias error (Bias); and (iv) coefficient of determination (R 2 ) [28][29][30][31].For BA and VOL, the statistics MAE, RMSE and Bias were back-transformed from logarithmic values.
The system of equations was also compared against other models reported in the literature for longleaf pine plantations using the same model evaluation dataset indicated above.The models compared were: (i) survival equations reported by Lohrey and Bailey [5], Brooks and Jack [7] and Lauer and Kush [32], (ii) dominant height equations reported by Farrar [4] and Brooks and Jack [7], and (iii) VOL OB equations reported by Lohrey and Bailey [5] and Brooks and Jack [7].The breakdown volume outside bark yield function was also compared against the functions reported for Pinus taeda [25] and Pinus elliottii [26] across different Dq and stand densities, and for three product classes (sawtimber: d = 29.2cm and t = 20.3cm; chip-and-saw: d = 21.6 cm and t = 11.2 cm; pulpwood: d = 11.4 cm and t = 5.1 cm), assuming stands with VOL OB of 100 m 3 •ha −1 .This value of VOL OB was selected to facilitate percent comparisons.The results of the breakdown volume yield function are independent of the value assumed.
An overall evaluation of the model was carried out for unthinned plots of the validation dataset.On each plot, for known initial stand age, N and SI, stand BA and VOL OB were estimated on the same ages where they were originally measured by using the final equations fitted to estimate N, H dom , BA and VOL OB .The same four measures of accuracy described previously were used to assess the agreement between observed and predicted values.
All of the summary, model fitting and model evaluation statistics were obtained using SAS 9.3 (SAS Inc., Cary, NC, USA) [33].When multiple linear regressions were carried out, the variance inflation factor (VIF) was monitored to detect multicollinearity between predicting variables, discarding all variables included in the model with VIF larger than 5, as suggested by Neter et al. [34].In the case of BA and VOL fitting, where multiple linear regressions were carried out, step-wise procedure was used with a threshold significance value of 0.15 as variable selection criteria to enter and to stay [34].For these responses, a logarithm transformation was preferred as it allows controlling for heterogeneity of variances, approximate to normality and uses the linear model framework to select among the large set of, potentially collinear, predicting variables.

Model Application Example
The system of equations developed was used to predict stand growth of unthinned and thinned (3 thinnings, removing 33% of living trees at ages 30, 40 and 50 yrs.)longleaf pine stands growing at sites with two different SI's: 20 and 30 m.The initial planting density used was 1400 trees ha −1 , the survival after the first year was assumed to be 95%, and the simulation length was 70 yrs.Here it was assumed that the percentage of removed trees was the same as the percentage removal of BA during thinning.

Model Fitting
The parameter estimates for the growth and yield predictive and projective equations (Equations 2, 4, 6, 8, 10 and 12) for longleaf pine plantations growing in Western Gulf Coastal Plain U.S. are reported in Table 3.All parameter estimates were significant at P < 0.05.Non-linear versions of the models presented in Equation 6(BA) and 8 (VOL) were also evaluated, but these resulted in no improvement in model performance (data not shown), therefore, natural logarithm-transformed response variables were used.Parameter estimates for the intercept in Equations 6 (c 1 ) and 8 (d 1 ) include the correction proposed by Snowdon [35].The correction factor proposed by Baskerville [36] was also evaluated, but it presented lower bias reduction (data not shown).
Parameter estimates for the model that projects dominant height (Equation 2) are shown in Table 3.For SI of 20 m, the model projects dominant height of 5.4 and 22.6 m at age 10 and 70 yrs., respectively.If SI increased to 29 m, the model projects dominant height of 7.9 and 32.8 m at age 10 and 70 yrs., respectively.For all 194 plots where SI was measured, the mean observed and predicted SI was 25.84 and 25.88 m, respectively.
The survival model (Equation 4) was dependent on stand age and H dom .The performance of the N model for the range of SI present on the dataset used for model fitting (i.e., between 20 and 29 m, see Tables 1 and 2) and using a planting density of 1500 trees ha −1 showed little mortality and only small differences in survival at age 10 yrs.(between 1450 and 1429 trees ha −1 , for SI 20 and 29 m, respectively).At age 70 yrs., however, the model estimated large differences in survival across SI's (between 493 and 300 trees ha −1 , for SI = 20 and 29 m, respectively).
After applying the step-wise procedure and checking variance inflation factors (VIF), the final selected model that predicts BA (Equation 6) was only dependent on N and H dom (Table 3).Although the variables 1/Age, ln(N)/Age and ln(H dom )/Age were significant after the step-wise variable selection procedure (P < 0.001, data not shown), their VIF's were high with values of 296, 173 and 35, respectively (data not shown).Therefore, these variables were discarded from the model and the goodness-of-fit of the final model was lower than the full model, having a CV of 3.6% and a R 2 of 0.944.Partial R 2 of H dom and N were 0.579 and 0.366, respectively (data not shown).
The final model that predicts stand volume (Equation 8), after the step-wise variable selection procedure and VIF checking, was dependent on N, BA, ln(BA)/Age and SI (Table 2).The variable H dom was discarded from the final model, even though they were selected after step-wise variable selection procedure (P < 0.001, data not shown), due to its high multicollinearity (VIF = 42.3,data not shown).The final models that predict stand V OB and V IB had a CV of about 1% and R 2 greater than 0.99.Stand BA explained most of the variability in V OB and V IB , with partial R 2 of about 0.912 and 0.867, respectively.Stand density presented partial R 2 of about 0.079 and 0.121, for V OB and V IB , respectively.Even though SI and ln(BA)/Age were significant, both explained less than 0.1% of changes in stand volume (data not shown).The model that projects the time trend of CI after thinning (Equation 10) was dependent on stand age (Table 3).The exponential of the coefficient for a 36 year-old stand (mean values of stand age reported in Table 1), represents an average annual decline rate of the CI as the stand ages after thinning of 4.3%.
Parameter estimates for merchantable volume yield breakdown function (Equation 12) for both, outside (VOL m-OB , m 3 •ha −1 ) and inside (VOL m-IB , m 3 •ha −1 ) bark total volume yield, are shown in Table 3.The models had a CV of about 11% and an approximate R 2 of 0.99 for both outside and inside bark volume yield breakdown estimates.
There was a good agreement between predicted and observed values of N (Figure 2a), H dom (Figure 2c) and BA (Figure 2e).The slope and the intercept of the relationship between predicted and observed values were not statistically different from one (P = 0.32) and zero (P = 0.14), respectively.If residuals are expressed as a percentage of the observed value, maximum absolute residuals observed represent about 17% and 16% of observed N and H dom , respectively.Residuals for predicting the BA model were larger, with maximum residuals of about 30% of observed BA, but centered around zero.There was no noticeable trend in residuals with observed values (Figure 2b for N; Figure 2d for H dom and Figure 2e for BA) or stand age (data not shown).
A growth model to project, or update, BA (BA j ) when some stand measurements are available, including current BA (BA i ), along with current (i) and future (j) Age, N and H dom , was derived from the fitted BA model (Equation 5), and is expressed as: As expected, this projecting growth model improved the estimations of BA, reducing the residuals as compared with BA predicting model (Figure 2g,h).
There was good agreement between predicted and observed values for VOL outside and inside bark (Figure 3a,c).In both cases, the intercept of those relationships was not statistically different from zero (P > 0.16).The slopes of the relationship between predicted and observed values (1.006 and 1.011, for VOL OB and VOL IB , respectively) were statistically different to one (P < 0.02).There was a good agreement between predicted and observed values of BA t (Figure 3e).The intercept and slope were not different from zero (P = 0.18) and one (P = 0.12), respectively.If residuals are expressed as a percentage of observed value, maximum residuals observed in Figure 3 represent about 15% and 17% of observed VOL OB and VOL IB , respectively.There was no noticeable trend in residuals with observed values (Figure 3b for VOL OB ; Figure 3d for VOL IB and Figure 3e for BA t ) or stand age (data not shown).
For the two combinations of t and d tested, there was good agreement between predicted and observed VOL m (Figure 3g,i).For the two examples shown in Figure 3, the slope and intercept of those relationships were not different from one (P > 0.41) and zero (P > 0.36), respectively.For all variables listed above, there was no noticeable trend of residuals with observed values.Only for sawtimber VOL m larger than about 500 m 3 •ha −1 (Figure 3i,j), there was a small tendency to increase residuals as VOL m increased, but the magnitude of that overestimation was less than 5% of observed values.All model performance tests showed that N, H dom , BA, BA t , VOL and VOL m estimations agreed well with measured values (Table 4).For all estimations, MAE and RMSE ranged between 2.4% to 10.3%, and 3.4% to 13.4% of the observed values, respectively.In all cases, BA estimations presented the larger differences between the observed and predicted values.The Bias ranged between 1.9% under-estimations for projected BA and 1.4% over-estimations for BA t , with no clear tendency to over-or under-estimate.Estimated and observed values were highly correlated, with R 2 values greater than 0.91.The performance of the BA model that included the variables with high collinearity was also tested, showing lower MAE and RMSE than the final model (8.8, 12.6%, respectively) and larger absolute Bias (−3.3%) (data not shown).When tested on the dataset used for model evaluation, predicted values of the models proposed in this study for N, H dom and VOL OB are within the range of variation of the estimations using other published growth and yield models.The effects of stand age on survival, H dom , and VOL OB were predicted using several models for longleaf pine (Figure 4).Across three stand age classes (<25, 25-49 and 50-75 yrs.), the models predict stand growth consistently, with no clear trend to over-or under-estimate.For example, N and H dom estimations of all models performed adequately with Bias less than 10% (Figure 4a) and RMSE less than 15% with no apparent trend to change across stand age classes (Figure 4b).The estimates of VOL OB were similar for all models for age less than 25 yrs., but for older stands the model reported by [7] Brooks and Jack (2006) over-estimated VOL OB by around 70 m 3 •ha −1 , while the model reported by [5] Lohrey and Bailey (1977) over-estimated VOL OB by about 50 and 10 m 3 •ha −1 for stand age between 25-49 and 50-75 yrs., respectively.

(V3).
Examples of merchantable yield breakdown function of VOL OB estimations for P. taeda [25], P. elliottii [26] and P. palustris (this study) are presented in Figure 5, showing, notably, similar results across species.For example, for sawtimber, defined as stem volume of trees with dbh larger than 29.2 cm outside bark (threshold dbh limit) to a top diameter of 20.3 cm outside bark (merchantability limit), when Dq was smaller than 20 cm there was no sawtimber volume production, but when Dq was 30 cm, sawtimber yield was about 55, 61 and 52% (N = 100 trees ha −1 ) or 57, 67 and 73% (N = 1000 trees ha −1 ) for P. taeda, P. elliottii and P. palustris, respectively (Figure 5a,b).In the case of chip-and-saw yield, when Dq was 10 cm, all models predicted no volume production for that product, which has a threshold dbh limit of 21.6 cm.Independent of N, when Dq was larger than 50 cm, sawtimber yield was larger than 95% of VOL OB and the production of chip-and-saw (Figure 5c,d) and pulpwood (Figure 5e,f) declined when the stands reached Dq larger than the merchantability limit for sawtimber.Overall, the merchantable yield breakdown functions presented in this study showed the expected behavior of product partitioning as Dq and N changed.The overall test of the model indicated that, if only initial (i.e., current) stand age, N and SI (reported at first measurement) are known, estimations of H dom , BA and VOL OB were not affected, in relative terms, by simulation length (Figure 6d,f,h).On the other hand, projections of N were sensitive to the length of the simulation (Figure 6b), with errors getting larger as simulation length increased.In all cases residuals were centered on zero.For all stand parameters simulated, there was no noticeable trend of residuals with observed values, and the slope and intercept of the relationships between observed and predicted values were not different from one (P > 0.24) and zero (P > 0.42), respectively.If initial stand age, N and SI are known, the overall test of the model system indicated that projections of N, H dom and predictions of BA and VOL OB for less than ~40 yrs.simulation length presented a bias that ranged between −7% and 10% (Table 5).Across simulation lengths, the overall bias of the model system for N, H dom , BA and VOL OB were over-estimations of about 6 trees ha −1 and under-estimations of about 0.2 m, 0.9 m 2 •ha −1 and 13.4 m 3 •ha −1 , respectively.The overall MAE and RMSE of the model system were about 12 and 18% for N, 3 and 3% for H dom , 12 and 16% for BA and 13 and 18% for VOL OB , respectively.The R 2 decreased as simulation length increased.The overall R 2 across simulation lengths were about 0.93, 0.97, 0.84 and 0.84, for N, H dom , BA and VOL OB , respectively.A trend of increasing error with simulation length was observed for Bias, MAE and RMSE (Table 5).Nevertheless, it is important to note that the number of observations decreases as the simulation length gets larger (Table 5), and therefore, the evaluation statistics on simulation lengths ~40 yrs.can be affected by the unbalanced sampling size and specific characteristics of the sampled plots for evaluations.An example of model behavior for a hypothetical longleaf stands planted with 1400 trees ha −1 is shown in Figure 7.The unthinned stands growing on a site with SI = 20 m reached at age 70 yrs.a survival of about 46% of initial planting density, H dom of 22.6 m, BA of 28.5 m 2 •ha −1 , SDI of 571 tress ha −1 and VOL OB of 338 m 3 •ha −1 .When SI was 30 m instead, at the same age the survival was 31% and H dom was 33.8 m.The unthinned stand reached a maximum BA of about 48.1 m 2 •ha −1 at age 58 yrs., and SDI peaked at about 870 trees ha −1 at age 49 yrs.and VOL OB was still increasing, reaching about 610 m 3 •ha −1 at age 70 yrs.When a scenario of three thinnings was applied to both stands, at age 70 yrs.the number of surviving trees was about 181 and 122 trees ha −1 and Dq was increased from 24.6 and 38.3 cm (unthinned), to 30.2 and 47.2 cm, for SI of 20 and 30 m, respectively.The harvested volume from thinnings was about 171 and 331 m 3 •ha −1 , and the final yield was 162 and 293 m 3 •ha −1 for SI of 20 and 30 m, respectively.

Discussion
Bringing existing longleaf pine stands under management and restoring longleaf pine stands from degraded or otherwise converted forest stands is a priority for a number of land management entities in the southeastern U.S. [37,38].Managers undertaking these tasks must have information about the response of growth and stand structure under alternative silvicultural scenarios.Growth and yield systems which incorporate long-term data from stands on a variety of sites and under a range of management regimes provide one of the best tools for exploring the possible outcomes of proposed management regimes.While a number of models predicting elements of longleaf pine plantation stand dynamics have been produced [4][5][6][7][8], this study represents, to our knowledge, the first comprehensive stand-level growth and yield model for longleaf pine plantations, including stand growth and merchantable volume estimations, that can be applied to plantations across a wide range of ages (from 7 to 73 years) and site quality (SI ranging from 20 to 29 m).All choices of model structure involve compromise.Whole-stand level models, as the one presented here, provide reliable prediction of stand variables, such as BA and N; on the other hand, they do not provide the level of detail that individual-tree level models produce, which could allow for more flexibility in modeling silvicultural practices.However, individual-tree models typically are unreliable in prediction of cumulated stand information, and often have issues with propagation of errors.In this study, we opted to fit a stand level model to be used as a baseline, and in future work, we will consider incorporating individual-tree level information.
Site index is the most widely used measure of forest productivity, particularly in plantations.Base age selection for SI can have significant implications for the accuracy of estimations, as bias increases as the stand age is further from the base age [13].For this study we decided to set SI at age 50 yrs., a widely used reference age in Southeast U.S. [4,5,32].The H dom model reported in this study, which behaved well for a wide range of stand age, performed similar or slightly better than the models reported by Farrar [4] and Brooks and Jack [7].However, larger bias using the Brooks and Jack [7] model could be a result of applying their equation out of the age range and geographic zone of inference, as it was fitted from stands in southwest GA, with ages between 2 and 19 years.Dominant height is a major component of yield prediction systems for southern pine plantations.The anamorphic model obtained by this study seems suitable for H dom estimations in the Western Gulf Coastal Plain, U.S.
In relation to the survival equations, the best model fitted was dependent on stand age and dominant height (a measure of site quality).Other models also incorporate the effect of site quality on survival, such as the models reported by Lauer and Kush [32] and Farrar and Matney [39] for naturally-regenerated longleaf pine stands, or the models reported for related southern pine species P. elliottii [40] and P. taeda [25,41].A model that was only dependent on stand age, similar to that reported by Brooks and Jack [7] or Pienaar et al. [26], was also fitted, but even though the resulting equation did perform well across all stands included in this study, the final model selected had a slightly better fit (data not shown), and at the same time allows the inclusion of the effect of site quality (reflected in H dom ) on resource competition: the larger the SI (and hence H dom ), the larger the mortality rate after canopy closure.This process of accelerated self-thinning has been well documented for southern pines [16,[42][43][44][45].For example, 25 yrs.old P. elliottii and P. taeda stands in fertilized plots (with higher SI) had an accumulative mortality of about 59 and 43%, respectively, while non-fertilized plots showed lower mortality of about 43 and 22%, respectively [45].Murphy and Farrar [46] reported that models that include H dom performed better than models that rely only on stand age to project survival, especially on prepared sites.
Other models that included SI [14][15][16][17] were also tested, but H dom was selected due to better predictive ability.An attempt was made to model survival as a two-step modeling approach [16], including an equation to predict the probability of survival of all trees in the stand over a measurement interval, but no improvement was observed.The model presented in this study performed similar to, or slightly better, than other reported models to estimate future survival, but performance can be influenced by the fact that, even though the validation plots are independent to the plots used for model fitting, they were located in the same geographic region and we are not using the model out of the geographic inference zone, as we did with the model of Brooks and Jack [7].Nevertheless, the model of Lohrey and Bailey [5] was fitted with a subset of the same dataset used for this study, but perhaps the lower range of ages and N in their dataset influenced the results presented here.On the other hand, the model of Lauer and Kush [32] performed well across different combinations of stand density and productivity in this study.Residual analysis for each model did not indicate any unusual trends, even though the stands cover a wide range of productivity, planting density and age.The models obtained by this study appear suitable for survival estimations within the Western Gulf Coastal Plain U.S.
The model that presdict BA for unthinned stands was only dependent on N and H dom .Other models reported for southern pine species include other simple or composite variables such as SI, Age, N/Age and H dom /Age as well.In this study, all of these variables were significant and could be included in the model, but the high multicollinearity between those variables indicated the need to drop them from the final model.Therefore, it was decided to discard those variables, obtaining a simpler model to predict BA for unthinned longleaf pine stands with similar goodness-of-fit which avoided over-fitting problems.Widely used models for other southern pines [25,26] and for longleaf pine [32] include other stand variables to predict BA, but the existing publications contain no information regarding multicollinearity and/or their relative importance.The model that projects BA reported in this study had the same structure of the model reported by Brooks and Jack [7].
In the case of thinned stands, the use of the approach proposed by Pienaar [20,23], that uses CI and the basal area of the unthinned counterpart, fitted the available data well allowing the estimation of BA growth after thinning.Modified versions of the model used were also evaluated, which included H dom or SI as dependent variables to modify the rate of decline of CI as stands ages after thinning, but none of those variables were significant into the model (data not shown) and the final model that estimates the decline of the CI as the stand ages after thinning was only dependent on stand age.For 25-and 50-year-old stands, the mean value of the annual rate of decline of the CI as the stand ages after thinning were 6.2% and 3.1%, respectively, reflecting the impact of stand age on the rate of decline of CI as stands ages after thinning.The value of the rate of decline of CI at age 50 yrs.is lower than the values reported for other southern pine species.For example, Pienaar [20] and Harrison and Borders [25] reported values of 9.3% and 7.6% for P. elliottii and P. taeda, respectively.Those authors reported CI's for stands thinned at younger ages of about 13 years and measured until age 30 in the most extended study [20].The dataset used in this study included plots thinned at ages ranging from 17 to 63 years, and measured, in average, for about 20 yrs.The slope of the relationship between observed and predicted BA after thinning was not different from one, supporting the robustness of the model that projects BA growth for thinned stands.
The models that predict VOL (outside and inside bark) did not depend on stand age or N, and were similar or slightly better than other models reported for longleaf plantations [5,7] or naturally-regenerated stands [32].In this study, as done previously with the BA model fits, variables that showed high levels of multicollinearity were dropped to obtain a parsimonious final model that was only dependent on BA and H dom and provides good prediction ability.The residuals of this model showed a tendency toward under-estimation for observed VOL OB greater than 600 m 3 •ha −1 , a condition which is rare in longleaf plantations, which typically include multiple thinnings [5,[47][48][49].
Nevertheless, the maximum residual found was less than 8% of the observed value, with residuals centered near zero at different age classes.For the dataset used for validation, the models reported by Brooks and Jack [5] and Lohrey and Bailey [7] showed very good prediction ability only for stands up to 25 yrs.old.This may be explained by the fact that the former was developed for stands younger than 20 yrs. in age and the latter includes a reduced range of ages and site quality as compared with the dataset used in the present study.This indicates that the new data used here, expanded the ability to accurately predict stand volume if compared with the model of Lohrey and Bailey [5].
One of the most important contributions of the system of equations presented in this study, in contrast to other longleaf models published, is the inclusion of the merchantable volume yield breakdown model.Similar to P. elliottii and P. taeda, it is now possible to estimate the merchantable volume for different combinations of threshold dbh and top diameter for longleaf pine stands by using equation 8. Although, the model presented in this study made similar predictions to some models currently used for P. elliottii and P. taeda, differences in diameter distribution between species could explain the disagreements observed in merchantable volume yield breakdown, especially in chip-and-saw and pulpwood products.
The overall evaluation, where N, H dom , BA and VOL OB , were calculated for all unthinned plots using the system of equations shown in Table 3 starting with a known initial age, N and SI (reported at first measurement), demonstrates the robustness of the model.As was expected, errors tended to get larger as the length of simulations increased, but overall, the residual centered on zero.Even though the number of observations was not balanced across simulation length classes, a fact that may explain the bias found for simulation lengths greater than 40 years, the results indicate that we can expect accurate estimations for simulation length of up to 40 years.Therefore, it is recommended that users update their stand inventories at least each 30-40 years, in order to improve the predictions from the use of the models presented here.
Despite the fact that the model system performed very well for the dataset used for validation, the functioning of the model outside the geographical range of the fitting data is uncertain.We strongly recommend using this system of equations only within the range of data used to fit (see Tables 1 and 2).In addition, the model does not include the effects of genetics, site preparation, weed control and fertilization.Nevertheless, the system presented here provides an important new tool for supporting present and future longleaf pine management decisions.Future research expanding the area of inference and including the effects of genetically improved material and intensive silvicultural management is needed to improve the predictive ability of the model and to address 21st century forest management approaches.

Figure 1 .
Figure 1.Location of 267 permanent plots measured within the Western Gulf Coastal Plain longleaf pine natural distribution range.

Figure 2 .
Figure 2. Validation of dominant height (Hdom) (a, b), surviving trees per hectare (N) (c, d) and basal area for unthinned stands (BA) (e to h) models based on 30 plots from the dataset used for model evaluation.Observed versus predicted (simulated) values (a, c, e and g) and residuals (predicted-observed) versus stand age (yrs.)relationship of Hdom (b) and N (d), and residuals versus observed values of predicting (f) and projecting (h) BA.Solid line represents linear fit between observed and predicted values and dotted lines for plots (a, c, e and g) correspond to the 1-to-1 relationship.Residuals are presented as a proportion of observed values.

Figure 3 .
Figure 3. Validation of total stem volume outside (VOL OB ) (a, b) and inside (VOL IB ) (c, d) bark, BA after thinning (BA t ) (e, f) and merchantable volume breakdown (VOL m ) (g to j) models based on 30 pots from the model evaluation dataset.Observed versus predicted (simulated) values (a, c, e, g, i) and residuals (predicted-observed) versus observed values of VOL OB (b) VOL IB (d), BA t (f) and VOL m (h, j).Two examples of VOL m outside bark are shown: using d = 10.16 cm and t = 20.32 cm (g, h) and d = 20.32 cm and t = 30.48cm (i, j).Solid line represents linear fit between observed and predicted values and dotted lines for plots a, c, e, g and i correspond to the 1-to-1 relationship.Residuals are presented as a proportion of observed values.

Figure 6 .
Figure 6.Overall simulation validation of survival (N) (a, b), dominant height (H dom ) (c, d), basal area (BA) (e, f) and stem volume outside bark (VOL OB ) (g, h) predictions.Observed versus predicted (simulated) values (a, c, e, g) and residuals (predicted-observed) versus simulation length (yrs.)(b, d, f, h) relationships for unthinned stands if initial age, N and SI are known, using the models to estimate N, H dom , BA and VOL OB for all unthinned plots in the dataset based only on knowing the initial stand age, N and SI.

Table 1 .
Summary statistics and stand characteristics for 267 permanent longleaf pine plots measured (thinned and unthinned).

Table 2 .
Distribution of observations by age, site index and surviving density for 267 permanent longleaf pine plots measured (thinned and unthinned).
*: SI reported for plots not measured at age 50 yrs.was estimated with model presented in this study.

Table 3 .
Parameter estimates and fit statistics of Western Gulf Coastal Plain U.S. longleaf pine plantation growth and yield equations.

Table 4 .
Summary of model evaluation statistics for N, H dom , BA, BA t , VOL, and VOL m estimations based on 30 plots.
N: trees per hectare (ha −1 ); H dom : average total height of dominant and codominant trees (m); BA: predicted basal area of unthinned stands (m 2 •ha −1 ); BA j : projected basal area of unthinned stands (m 2 •ha −1 ); VOL: total stem volume (m 3 •ha −1 ); BA t : basal area of thinned stands (m 2 •ha −1 ); VOL m : merchantable stem volume for trees d cm and above to a t cm top diameter limit (m 3 •ha −1 ); : mean observed value; : mean predicted value; n: number of observations; MAE: mean absolute error; RMSE: root of mean square error; Bias: absolute bias; R 2 : coefficient of determination.OB : outside bark; IB : inside bark; t: top diameter (outside bark) merchantability limit (cm); d: dbh threshold limit (cm).Values in parenthesis are percentage relative to observed mean.

Table 5 .
Summary of overall model evaluation statistics for N, H dom , BA and VOL OB estimations using different reference age for SI for different simulation lengths.
N: trees per hectare (ha −1 ); H dom : average total height of dominant and codominat trees (m); BA: stand basal area (m 2 •ha −1 ); VOL OB : total stem volume outside bark (m 3 •ha −1 ); : mean observed value; : mean predicted value; n: number of observations; MAE: mean absolute error (m); RMSE: root of mean square error (m); Bias: absolute bias (m); R 2 : coefficient of determination.Values of MAE, RSME and Bias are percentage relative to observed mean.