SITE INDEX IN EUCALYPTUS STANDS APPLYING ORDINARY KRIGING: AN APPROACH WITH DIFFERENT MODELS AND METHODS OF CLASSIFICATION

The objective of this work was to evaluate the efficiency of models and methods to obtain the site index, associated with ordinary kriging, to classify productive capacity in eucalyptus stands. Thus, the site quality was performed considering the traditional modeling in clonal stands (2,119 hectares) located in Minas Gerais state, Brazil. 170 plots of 400m were randomly allocated, representing a sampling intensity of 0.32%. The dominant height of trees (Assmann) was measured at 24, 36, 48, 60, 72, and 84 months. The site index (S) was estimated by the guide curve and algebraic difference methods, using the models of Schumacher, Chapman and Richards, and Bailey and Clutter. 136 plots were used in the fit and 34 plots in the predictive validation. The spatial dependence of site index was evaluated by experimental semivariogram and adjustment of exponential, spherical, and gaussian models. After confirming the spatial dependence, ordinary kriging was performed to spatialize the site index. For the predictive validation, the dominant height values at 72 months were used. The algebraic difference method provided excellent estimates of site index, which showed spatial dependence in all adjustments, from moderate to strong. In most cases, the gaussian model was the most accurate. It is concluded that the algebraic difference method was more efficient and the site index showed strong spatial dependence at all ages, regardless of the model used. Thus, regression models for site index estimation can be used in combination with ordinary kriging techniques.


INTRODUCTION
The productive capacity refers to the productive potential of wood of a certain area and species, and it can be performed using indirect or direct classification methods (CAMPOS; LEITE, 2017). Direct methods make it possible to represent the productive potential through a quantitative value, called the site index. The site index is a robust indicator of site quality, used to help apply management techniques and assess the economic viability of forest plantations, in addition to being a fundamental variable in growth and production models (LEITE et al., 2011;CAMPOS;LEITE, 2017 The application of direct methods to classify productive capacity is considered efficient (LEITE et al., 2011), based on the use of regression models that relate dendrometric variables with the age of the population. Among the variables that can be used as an indicator of the site quality, the volume, basal area, total height, and dominant height of eucalyptus stand out (SCOLFORO, 2006). However, the dominant height is considered the variable that provides the most consistent results once it is low influenced by the density of the population and presents a significant correlation with the volume of wood.
The methods generally used to classify the productive capacity do not consider the spatial dependence of the dominant height variable, not even the spatialization of the site index. Some studies have confirmed the strong spatial dependence of the dominant height in commercial stands of Eucalyptus sp. (MELLO et al., 2005;GUEDES et al., 2015;LUNDGREN et al., 2017;ATAÍDE et al., 2020). In addition, the regression models provide the average estimate of growth in dominant height and allow site index estimates only in places where there was a measured sample unit. In this case, for a more detailed mapping of the productive capacity, it would be necessary to increase the sampling intensity to represent different locations, increasing sampling costs.
Ordinary kriging is the option commonly used in the spatialization of dendrometric variables in nonsampled locations. It can be cited studies carried out with diameter at 1.30 m from the ground, basal area, the total and dominant height, volume, and average annual increment, being applied for Eucalyptus sp. stands, as well as to other species, for example, Pinus sp. and Tectona grandis (MELLO et al., 2005;PELISSARI et al., 2014;GUEDES et al., 2015;PELISSARI et al., 2015;LUNDGREN et al., 2017;ZECH et al., 2018;ATAÍDE et al., 2020). Even with numerous studies evaluating the spatial structuring of dendrometric variables, few use the site index variable, mainly when obtained by mathematical models. Thus, the results this variable obtained in different adjustment methods and models are still unknown in the spatial estimates of the variable.
Considering that the dominant height has proven spatial dependence, it is possible to assume that the site index will show strong spatial dependence, as the site index is the dominant height at the reference age. If accepted, this hypothesis would make it possible to spatialize the classes of productive potential in forest stands, generating relevant spatial information for management. As there are different procedures to estimate the site index and different geostatistical models to obtain spatial estimates of a variable, it is important to evaluate which ones provide satisfactory results in estimating the site index and spatialization of the productive potential classes.
This work hypothesizes that the choice of the regression model and the method for constructing the site curves influence the spatial continuity of the site index and the spatial estimates through ordinary kriging. Therefore, the objective of this work was to evaluate the effect of different models and methods on the spatial dependence structure of the site index and the application of ordinary kriging to spatialize the productive capacity in eucalyptus stands.

Study area
The study was carried out in 2,119 hectares with clonal Eucalyptus sp. stands, located at coordinates 17º17'26.9" S, 43º42'33.7" W, in the municipality of Bocaiúva, Minas Gerais state. The local climate is of the Aw type, according to the Köppen climate classification, characterized as being tropical wet savanna, dry winters, and rainy summers (ALVARES et al., 2013). The region has an average altitude of 820 m, with rainfall and an average annual temperature of 1,246 mm and 24ºC, respectively. The predominant soils are Dystrophic Dark Red Latosol or Dystrophic Red Yellow Latosol.

Data collection and description
Information was obtained from plots allocated in 62 blocks of the stands, with the same genetic material and initial spacing. 170 square plots were selected with a fixed area of 400m², representing a sampling intensity of 0.32% of the total area (one plot for every 12.46 hectares). The distances between the plot ranged from approximately 200 to 44,354 m. In all plots, the dominant heights of trees (Hd) were obtained following the concept of Assmann (1970), at ages 24, 36, 48, 60, 72, and 84 months. The plots were georeferenced based on their central geographic coordinates. The plots were divided into two groups: 136 (1 plot for every 15.58 hectares) were used to adjust the models and evaluate the productive capacity classification procedures; 34 (1 plot for every 62.32 hectares) were used for the predictive validation of the estimates of the best selected procedures. Validation plots were selected throughout the study area, considering at least one validation plot per stand, when possible.

Classification of the productivity capacity
The site index was estimated considering the guide curve and algebraic difference methods, with a reference age of 72 months. In each method, three regression models were adjusted to estimate the dominant height as a function of age (Table 1). The adjustments were realized using the nlstools package (BATY et al., 2015) (R CORE TEAM, 2015). Six equations were obtained to estimate the site index (Table 1) after rearranging the three fitted models, using the guide curve and algebraic difference methods.
The adjusted equations were evaluated using the adjusted coefficient of determination (R²adj.), standard error of the estimate in percentage (Syx%), root mean square error (RMSE), the mean error of prediction in percentage (bias%), and graphical analysis of the normalized residuals. The significance of the coefficients was also verified by the t test and the normality of the residuals using the Shapiro-Wilk test, both with 5% probability (α=0.05). Table 1. Models used to estimate site index through guide curve and algebraic difference methods. Tabela 1. Modelos utilizados para estimar índice de sítio pelos métodos da curva-guia e diferença algébrica.

Original model Author Model by the Guide Curve method
Model by the Algebraic Difference method Subtitle: Ln = neperian logarithm; Hd1 = dominant height of trees (m); β0, β1, β2 and β3= model parameters; I = age of population (months); S the site index (m); Iref = the reference age (months); e = exponential; εi = residue.
The spatial continuity structure of the site index, estimated in all combinations of methods and models, was evaluated using the experimental semivariogram. The anisotropy was evaluated by constructing directional semivariograms (0, 45, 90, and 135º), which confirmed the isotropic behavior of the variable. The exponential, spherical, and gaussian models were fitted to the omnidirectional experimental semivariogram by the Maximum Likelihood Method, using the geoR package (RIBEIRO JÚNIOR; DIGGLE, 2001) in the R software (R CORE TEAM, 2015).

Gaussian model
Where: ̂(ℎ) the semivariance, ( ) the value of the regionalized variable at point x, ( + ℎ) the value of the variable at point, + ℎ, (ℎ) the number of pairs separated by the distance h, 0 the nugget effect, the contribution, the reach.
The criteria used to select the best semivariance model were: Akaike's information criterion (AIC), reduced mean error (RME) and standard deviation of the reduced error (SDRE) obtained by cross-validation (ARAÚJO et al., 2019). The spatial dependence index (SDI) was obtained for each of the theoretical models adjusted to semivariograms, classified as weak (SDI ≤ 0.25), moderate (0.25 < SDI ≤ 0.75), and strong (IDE > 0.75) spatial dependence (ZIMBACK, 2003). Once the spatial dependence of the site index was confirmed, ordinary kriging was applied to obtain the spatial estimates of the site index in the non-sampled locations.
Ordinary kriging was used to predict the spatial estimates of site index obtained in each regression model adjusted by the guide curve and algebraic difference methods. The spatial estimates of the ordinary kriging in the Where: SDI the spatial dependence index, C the contribution, C0 the nugget effect, c the performance index, r the Pearson correlation coefficient and d the Willmott agreement index.

Predictive validation
The 34 plots not used in the fit of the regression models were used to perform the predictive validation of the site index estimated from the regression models and ordinary kriging. The site index estimates by regression and ordinary kriging were compared with the observed values of the mean height of the dominant trees in the validation plots, at the reference age of 72 months. In ordinary kriging, the maps were converted to raster format with pixel dimensions of 20 m x 20 m, presenting an area equal to the size of the plot. In the maps, the 34 plots separated for validation were inserted and, later, the estimated values of the site index were extracted. The extracted values were compared with the dominant height of the trees at the reference age of 72 months, through the difference between the dominant height observed in the plot and the limit of the respective class where it was positioned.

Adjustment of models for classification of productive capacity
The coefficients of the regression models were significant by the t-test (p < 0.05) for all adjustments made ( Table 2). The adjusted coefficient of determination (R adj. 2 ) was higher than 0.80 for the two methods of evaluation of the productive capacity. However, the other statistical criteria showed different results between the methods and models, indicating the superiority of the Bailey and Clutter model in the algebraic difference method. The standard error of the estimate in percentage (Syx%) was less than 10% in all adjustments, but with values less than 7% only when the algebraic difference method was used. A similar result was obtained for the RMSE estimates, in which the algebraic difference method provided improvements of about 25% concerning the guide curve, with values around 1.5 m. Regarding the bias% there was a distinct trend: while the guide curve method tended to overestimate, the algebraic difference provided underestimates. In the graphical analysis of the normalized residues (Figure 1), it is observed that the models fitted in the algebraic difference method showed greater homogeneity in the dispersion of residues.

Spatial analysis of the site index
The site index presented spatial dependence ranging from moderate (0.25 < SDI ≤ 0.75) to strong (SDI > 0.75) in all regression models used (Schumacher, Chapman and Richards, Bailey and Clutter) in the two methods of construction of site curves (guide curve and algebraic difference) ( Table 3). The SDI ranged from 63.42% to 90.75%, indicating a predominantly strong spatial dependence of the site index variable. The spherical, exponential, and gaussian models presented good fits to the experimental semivariogram, in all the productive capacity classification procedures, with RME values close to 0 and SDRE close to 1. Subtitle: C0 the nugget effect; C the contribution; A the reach; RME the reduced mean error; SDRE the standard deviation of the reduced error; AIC Akaike's information criteria; SDI% the spatial dependence index in percentage.
To select the best models at each age, the lowest AIC value was used as the main selection criteria. The models are considered statistically different if the difference between the AIC values is greater than 2. Thus, it was found that the gaussian model presented the best fits in most ages and productive capacity classification procedures (combination of 3 models with 2 methods). Exceptions occurred for the adjustments made at the age of 48 months, where the exponential model stood out as the best, and at 24 months using the Schumacher model by the guide curve method, in which the spherical semivariance model presented the best fit.
The nugget effect values ranged from 1.11 to 3.35, while the contribution from 3.21 to 16.81, and the range of spatial dependence was between 968 to 4,888. The highest contribution values in relation to the nugget effect prove the spatial dependence of the variable, making it possible to apply ordinary kriging to spatialize the site index in the non-sampled locations. Figure 2 shows the comparison statistics between the estimated site index values using the adjusted regression equations and applying the ordinary kriging. The RMSE decreased exponentially from 24 to 84 months, indicating an increase in the accuracy of ordinary kriging in estimating the site index over age. The bias% showed negative values, demonstrating a mean tendency to overestimate the site index applying ordinary kriging. Pearson's correlation values, Willmott's agreement index, and performance showed a tendency to decrease up to 48 months and, from that age onwards, they started to increase, for the different procedures to classify the productive capacity.

Spatialization and validation of productive capacity
The site index estimated at 72 months, based on the algebraic difference method and ordinary kriging, for the three adjusted regression models were compared with the dominant height observed at 72 months, in the 34 plots of validation (Figure 3). The site index estimated by regression models in the algebraic difference method showed greater precision and accuracy than those estimated via ordinary kriging. Determination and angular coefficients greater than 0.90 demonstrate the predictive power of these models. In the spatial estimates of the site index, applying the ordinary kriging, these coefficients were reduced to values between 0.6 and 0.7, demonstrating that in the spatialization process there is a reduction in the accuracy.  Figure 4 shows the spatialization through ordinary kriging using 3 and 7 site index classes (amplitude of 5 and 2 m, respectively) obtained by the models of Schumacher, Chapman and Richards, and Bailey and Clutter, in the method of the algebraic difference. The model by Chapman and Richards showed greater accuracy, with 20 plots (58.83%) classified in their respective classes and 14 plots (41.17%) classified in other classes. The Schumacher and Bailey and Clutter models showed lower accuracy, with 15 plots (44.11%) and 16 plots (47.05%), respectively, classified outside the correct site classes. However, the greatest difference found among all models was 3.72 m, equivalent to 14.54% of the observed average. Most of the errors obtained were not greater than 1 m (3.90%), demonstrating that even using smaller amplitude classes, the spatialization of the site index is performed with high accuracy.
By increasing the amplitude of the class and reducing the number of classes to 3, the following classification was established for the classes: good (Class I) ranging from 27 to 32 m, regular (Class II) from 22 to 27 m and poor (Class III) from 18 to 22 m. Under these conditions, the accuracy of the productive capacity mapping increases considerably for all models, as 28 plots of validation (82.36%) fit the correct site class and only 6 plots (17.64%) did not fit, but with errors less than 1 m.
Considering the same class amplitude, the differences between the ordinary kriging maps obtained by the Schumacher, Chapman and Richards, and Bailey and Clutter models are imperceptible, indicating that the site mapping is little influenced by the regression model used. In this case, the choice of the method of construction of the site curves will have a more significant influence, with superiority for the algebraic difference method.

DISCUSSION
The Schumacher, Chapman and Richards, and Bailey and Clutter models presented satisfactory fits, with superiority of the algebraic difference method compared to the guide curve method. The superiority of the algebraic difference method in the construction of site index curves was also observed by Castro et al. (2015), in a eucalyptus stand in Minas Gerais, with the standard error values of the estimate close to 6%. Pissinin and Schneider (2017) observed that the algebraic difference method presented the smallest estimation errors compared to the guide curve, using the Chapman and Richards model. Kitikidou et al. (2011), applying the algebraic difference method to construct site index curves in a Pinus brutia stand, observed that the Bailey and Clutter model was superior to Chapman and Richards, corroborating the results obtained in this study. In the present work, the Bailey and Clutter model in the algebraic difference method presented the best fit statistics, with values (in module) of Syx%, RQEM, and bias% of 6.363, 1.479, and 0.158, respectively. The superiority of the algebraic difference method was confirmed from the graphic analysis of the normalized residuals, which presented a more homogeneous distribution compared to the guide curve method (Figure 2). A similar characteristic was found by Pissinin and Schneider (2017) when they adjusted different models in the guide curve and algebraic difference methods in eucalyptus stands, in which the algebraic difference method presented a homogeneous distribution of residuals and smaller errors.
In the spatial structuring of the site index, there was moderate to strong spatial dependence in all analyzed scenarios, with SDI values above 60% (Table 4). This is due to the strong spatial structure of the dominant height in species of the genus Eucalyptus sp. (MELLO et al., 2005;GUEDES et al., 2015;ATAÍDE et al., 2020). This characteristic was also observed in stands of other forest species, such as Tectona grandis (PELISSARI et al., 2015) and Pinus sp. (ZECH et al., 2018). The most accentuated nugget effect values were obtained in the initial growth phase of the stands, between 24 and 48 months, ages in which a greater variability was observed in the