SITE QUALITY CURVES FOR AFRICAN MAHOGANY PLANTATIONS IN BRAZIL

Site quality estimation is an important tool in forest management since it is useful for modeling growth and yield for even-aged stands. Data from African mahogany (Khaya ivorensis A. Chev.) Brazilian plantations were used to develop a model to predict dominant height growth, comparing dynamic base-age invariant site index models with the guide curve method (static models). For the evaluation of the candidate models qualitative and quantitative criteria were used. We also verified the stability of the candidate models, preferring a model providing fewer site class changes when predicting site index from different ages. The Lundqvist-Korf function fitted with the guide curve method proved to be effective and accurate for site classification and dominant height predictions of African mahogany stands. The range of observed site index, at a reference age of 15, was between 17 and 33 meters. CURVAS DE CLASSIFICAÇÃO DE SÍTIO PARA PLANTIOS DE MOGNO AFRICANO NO BRASIL RESUMO: A estimativa da qualidade de sítio é uma importante ferramenta para o manejo florestal, uma vez que auxilia na modelagem do crescimento e produção de florestas equiâneas. Dados de plantios brasileiros de mogno africano (Khaya ivorensis A. Chev.) foram utilizados para avaliar modelos de crescimento em altura dominante comparando modelos dinâmicos invariantes com a idade padrão e modelos estáticos (método da curva guia). Para a avaliação dos modelos foram utilizados critérios quantitativos e qualitativos. Também foi verificada a estabilidade dos modelos candidatos, sendo preferido o modelo que resultou em menos mudanças de classes na predição da mesma parcela em diferentes idades. O modelo de Lundqvist-Korf com o método da curva guia se mostrou eficiente, sendo recomendado para a classificação de sítios e para a predição de altura dominante em povoamentos de mogno africano. A amplitude dos índices de sítio, na idade de referência de 15 anos, foi de 17 a 33 metros.


INTRODUCTION
African mahogany, that includes the species Khaya ivorensis A. Chev., is a tree species that belongs to the botanical family Meliaceae and is well known in the international wood market (PINHEIRO et al., 2011).The recent cultivation of the species outside its native location is relatively new, starting in Brazil in 1976 at Pará state with five seedlings (FALESI; BAENA, 1999).These specimens yielded seeds which were used to start up plantations in other regions of the country (FRANÇA et al., 2016).The definition of an accurate prediction of site productivity is fundamental to predict timber yields and for meaningful simulation studies (VANCLAY, 1992;BRAVO;MONTERO, 2001).This concern is not new and the first land classification system of productivity is attributed to the roman Cato in 234-139 B.C. (TESCH, 1981) and several methods for evaluating forest site productivity have been studied (BURKHART;TOMÉ, 2012).
Many are the models based on the dominant heightage relationship applied to site index curves construction (SCOLFORO; MACHADO, 1988;DIÉGUEZ-ARANDA et al., 2006, DAVID et al., 2015), and those are the basis of most traditional stand-level forest management regimes (RENNOLLS, 1995).Among the techniques available to fit site index curves, most are based on three general methods (CLUTTER et al., 1983): guide curve; parameter prediction, and difference equation methods.Although there are many techniques available for site classification no, research has been published until now for African mahogany plantations in Brazil or in other parts of the world.
The guide curve represents the height development for the average site index in the data.Heights at all ages for all other site classes are typically assumed to be proportional to that of the guide curve (BURKHART; TOMÉ, 2012).
The difference equation approach can be applied to any height-age equation to produce families of anamorphic or polymorphic curves (BURKHART; TOMÉ, 2012).The height (hdom 2 ) at a future age must be expressed as a function of future age (t 2 ), current age (t 1 ) and current height (hdom 1 ), that is: hdom 2 =(t 2 , hdom 1 , t 1 ).The formulation of a growth function as a difference equation can be obtained by substituting one parameter by its value computed from the initial condition (t 1 , hdom 1 ).The substitution of the asymptotic parameter in the growth model produces anamorphic curves with multiple asymptotes and the substitution of any other parameter produces polymorphic curves with a common asymptote (PALAHÍ et al., 2004;TAMARIT-URIAS et al., 2014;ERCANLI et al., 2014, TEWARI et al., 2014).
The objective of this study was to develop a site index model for Khaya ivorensis A. Chev.stands in Brazil comparing the difference equation approach with the guide curve methodology.

MATERIAL AND METHODS
Data were gathered in plantations located in three Brazilian regions (Southeast, Midwest and North) comprising 150 permanent sample plots with different ages and remeasurement times, usually around one year interval (Table 1 and Figure 1).
Forest management and genetic material of the plantations across Brazil are similar, planted with a wide spacing (4x5 to 6x6 for plots with less than 8 years and 12x12 m for older plots), intensive fertilization regimes and in some cases irrigation (areas with precipitation less than 800 mm per year).The plot area varies with the Given the low density of the plantations (the majority of the plots have less than 300 trees per hectare), the traditional definition of dominant height (Assmann's mean height of the 100 thickest trees per hectare) would require height measurement of an excessive number of trees per plot.Also, some of plots have less than 100 trees per hectare, making Assmann's dominant height definition inviable.To avoid these problems, we defined the dominant height as being the mean height of the 30 thickest trees per hectare.Similar strategies have been adopted by Paulo et al. (2011) for Quercus suber stands in Portugal and Danquah (2012) for African mahogany species planted in Ghana, selecting 25 and 40 thickest trees per hectare, respectively, to calculate the average of dominant height.
Several candidate functions to estimate the dominant height (Table 2) were selected, a total of 434 height-age pairs were used to fit the dynamic functions and 584 observations to fit the base functions (static) for the guide curve method.All candidate dynamic functions have the property of being time invariant, ensuring that projections using different starting ages but with equal final ages are equal, as described by Palahí et al. (2004).
The algebraic difference approach (ADA) presents the structure of equations as base-age invariant, and the generalized algebraic difference approach (GADA) improves the traditional ADA, allowing more flexible dynamic equations which can be polymorphic and with multiple asymptotes (TAHAR et al., 2012).When only one parameter of the base model is related to a theoretical site quality measure, GADA is equivalent to ADA (NUNES et al., 2011).
The derivation of difference equation follows the methodology proposed by Amaro et al. (1998) for Richards and Lundqvist-Korf functions, consisting in equating subdefined ratios of base equations for the ADA method (BAILEY; CLUTTER, 1974).GADA (CIESZEWSKI; BAILEY, 2000) method was also used.The McDill-Amateis function (McDILL;AMATEIS, 1992) has the Hossfeld IV as integral form and has been used with success in several growth studies (WANG;PAYANDEH, 1996).
The parameter estimation was done using the nls function of the statistical environment R (R Core Team, 2015), performing a nonlinear regression analysis using the Gauss Newton algorithm.The model selection was based in qualitative and quantitative criteria (goodnessof-fit statistics).The root mean square error (RMSE) [1] in meters, the relative error in predictions (RE) [2] in percentage and the bias ( e ) [3] were calculated as follows, where: p is the number of parameters in the model, n is the number of observations used to fit the model, y i is an observed dominant height, ˆi y is the LG A,k Rk Richards ( ) ) corresponding predicted dominant height and is the average of the observed dominant height.

RESULTS AND DISCUSSION
The model with the best performance using the guide curve method (integral form) was the Lundqvist-Korf.For the dynamic functions the models Lk, LG and MD had the best statistics (Table 3).The nonlinear regression analysis using the Richards dynamic function either did not converge (Rk) or didn't have all the parameters significantly different from zero (RA), exception for model Rm with the m parameter as site specific.
The Lundqvist-Korf model had a rather high asymptote (A=58 m), more realistic values were found in equations Lk and Richards (A≈35 m), although Lemmens (2008) reported growth for K. ivorensis reaching 60 meters tall.Considering the goodness-of-fit statistics and the behavior of the asymptotic values of the model parameters we chose Richards, Lundqvist-Korf, McDill-Amateis for the guide curve method and Lk, LG and MD for the dynamic base-age invariant as the best candidate equations for site classification and discarded the others from further analysis.
Of the six selected models, the models fitted with the guide curve method presented a tendency towards heteroscedasticity (Figure 2).However, only the Lundqvist-Korf model did not present some tendency to underestimate values for smaller dominant height estimates (below 7 m) that could be observed for the other candidate models.
Before proceeding to the actual site classification, we plotted the site class limits of the best models to evaluate their behavior (Figure 3). ) For the qualitative criteria an evaluation of the different candidate models involved visual examination of the residuals against the estimated values and the fitted curves for different site index classes overlaid on the profile plots of the stands.Visual inspection is an essential point in selecting the most accurate model because curve profiles may differ considerably, even though goodness-of-fit statistics are similar (DIÉGUEZ-ARANDA et al., 2006).The models were chosen according to the goodness-of-fit, predictive ability, biological sense and compliance with the assumptions of nonlinear least square method (homoscedasticity, lack of autocorrelation and normality of residuals).
Due to the absence of continuous forest inventory data till the final rotation and the recent implementation of African mahogany culture in Brazil, an independent database for validation is not available, so the consistency of the models was analyzed by observing the stability of the site classification at the plot level, preferring classifications that were more stable, with fewer class changes for a single plot throughout the remeasurements.While this approach has been applied before (e.g.SCOLFORO, 1992;SELLE et al., 1994;MACHADO et al., 2011) it has been used to determine if the selected model produced an adequate classification using qualitative criteria.Here we propose a quantitative test to enable the comparison of the stability of the site classification generated by different models.This is done by quantifying the number of times each permanent plot changed site index over its life span.We then attributed weights for the number of changes with lower values for fewer changes (0 for no changes, 1 for one change and so on).The sum of the number of changes multiplied by this weight permitted the comparison of different classifications, identifying the model that presented the most stable classification as being the one with the lowest weighted sum.Although the values of the goodness of fi t for functions Lk and MD were good, the curves showed unrealistic shapes.For both models even fi xing the asymptote parameter still yielded curves with unrealistic classifi cations, but with better shape (Figure 4).2004) modelling site index for Pinus sylvestris in Spain also found that the McDill-Amateis model (MD) did not perform well resulting in poor estimates when extrapolated outside the age and site index range data, showing a similar curve shape as we found in our study to the highest site index before fi xing the asymptote parameter.Castro et al. (2015) evaluated the efficiency and reliability of two methods, guide curve and differences equations, concluding that the model developed by the difference equations method was more efficient to perform the prognosis of the studied Eucalyptus stand in Brazil.Arias-Rodil et al. (2015) comparing the agedependent method with age-independent alternatives for estimating site index curves for Maritime Pine in Spain found that the age-independent equation performed best.
To check if our selected models presented autocorrelation, we used the plotted residuals versus lag residuals (Figure 5).The Lundqvist-Korf and Richards equations fitted with the guide curve method present a positive autocorrelation, confirmed by the significant values of the Durbin Watson test (d= 0.797 for Lundqvist-Korf and d= 0.815 for Richards).Equation LG did not show a strong temporal autocorrelation (d= 2.197) at P<0.05.The Q-Q plot of the standardized residuals for all equations showed that the residuals follow approximately the normal distribution.Due to the characteristics of the modelling data, the assumption of the independence of the error term is likely to be violated (PALAHÍ et al., 2004), however, hypothesis test may be quite unnecessary when the form of the model being fitted is limited, as the case of the traditional Lundqvist-Korf and Richards models.
We are aware of the restrictions of our data, since few continuous forest inventories are available in Brazil and only three plots with measurements over the age of 10 compose the data base.Model selection was viewed as a compromise between statistical and biological considerations (reasonable values for prediction of dominant height for future ages as presented in Table 6), so the Lundqvist-Korf in its integral form is more For the final decision of the best model for site classification the classification stability (Table 5) was analyzed for the three more realistic models (Lundqvist-Korf and Richards using the guide curve method and LG using dynamic base-age invariant method).
The Lundqvist-Korf function was superior for the stability of the site classification.Model LG is a dynamic site index equation with 2 site-specific parameters (a representing the asymptote and b the inflexion point), providing two important and desirable attributes for site index classification, polymorphism and multiple asymptotes (CIESZEWSKI, 2002).In developing site index prediction systems, some desirable characteristics of site index models are polymorphism, multiple horizontal asymptote, one inflection point, and base-age invariance (ELFVING; KIVISTE, 1997;CIESZEWSKI;BAILEY, 2000).indicated for site classifi cation, but the polymorphic equation derived by the GADA approach from Lundqvist-Korf function (LG) resulted in a more adequate model for prediction of dominant height.Although advantages of dynamic site equations have been pointed out in several studies (AMARO et al., 1998;CALAMA et al., 2004;LOPEZ-SENESPLEDA et al., 2014) we decided to choose the more parsimonious model of Lundqvist-Korf using the guide curve method for site classifi cation.
The site quality curves obtained in this study clearly show that Brazil is an excellent territory for African mahogany growth.The mean annual diameter increment for the different site quality classes (Figure IV ( 21) III ( 25) II ( 29) I ( 33) 1 1.6 -2.0 2.0 -2.5 2.5 -2.9 2.9 -3.3 3.3 -3.8 5 7.7 -9.8 9.8 -11.8 11.8 -13.9 13.9 -15.9 15.9 -18.0  4 -20.7 20.7 -25.1 25.1 -29.4 29.4 -33.8 33.8 -38.2 25 17.1 -21.6 21.6 -26.2 26.2 -30.7 30.7 -35.3 35.3 -39.8 30 17.4 -22.1 22.1 -26.7 26.7 -31.4 31.4 -36.0 36.0 -40.7 35 17.6 -22.3 22.3 -27.0 27.0 -31.7 31.7 -36.4 36.4 -41.1 6) showed high values compared with other places that plant African mahogany, e.g.Heryati et al. (2011) comparing the growth of K. ivorensis A. Chev. in three different soil series in Malaysia reported values of mean annual increment for diameter from 2.32 to 2.88 cm year -1 .Erskine et al. (2005) conducted a trial with exotic species for fi ve years in Australia, including K. nyasica, related a mean annual increment in diameter for Khaya of 1.58 cm .year -1 .A more detailed study on quality of the proposed site classifi cation should be carried out at the end of the fi rst rotation of the plantations, with the aim of confi rming whether the chosen model and reference age limits remain the most appropriate for the Khaya ivorensis A. Chev.plantations.Using the model selected on this study (Lundqvist-Korf integral form) for the fi ve site indices, we calculated the limits for the dominant height by age and classes for the different site classes (Table 6).

CONCLUSIONS
For site classifi cation purposes, we found that the Lundqvist-Korf model, using the guide curve method, is the one with the most stable classifi cation.
[4] [5] This model leads to the following equation for S prediction.
This was confi rmed by the higher mean annual increment in diameter for the higher productivity classes.

FIGURE 2
FIGURE 2 Standardized residuals plotted over the fi tted values for the best equations to estimate dominant height.

FIGURE 3
FIGURE 3Comparison of height curves produced by the best models selected for site classifi cation of African mahogany plantations in BrazilPalahí et al. (2004) modelling site index for Pinus sylvestris in Spain also found that the McDill-Amateis model (MD) did not perform well resulting in poor estimates when extrapolated outside the age and site index range data, showing a similar curve shape as we found in our study to the highest site index before fi xing the asymptote parameter.

FIGURE 4
FIGURE 4 Height curves for equation Lk and MD with fixed asymptote parameter (A=45 and A=40, respectively).

FIGURE 5
FIGURE 5 Autocorrelation and Q-Q plot for the selected equation Lundqvist-Korf (a), Richards (b) and LG (c).

FIGURE 6
FIGURE 6 Mean annual increment in diameter of plantation grown Khaya ivorensis A. Chev. in Brazil considering different site index classes obtained by the Lundqvist-Korf integral form model.

TABLE 1
Characteristics of the sample plots of African mahogany plantations used in the study.Observed dominant height development of Khaya ivorensis A. Chev.stands in Brazil spacing used (800 to 4,400 m²).The oldest plot (in the North region) is composed of three trees planted in 1976 (about 20 meters from each other), so the dominant height was computed as the average of those trees.It was decided to keep these trees in the data base because they are the oldest and highest African mahogany trees in Brazil, thus allowing a more realistic estimation of the asymptote parameter of the models.

TABLE 2
Candidate models for dominant height estimative in function of age for Khaya ivorensis A. Chev.stands in Brazil.

TABLE 3
Goodness-of-fit statistics and parameter estimations of candidate functions, where ns indicates a nonsignificant parameter at 95% of probability.

TABLE 5
Classification stability of the best models selected for site classification.

TABLE 6
Dominant height estimative limits in meters by age and site index classes for Khaya ivorensis A. Chev. in Brazil.