Approaches to Forest Site Classification as an Indicator of Teak Volume Production

Abstract

We conducted a study on the dominant height growth of clonal teak (Tectona grandis Linn F.) plantations in the Brazilian Amazon to assess their potential and its agreement with volumetric production. We employed two approaches, ADA (algebraic difference) and GADA (generalized algebraic difference), and analyzed data from 58 permanent plots collected over a 10 year period. To classify the sites, we developed equations and evaluated their accuracy using various criteria, including correlation coefficient, mean square of residual, Akaike’s criterion, distribution of residuals, and validation through equivalence testing (TOST). We also assessed the biological realism of the constructed curves. We used cluster evaluation and dendrogram comparison to assess the agreement between site index and volumetric production for each approach. The Lundqvist–Korf baseline models (M1–ADA and M4–GADA) proved to be accurate and realistic in estimating dominant height in both approaches. Our findings indicate that the approaches utilizing dynamic equations and generating polymorphic curves effectively represent the sites and indicate the volumetric production of the plantations, with 98.3% of agreement rate. Based on our results, we recommend the use of ADA and GADA approaches for estimating the dominant height of clonal teak plantations in the Eastern Brazilian Amazon.

1. Introduction

Teak (Tectona grandis Linn F.) is one of the most cultivated timber forest species in the world, which is due to its high commercial timber value, given its unique mechanical, structural, and aesthetic qualities. This species is originating from several tropical and subtropical climate Asian countries, such as India, Myanmar, Laos, and Thailand, where natural forest stands still exist. Furthermore, teak is currently observed in Asia, Africa, and Latin America, with more than 4 million hectares of planted stands [1,2,3]. In Brazil, it is among the most emerging cultivated tree species, and its cultivation in planted areas has grown substantially, along with other species, such as Schizolobium parahyba var. amazonicum (Huber × Ducke) Barneby and Acacia mangium [4]. Due to the characteristics of its wood, teak plantations are considered an important source of high-quality wood, which has contributed to mitigating human pressure on native forests [5]. The worldwide expansion of this tree species indicates teak’s high ecological versatility, especially regarding its adaptability to diverse soil and climate conditions [6].

Despite the great aptitude that teak has for large-scale timber production, it is crucial to adopt advanced silvicultural management tools and techniques aiming to ensure successful establishment, growth, and production [7]. In this regard, information on the silvicultural behavior of this tree species, coupled with advanced mathematical modeling techniques, has been extremely important for strategies to optimize teak production [8,9]. It requires an advancement of those techniques, especially for clonal stands of these tree species, given the greater potential of those stands to leverage productivity [10] compared to seminal stands.

In this context, the site classification through dominant height modeling, aiming to determine the productive capacity of teak stands, has been one of the most promising tools to support management practices for this species [11,12]. There are several available techniques to optimize site index curves, but most of them are based on the following methods [13]: Guide curve, parameter prediction, algebraic difference approach (ADA), and generalized algebraic difference approach (GADA) [14].

In Brazil, commercial forests occupy about 10 million hectares, and these forests are responsible for the production of a significant amount of wood, with part of this production intended for export and being able to meet the international market demand [4]. Given the crucial role that Brazilian forests play, it is essential to determine an accurate and realistic method to assess the productive capacity of these areas. Although the guide curve method is still widely used in the country for several species and management practices, such as Eucalyptus spp. [15], Pinus spp. [16], Khaya spp. [17], and teak [18], it has some restrictions and shortcomings that may affect the accuracy of projections in volumetric production of Brazilian planted forests, as reported by [19].

As an alternative to the guide curve method, the authors of [20] proposed an algebraic difference-based approach (ADA). This method consists of deriving dynamic equations, allowing for the substitution of equation parameters and generation of anamorphic or polymorphic curves. Polymorphic curves can describe the dominant height growth of hardwood forest species, such as teak, which, as a rule, presents dynamic growth that changes in response to the environment and silvicultural practices [21]. Several studies have used ADA as a way to determine the productive capacity of forest stands [22]. Despite the benefits offered by this method, there is the possibility of obtaining curves with one single asymptote, depending on the parameters replaced in the equation. This may lead to inconsistencies in the dominant height estimation [23].

Some studies have used an adaptation of the ADA method, which consists of a generalization of the algebraic difference approach (GADA). This adaptation of the ADA method aims to predict polymorphic curves more efficiently, using more than one parameter and multiple asymptotes, which can be based on characteristics of the environment where the forest stand is located [23,24,25]. For clonal teak stands located in tropical regions, such as the Eastern Amazon, the use of the GADA-based approach can be an interesting management tool to make a prediction of volumetric production [26] and/or to support thinning practices [12] with ADA, for example.

The comparison among different approaches to determine the productive capacity of teak is indispensable in forestry studies. This comparison should take into account not only the accuracy of the estimates, the biological realism of the dominant height growth curves, and their generated productive capacity classes, but also their indicative yield. Therefore, the estimates performed by each site classification approach may generate a distinct degree of responsiveness in volumetric production. They affect the projections of forest growth and production models since height is strongly correlated with potential volume growth [27].

This study investigates the relationship between dominant height growth, productive potential, and volumetric production in teak clonal plantations. Field data collected from permanent plots in teak forest stands in the Amazon region are used for this analysis. The main research question is whether the ADA (algebraic difference) and GADA (generalized algebraic difference) approaches are suitable for accurately and realistically representing the relationship between the productive capacity (site) and volumetric production in teak clonal plantings. Our hypothesis proposes that dynamic equations employed in these approaches offer an accurate and realistic representation of the productive potential, thereby providing an indication of the volumetric yield. The objectives of this study were to determine and compare the productive capacity of teak clonal plantations in the Eastern Brazilian Amazon using the ADA and GADA approaches. Additionally, this study aimed to evaluate the agreement of each approach with the volumetric production to investigate their effectiveness in confirming the yield.

2. Materials and Methods

2.1. Study Area and Silvicultural Practices

We conducted the present study in clonal teak plantations, located in the municipality of Capitão Poço, State of Pará, Eastern Brazilian Amazon (central coordinates: 2°30′0″ S; 47°20′00″ W and 2°20′0″ S; 47°30′0″ W) (Figure 1). The region presents a slightly undulating and flat relief, originally covered by Dense Ombrophylous Forest [28]. The predominant soils are Petroplintic Dystrophic Yellow Latosol, Typical Dystrophic Yellow Latosol, and Petric Concretionary Plintosoil [29]. According to the Köppen classification, the climate is Am type and it is characterized as hot and humid, with a short dry season [30]. The average annual rainfall and temperature are 2256 mm and 26.1 °C, respectively [31].

The monoclonal teak stands are between 7 and 12 years old and were manually planted at 3.5 × 3.5 m, 3.75 × 3.75 m, and 4 × 4 m spacings. The same silvicultural practices were applied in all plantings established between 2010 and 2015, as well as in the sampling plots established in these plantings, in identical calendars. These silvicultural practices consisted of combating leaf-cutting ants with ant bait; cleaning the area with a bulldozer; liming with dolomitic limestone (3 t ha⁻¹); planting fertilization in the planting hole (200 g plant⁻¹ of NPK 8-28-16 and 100 g plant⁻¹ of KCl); control of weeds through crowning with hoe; mechanized and semi-mechanized weeding with a hydraulic tractor; maintenance fertilization with application of Boron (7 g plant⁻¹) and KCl (100 g plant⁻¹) and artificial pruning with a saw and motor pruner [12]. Thinning was systematically performed at 4.5 (1st) and 8.5 (2nd) years, reducing the basal area by 50% in both interventions [12,26].

2.2. Forest Inventory

We performed periodic inventories for 10 years in 58 circular plots (500 m²), using the systematic sampling process. The plots were distributed on a regular grid of 320 m × 320 m (Figure 1). The measured variables were: diameter at 1.3 m from the ground (dbh) measured with a diametric tape and total height (h) with a Vertex IV hypsometer [32]. We determined the dominant height (Dh) according to the criteria established by [33], which consisted of the average height of the 100 trees with the highest dbh per hectare; therefore, measuring five dominant trees per plot. We also calculated the volume of a single tree with bark (v) using the Takata model (Equation (1)), constructed for the same stand by [12,34], and then extrapolated to the equivalent measurement of one hectare (m³·ha⁻¹).

$$v=\frac{{\mathit{dbh}}^{2}h}{(\mathrm{20,510.8}+286.7\mathit{dbh})}$$

(1)

where h: total height (m); dbh: diameter at 1.3 m from the ground (cm); v: volume of single tree with bark (m³).

2.3. Modeling the Production Capacity

Site indices (SI) were obtained using the algebraic difference approach (ADA) and generalized algebraic difference approach (GADA), with the selection of the Lundqvist–Korf (M1), Chapman–Richards (M2), and Hossfeld (M3) models for ADA and Lundqvist–Korf (M4), also known as the Bailey–Clutter model (which is widely used in forestry modeling), Cieszewski (M5), and Schumacher (M6) for GADA [7,23,24,35] (Table 1). These models were fitted to data pairs of current Dh (Y₀) as a function of current age (t₀) and future Dh (Y) as a function of future age (t₁). These constructed models allowed us to generate anamorphic (M2) and polymorphic (M1, M3, M4, M5, and M6) curves, which can accurately represent the growth of major forest variables [7]. Candidate dynamic functions have the property of temporal invariance, which means that projections that use different initial ages but have the same final ages are equivalent [36].

The generalized algebraic differences approach (GADA) is an improvement of the traditional ADA method, which allows for more flexibility of dynamic equations, making them polymorphic and with multiple asymptotes. When only one parameter of the base model is associated with a theoretical measurement of site quality, GADA is equivalent to ADA [37]. The development of GADA models is based on modifying parameters of the baseline model by explicit functions of X, and it can be constructed to account for the expected temporal correlation for longitudinal data of dominant height evolution [7]. It consists of an unobservable independent variable, in which it can describe the bioedaphoclimatic factors from the site [26]. As defined by [26] and [12], the plantations were stratified into three productive capacity classes, namely, low (class 3), medium (class 2), and high (class 1) productivities, with a reference age (RA) of 12 years, consisting of the age closest to the rotation age [38], which, in accordance with [26], can range between 14 and 17 years in the Eastern Amazon.

Table 1. Dominant height growth models selected to fit data from teak (Tectona grandis Linn F.) clonal plantations in the Eastern Brazilian Amazon.

Model No.	Base Model	Parameters Related to the Target Variable X	Initial Solution for X with Y0 and t0	Dynamic Equation
M1	Lundqvist–Korf ${Y=A\mathit{exp}}^{\left(-B{t}^{-\Gamma }\right)}$	$B=X$	$X_0=-Ln\left(\frac{Y_0}{b_1}\right){t_0}^{b_3}$	${Y=b}_1{\mathit{exp}}^{\left[\mathit{Ln}\left(\frac{Y_0}{b_1}\right){\left(\frac{t_0}{t_1}\right)}^{b_3}\right]}$
M2	Chapman–Richards $Y=A[1-{{\mathit{exp}}^{(-Bt)}]}^{\Gamma }$	$A=X$	$X_0=\frac{Y_0}{{[1-{\mathit{exp}}^{\left({-b}_1{t}_0\right)}]}^{b_2}}$	${Y=Y}_0{\left[\frac{1-{\mathit{exp}}^{(-b_1{t}_1)}}{1-{\mathit{exp}}^{(-b_1{t}_0)}}\right]}^{b_2}$
M3	Hossfeld $Y=\frac{A}{1+B{t}^{-\Gamma }}$	$B=X$	$X_0={t_0}^{-b_3}\left(\frac{b_1}{Y_0}-1\right)$	${Y=b}_1/\left[1-\left(1-\frac{b_1}{Y_0}\right){\left(\frac{t_0}{t_1}\right)}^{b_3}\right]$
M4	Lundqvist–Korf ${Y=A\mathit{exp}}^{\left(-B{t}^{-\Gamma }\right)}$	$A={\mathit{exp}}^{\left(X\right)}$ $B={(b}_1+b_2)/X$	$X_0=\frac{1}{2}{t_0}^{-b_3}\left[b_1+{t_0}^{b_3}\mathit{Ln}{(Y}_0)\text{±}\sqrt{{{4b_{2}t}_0}^{b_3}+L_0}\right]$ With, $L_0=(-{b_1{-t_0}^{b_3}\mathit{Ln}(Y_0\left)\right)}^2$	${Y=\mathit{exp}}^{{(X}_0)}{\mathit{exp}}^{\left[-\left(\frac{b_{1+}{b}_2}{X_0}\right){t_1}^{{-b}_3}\right]}$
M5	Cieszewski $Y=\frac{B{t}^{\Gamma }}{t^{\Gamma }+A}$	$A=B+X$ $B=\frac{A}{X}$	$X_0=h_0-b_1+\sqrt{{\left(h_0-b_1\right)}^2+2h_0{\mathit{exp}}^{\left(\frac{b_1}{t_0^{b_3}}\right)}}$	${Y=Y}_0\left[\frac{{t_1}^{b_3}({t_0}^{b_3}{X}_0{\mathit{exp}}^{{(b}_2)})}{{t_0}^{b_3}({t_1}^{b_3}{X}_0{\mathit{exp}}^{{(b}_2)})}\right]$
M6	Schumacher ${Y=\mathit{exp}}^{\left(A+B{t}^{-1}\right)}$	$A=X$ $B=b_{1}X$	$X_0=\frac{\mathit{Ln}{Y}_0}{\left(\frac{t_0+b_1}{t_0}\right)}$	${Y=\mathit{exp}}^{\left[X_0+X_0\left(\frac{b_1}{t_1}\right)\right]}$

where Y and Y₀: variables of interest at age t₀ and t₁, respectively; t, t₀, and t₁: stand ages (months); X: unobservable and unquantifiable theoretical variable; X₀: solution of X for initial height and age; A, B, and Γ: parameters of the base model; b₁, b₂, and b₃: global parameters of the dynamic equations. Source: [7,19,23,39].

Table 1. Dominant height growth models selected to fit data from teak (Tectona grandis Linn F.) clonal plantations in the Eastern Brazilian Amazon.

Model No.	Base Model	Parameters Related to the Target Variable X	Initial Solution for X with Y0 and t0	Dynamic Equation
M1	Lundqvist–Korf ${Y=A\mathit{exp}}^{\left(-B{t}^{-\Gamma }\right)}$	$B=X$	$X_0=-Ln\left(\frac{Y_0}{b_1}\right){t_0}^{b_3}$	${Y=b}_1{\mathit{exp}}^{\left[\mathit{Ln}\left(\frac{Y_0}{b_1}\right){\left(\frac{t_0}{t_1}\right)}^{b_3}\right]}$
M2	Chapman–Richards $Y=A[1-{{\mathit{exp}}^{(-Bt)}]}^{\Gamma }$	$A=X$	$X_0=\frac{Y_0}{{[1-{\mathit{exp}}^{\left({-b}_1{t}_0\right)}]}^{b_2}}$	${Y=Y}_0{\left[\frac{1-{\mathit{exp}}^{(-b_1{t}_1)}}{1-{\mathit{exp}}^{(-b_1{t}_0)}}\right]}^{b_2}$
M3	Hossfeld $Y=\frac{A}{1+B{t}^{-\Gamma }}$	$B=X$	$X_0={t_0}^{-b_3}\left(\frac{b_1}{Y_0}-1\right)$	${Y=b}_1/\left[1-\left(1-\frac{b_1}{Y_0}\right){\left(\frac{t_0}{t_1}\right)}^{b_3}\right]$
M4	Lundqvist–Korf ${Y=A\mathit{exp}}^{\left(-B{t}^{-\Gamma }\right)}$	$A={\mathit{exp}}^{\left(X\right)}$ $B={(b}_1+b_2)/X$	$X_0=\frac{1}{2}{t_0}^{-b_3}\left[b_1+{t_0}^{b_3}\mathit{Ln}{(Y}_0)\text{±}\sqrt{{{4b_{2}t}_0}^{b_3}+L_0}\right]$ With, $L_0=(-{b_1{-t_0}^{b_3}\mathit{Ln}(Y_0\left)\right)}^2$	${Y=\mathit{exp}}^{{(X}_0)}{\mathit{exp}}^{\left[-\left(\frac{b_{1+}{b}_2}{X_0}\right){t_1}^{{-b}_3}\right]}$
M5	Cieszewski $Y=\frac{B{t}^{\Gamma }}{t^{\Gamma }+A}$	$A=B+X$ $B=\frac{A}{X}$	$X_0=h_0-b_1+\sqrt{{\left(h_0-b_1\right)}^2+2h_0{\mathit{exp}}^{\left(\frac{b_1}{t_0^{b_3}}\right)}}$	${Y=Y}_0\left[\frac{{t_1}^{b_3}({t_0}^{b_3}{X}_0{\mathit{exp}}^{{(b}_2)})}{{t_0}^{b_3}({t_1}^{b_3}{X}_0{\mathit{exp}}^{{(b}_2)})}\right]$
M6	Schumacher ${Y=\mathit{exp}}^{\left(A+B{t}^{-1}\right)}$	$A=X$ $B=b_{1}X$	$X_0=\frac{\mathit{Ln}{Y}_0}{\left(\frac{t_0+b_1}{t_0}\right)}$	${Y=\mathit{exp}}^{\left[X_0+X_0\left(\frac{b_1}{t_1}\right)\right]}$

where Y and Y₀: variables of interest at age t₀ and t₁, respectively; t, t₀, and t₁: stand ages (months); X: unobservable and unquantifiable theoretical variable; X₀: solution of X for initial height and age; A, B, and Γ: parameters of the base model; b₁, b₂, and b₃: global parameters of the dynamic equations. Source: [7,19,23,39].

2.4. Model Selection and Validation

For each approach, the model that showed the highest accuracy, as well as logical behavior and biological realism in the dominant height growth curves, was selected by evaluating the inflection point, growth rate, and asymptotic point [40]. The model selection adopted the following statistical criteria: highest Pearson’s linear correlation coefficient between the observed and predicted values ($r_{\widehat{y}y}$), lowest root mean square error (RMSE), lowest root mean square error percentage (MSPE), and lowest Akaike information criterion (AIC) [41]. We also performed the graphical analysis of the absolute residues, distribution of the observed and predicted values, histogram of the relative error frequency with residual classes that represent ranges of values in which these residuals are grouped, significance of the regression parameters, and the normality of the residuals by the Shapiro–Wilk test at 95% probability [42].

Complementarily, we performed fittings of the dominant height growth models using the least squares method—generalized nonlinear, from the “gnls” function of the “nlme” package, from R^® studio software, 4.3.1. version [43]. For assessing potential autocorrelation, we modeled the error term using a first-order continuous autoregressive error structure, a method that allowed the models to be applied to longitudinal and unbalanced data [44]. We assumed an autoregressive structure of residuals, whose models were fitted separately [45].

We divided the database of dominant heights into two random groups, considering all ages, with one group for model fitting (80%) and the other one for model validation (20%) (

Table 2. Descriptive statistics of variables of teak (Tectona grandis Linn F.) clonal plantations in the Eastern Brazilian Amazon.

Variables	Fit Data				Validation Data
	Min.	Max.	Mean	SD	Min.	Max.	Mean	SD
t	1.00	12.00	5.31	2.65	1.00	12.00	5.68	2.51
Dh	3.14	21.18	13.84	3.49	5.22	20.31	14.34	2.95
Plots	1.0	46.0	23.5	13.42	1.0	12.0	6.5	3.61
n	1.0	427.0	214	123.41	1.0	107.0	54.0	31.03
V	68.35	227.53	158.84	31.23	-	-	-	-
dbh	17.51	24.44	20.93	2.04	-	-	-	-

where t: age (years); Dh: dominant height (m); n: number of observations (pairs of dominant height and age); V: total volume per plot (m³ ha⁻¹); dbh: diameter at 1.3 m from the ground (cm); SD: standard deviation.

). For validation of the selected models, we adopted the regression-based TOST (two one-sided test) equivalence test using bootstrap, since this is one of the most suitable methods for assessing equivalence between estimated and observed values [46,47]. The dissimilarity hypothesis is rejected or not rejected according to the equivalence regions for the regression parameters (intercept and slope, with 25% at the 99% probability level, with 1000 bootstraps). This test is a statistical procedure used to evaluate whether a model is equivalent to a given standard or reference [12]. We also performed a linear regression between the observed and predicted values to calculate two confidence limits for the parameters and their respective comparison with the estimated equivalence region. Part of the hypothesis of the present article, which is about the accuracy of the models and, consequently, the approach to estimating the dominant height, was tested through statistical criteria for model selection, model validation, and evaluation of biological realism in the estimation of site index curves. This validation test was also applied to clonal teak stands in the studies by [12,26].

2.5. Hierarchical Cluster Analysis

To conclude the hypothesis test of the present article, we used the extent of agreement between the site index (SI) and the volumetric production for each approach method. This procedure was performed through hierarchical clustering analysis, using the “dendextend” package and applying the “tanglegram” function [48]. This test was conducted to assess the extent of agreement between the volume production and site indices of each classification approach. To achieve it, dendrograms were created for ADA and GADA to support comparisons of the similarity among the inventoried plots. We performed the grouping comparison (cluster) from matrices, comparing in pairs, dendrograms consisting of SI data and total volume (m³ ha⁻¹). The total volume with bark consisted of the net production, which was calculated by the sum of the remaining production with the accumulated production in the thinning practices at 7 years of plantation age. The age of 7 years was chosen for comparison purposes, since at that age we had volume monitoring available in all 58 sampling plots of the forest inventory.

We compared each pair of dendrograms (volume × SI method) using the tanglegram method, which checks the existing correlation for each plot individually, comparing its position in each of the two dendrograms to be compared. In the volume dendrogram, three groups (clusters) were expected to be formed, aiming to compare them with the three classes of productive capacity, generated by each method. The distances between clusters were recalculated by using the Lance–Williams dissimilarity update formula, and the unweighted pair group method with arithmetic mean (UPGMA) was selected. The degree of agreement between the two dendrograms was expressed as the entanglement metric (based on the positioning of each plot in the two dendrograms), where zero indicates no entanglement and one is full entanglement.

3. Results

3.1. Modeling of Production Capacity

Dominant height modeling using ADA and GADA revealed the following order of accuracy: M3 > M1 > M4 > M5 > M2 > M6, with estimation errors of less than 1 m (RMSE) and 6% (RMSPE) (

Table 3. Estimators and precision statistics of dominant height growth models fitted using ADA (M1, M2, and M3) and GADA (M4, M5, and M6) approaches in clonal teak (Tectona grandis Linn F.) plantations in the Eastern Brazilian Amazon.

Model No.	Parameters	Standard Error of the Parameters	${\mathit{r}}_{\widehat{\mathit{y}}\mathit{y}}$	RMSE	RMSPE	AIC
M1	b1 = 24.958604	1.2115007	0.9678	0.6683	4.5251	875.63
	b3 = 0.605413	0.0382228
M2	b1 = 0.0232586	0.00201532	0.9596	0.7687	5.2047	995.12
	b2 = 0.7577856	0.03424046
M3	b1 = 20.870433	0.5513357	0.9678	0.6676	4.5199	874.65
	b3 = 1.102182	0.0416282
M4	b1 = −42,437.13	4286.562	0.9677	0.6684	4.5254	877.67
	b2 = 136,544.68	14,206.718
	b3 = 0.6054006	0.035
M5	b1 = 10.707307	3.184992	0.9648	0.7029	4.7595	920.76
	b2 = 5.031924	0.507783
	b3 = 0.751607	0.026144
M6	b1 = −4.44418	0.08000048	0.9589	0.8324	5.6363	1059.93

where b_i: model parameters; RMSE: square root of the mean error; RMSPE: square root of the percentage mean error; $r_{\widehat{y}y}$: correlation coefficient between observed and predicted values; and AIC: Akaike Information Criterion.

). The M3 and M4 models showed the highest accuracy in each approach for predicting the dominant height according to the statistical criteria, with estimation errors of less than 1 m (RMSE), errors between 4 and 5% (RMSPE), highest correlation ($r_{\widehat{y}y}$), and the smallest value of Akaike’s criterion (AIC). The ADA models, M1 and M3, stood out when compared to the GADA models, except for M2, which presented the second highest error (0.7687 and 5.2047%).

The distributions of the residuals for the ADA models showed similarities, especially for models M1 and M3, which revealed less bias in the estimation of the dominant height both in the constructed statistics and in the residue patterns. In both models (Figure 2a), the amplitude of ±2.4 m of absolute residual was concentrated in the central error classes, indicating a greater number of observations close to zero, which was confirmed by the Shapiro–Wilk test, concluding that the residuals follow a normal distribution (Figure 2c), with significant and high correlations between observed and predicted values (>0.96) (Figure 2b). When evaluating model M2, a particular trend was noted, with residue amplitudes of ±2.8 m, evidencing underestimation and overestimation in the smaller and larger dominant heights, respectively.

On the other hand, model M6 of the GADA approach presented underestimations of up to 2.91 m in the prediction of the dominant height at early ages, while it showed a trend of overestimation of up to 2.8 m (Figure 3a). Model M4, in turn, showed no bias in the estimations, presenting residual amplitudes of ± 2.3 m (Figure 3a). Additionally, this model showed a high correlation between the estimated and observed values (>0.96) (Figure 3b), and the residuals of this model followed a normal distribution, concentrating on the central error classes (Figure 3c).

The ADA and GADA models indicated accurate estimates for dominant height and were validated by the equivalence test, except for model M6, in which the intercept parameter of test (14.44 ± 14.81), corresponding to the similarity region, extrapolated the acceptance margin (14.14 ± 14.64) (

Table 4. Validation of the ADA and GADA models for estimating the dominant height in teak (Tectona grandis Linn F.) clonal plantations in the Eastern Brazilian Amazon.

Model No.	Parameters	Confidence Interval	Similarity Region	Dissimilarity
M1	Intercept	14.34 ± 14.84	14.48 ± 14.78	Reject
	slope	0.75 ± 1.25	0.95 ± 1.07	Reject
M2	Intercept	14.33 ± 14.83	14.46 ± 14.79	Reject
	slope	0.75 ± 1.25	0.90 ± 1.04	Reject
M3	Intercept	14.35 ± 14.85	14.47 ± 14.78	Reject
	slope	0.75 ± 1.25	0.96 ± 1.07	Reject
M4	Intercept	14.34 ± 14.84	14.47 ± 14.79	Reject
	slope	0.75 ± 1.25	0.96 ± 1.08	Reject
M5	Intercept	14.31 ± 14.81	14.47 ± 14.79	Reject
	slope	0.75 ± 1.25	0.95 ± 1.08	Reject
M6	Intercept	14.14 ± 14.64	14.44 ± 14.81	Do not reject
	slope	0.75 ± 1.25	0.89 ± 1.06	Reject

). For model M1, the test generated a reliability interval for the intercept of 14.34 ± 14.84, and the similarity region (14.48 ± 14.78) was within this interval. The same pattern was observed for the slope parameter, with the region of similarity (0.95 ± 1.07) contained within the confidence interval (0.75 ± 1.25). Therefore, the test of the straight-line equation parameters estimated values within the confidence region, indicating reliability of the results. The validation confirmed the accuracy of models M1, M2, M3, M4, and M5, to re-estimate the dominant height, indicating no statistical difference between observed and predicted values. The hypothesis of dissimilarity was rejected for the models fitted by ADA and GADA, since the confidence intervals of their parameters were found within the region of similarity.

When analyzing the density of the dominant height estimation performed by the approaches, we found a greater overlap between observed and predicted values for models M1 and M4, followed by M3 (Figure 4). On the other hand, models M2, M5, and M6 presented larger oscillations between the estimated and observed areas. These results corroborate the greater errors obtained for these models in estimating the dominant height, as also evidenced by the equivalence test, which did not validate model M6.

For model M6, larger discrepancies were obtained between the observed and estimated values, especially in the range from 3 to 15 m of dominant height. The lack of agreement between these heights indicates a lower accuracy of model M6 with regard to the observed data. This probability density analysis between observed and estimated values from modeling provided a more detailed understanding of the performance of each model with respect to the estimation of the dominant height. Models M1 and M4 showed better adherence to the observed data, while M2, M5, and M6 showed greater discrepancies and oscillations, indicating their lower accuracy.

The site index curves obtained in each approach highlight the existence of sites with distinct productive characteristics, reflecting the variation in teak productive potential. This variation was represented by site indexes at 16, 18, and 20 m (Figure 5), indicating different levels of tree growth quality. Both approaches yielded equivalent results, revealing that most plots were concentrated in the intermediate yield capacity classes, corresponding to 63% of the total plots. On the other hand, the lower and upper classes showed smaller proportions, representing 23% and 13% of the plots, respectively. This distribution indicates that the study area has a significant number of plots with moderate production capacity, while plots with lower or higher productive potential are less common.

Although the quality of fitting values indicated that model M3 had greater accuracy in the ADA approach, the site index curves generated by this model presented unrealistic shapes, as illustrated in Figure 5—M3. These inconsistencies were observed when projecting curves that exceeded the real limits of the dominant height values, especially in the most productive class. This projection compromised the adequate representation of the productive characteristics of the sites. On the other hand, model M2 generated anamorphic curves (Figure 5—M2), which were not able to include the dominant height values among the classes, resulting in a distorted representation of the forest sites’ characteristics. The GADA models, M5 and M6, in turn, projected curves that could not include the dominant height values at ages lower than 7 years, not adhering to the actual observed values. This limitation compromised the accuracy of the representation of forest site characteristics at younger ages.

The best models, M1 from ADA (Figure 5—M1) and M4 from GADA (Figure 5—M4), generate polymorphic site index curves that consistently follow the biological behavior of the dominant height. The actual points of the dominant height variable were plotted broadly within the boundaries of the protected site classes. These curves provide an accurate representation of the productive characteristics of the forest sites. Both the ADA and GADA approaches provided equivalent heights for the estimation of dominant height and the construction of site curves. However, the ADA approach stood out by providing a more accurate classification of the forest sites. This means that the ADA approach was able to capture site characteristics more accurately, generating site index curves that are more faithful to reality.

We provide a practical example of how to use the equation developed in this study, specifically focusing on the algebraic difference approach represented by model M1. Equation (2) enables the projection of the future dominant height (Dh₁) based on the current dominant height (Dh₀), current age (t₀), and future age (t₁). When the future age (t₁) is equal to the reference age (RA), which is 144 months (or 12 years) in this case, the future dominant height (Dh₁) equals the site index (SI). This process culminates in the construction of Equation (3), which can be employed to estimate the dominant height (Dh). Furthermore, by inverting Equation (3), we obtain Equation (4), which allows for the prediction of the site index based on the dominant height (Dh), demonstrating consistency of the methodology employed.

$${Dh}_1=b_1·{\left(\frac{{Dh}_0}{b_1}\right)}^{{\left(\frac{t_0}{t_1}\right)}^{b_3}}=24.958604·{\left(\frac{{Dh}_0}{24.958604}\right)}^{{\left(\frac{t_0}{t_1}\right)}^{0.605413}}$$

(2)

$$Dh=b_1·{\left(\frac{SI}{b_1}\right)}^{{\left(\frac{144}{t}\right)}^{b_3}}=24.958604·{\left(\frac{SI}{24.958604}\right)}^{{\left(\frac{144}{t}\right)}^{0.605413}}$$

(3)

$$SI=b_1·{\left(\frac{Dh}{b_1}\right)}^{{\left(\frac{t}{144}\right)}^{{\mathrm{b}}_3}}=24.958604·{\left(\frac{Dh}{24.958604}\right)}^{{\left(\frac{t}{144}\right)}^{0.605413}}$$

(4)

3.2. Relationship between Site Index and Volume Production

The analyses of the ADA and GADA approaches, despite revealing differences in the estimation of dominant height by models M1 and M4, indicated the same results in the dendrograms. These approaches generated the same levels of agreement (both for Volume × IS-ADA and Volume × IS-GADA) when comparing the respective site indices with the volumetric yields for each sample plot (Figure 6). The resulting dendrogram was stratified into three distinct volumetric yield groups. Group 1 (represented in green) comprised volume values between 189.44 and 227.53 m³ ha⁻¹, Group 2 (represented in orange) comprised values between 131.98 and 188.07 m³ ha⁻¹, and Group 3 (represented in blue) encompassed values between 68.35 and 118.99 m³ ha⁻¹. This volumetric stratification allowed for a more congruent comparison with the three groups represented by the 16, 18, and 20 m site indexes.

We observed that the ADA and GADA approaches represented by models M1 and M4, respectively, besides faithfully representing the productive capacity, presented a strong cause-and-effect relationship with the volume, allowing us to observe a greater in-fluence of the productive capacity on the volumetric production of teak in a more explicit way. The approaches with dynamic equations generated equal site indexes for each plot and, consequently, the same level of agreement (mismatch) between volume and site index (0.1764). As ADA presents greater simplicity of model application, slight superiority in accuracy and the same level of agreement when compared to GADA, we recommend the ADA approach for classification of production units as indicative of the volumetric production of the evaluated teak clonal stands.

By analyzing Figure 7, it is possible to observe the distribution of site classes regarding volumetric production, with classes presented in descending order. We found that the volumetric production hierarchy agrees with the site classes, which indicate higher, intermediate, and lower productive capacities, according to the ADA and GADA approaches. This means that plots with higher volume production are associated with site classes that have a higher productive capacity, while plots with lower volume production correspond to site classes with lower productive capacities. The analysis revealed a 96.5% accuracy rate in the volume ranking order with regard to site classes for both approaches. There were discrepancies only in two pair cells, whose volumes were 143.58 and 155.72 m³ ha⁻¹, with regard to the predicted site classes. These results highlight the relationship between volume and site index, demonstrating that volumetric production has a strong responsiveness with the productive capacity of the different site classes.

4. Discussion

The constructed models from the ADA and GADA approaches provided accuracy in estimating the dominant height, especially for models M1 and M4, which are Lundqvist–Korf base models. These results may be related to the fact that ADA and GADA are dynamic auto reference equations, which express the variable Dh at future age as a function of future age and Dh at initial age as a function of initial age [7,19,26], considering a base age of 12 years for our study. Conventional biological models of dominant height growth are widely used to classify teak production units with anamorphic curves [49]. However, errors above 1 m in height are common, and trends in the distribution of residuals can be observed in estimates [11,49]. Several studies have achieved efficient results for modeling the productive capacity for teak stands using dynamic equations, both for ADA, e.g., [12] and GADA, e.g., [7,26,50], indicating the predictive potential of these approaches for estimating the dominant height of teak stands.

The selection of mathematical models should also consider the biological interpretation of their estimates, and not only conventional accuracy by statistical criteria [19] and validation by comparison tests. Although model M4 proved to be more accurate, based on statistical criteria, it was not enough for its implementation in the construction of the site index curves, since this model showed inconsistency related to the growth trend of teak, which was not able to represent its biological behavior realistically. It is highlighted that statistical accuracy does not always represent an assurance of biological consistency, especially when it comes to the growth dynamics of planted forest stands [17].

Regarding the analysis of the residuals’ dispersion, we observed an unstable distribution in the early ages of the plantings for model M6, and this model was the only one not validated. The plots of the generated residuals, from each approach, showed lower trends for the M1–ADA and M4–GADA dynamic models. This behavior shows an association with the nature of the data used in the dynamic equations, as they consider simultaneous comparisons of the change rate of dominant height with age [51]. The characteristic of differential equations enables them to obtain the evolution pattern of the variable and, consequently, to generate more accurate estimates that are consistent with the biological behavior of the species, as previously described. Age-invariant equations, such as ADA and GADA, polymorphic and with one or multiple asymptotes, allow for the derivation of dynamics when appropriate, which gives more realistic predictions of the variable of interest [19].

Among the models assessed by the present study, model M2, which generated anamorphic curves, was not able to describe the behavior of the real values of Dh. According to [52], this can be justified by the fact that the proportionality behavior (anamorphic) is not able to describe the biological realism of the dominant height as a function of stand age in a reliable way, even when we use the algebraic approach (ADA). Those authors also claim that the curve shapes vary between sites and are dynamic over time, and for this reason, the dominant height structure has a polymorphic character. According to [53], the height growth curve assumes a more accentuated sigmoid shape in places with higher productivity, indicating that the proportionality between the local curves is flawed. On the other hand, in less productive sites, the height growth pattern is smoother, with the inflection point being reached later than in more productive locations. Therefore, the use of polymorphic models results in growth curves which are more representative of biological reality.

The meristematic activity of the studied species is responsive to the environment, resulting in different height growth rates for different forest sites, which justifies the polymorphic character of the dominant height growth curves [12]. The polymorphic curves, generated by ADA and GADA, reflected the cloud scattering trend of the observed dominant height values, suggesting that the algebraic difference method is efficient and accurate for determining teak stand site indices. The site indices generated by model M1 (ADA) provided site class centers ranging from 16 to 20 m, and it was the most accurate model among the approaches.

The existence of different sites, represented by the SI of 16, 18, and 20 m, was verified, representing distinct bioedaphoclimatic conditions. One of the justifications for the differences in productive potential observed for these stands lies in the fact that there are differences in the concentration of cationic macronutrients in the soil, such as Ca and Mg, which are important for cell structure and formation of bark and wood and composition of chlorophyll molecules, respectively [54,55,56]. This influence of soil exchangeable Ca and Mg content on the dominant height of teak stands was found by [57], in which higher content of these nutrients indicated sites with higher productive capacity in a significant way.

To predict the dominant height growth and, consequently, design the site index curves of the teak clone plantations studied, the ADA (M1) modeling framework was superior to GADA (M4) when assessing accuracy. However, modeling using GADA (M4) was also shown to be adequate; therefore, this model can be recommended, as also observed in the study of [26], which reported an error of 0.67 m when using the same Lundqvist–Korf base model, also called the Bailey–Clutter model. This approach is suitable to explain the growth rate in such a way that more than one parameter of the model has dependence on the production capacity of the site, making the curves more flexible, and thus able to portray polymorphism with multiple asymptotes [24].

With the analysis of the tanglegram, in which the dendrograms of the total volume and site index for each approach method were compared, we observed that the classification approach of the productive units of the investigated plantation can serve as an indicator of volumetric production, as well as being the response of the volume in relation to the estimated productive capacity. This result demonstrates that approaches that encompass the growth dynamics of dominant height are more accurate and efficient in indicating productivity, ratifying the assumption that dominant height is an effective integrator in indicating and responding to biological determinants of growth [44].

When sites are properly classified, there is a smaller tendency for systematic errors to influence the prognosis of stands’ growth and production, in which the site index is a significant factor. The expected biological behavior in areas with greater productive capacity is that the trees will present greater volume, and the opposite occurs in areas with lower productive capacity, as shown in Figure 7. In this regard, when evaluating the accuracy, validation, and biological realism of the approaches and their correspondence with the volumetric production, we cannot reject the formulated hypothesis that the ADA and GADA approaches offer an accurate and realistic representation of the volumetric production, through the site index. Therefore, it is inferred that these dynamic equations and their polymorphic curves can serve as indicators of the volumetric production of the teak clonal stands that were evaluated.

5. Conclusions

This study aimed to determine the most accurate approach for assessing the productive capacity of clonal teak plantations and its relationship with volumetric production. Among the six models evaluated, the Lundqvist–Korf equations for ADA and GADA yielded the highest accuracy, with errors lower than 0.67 m and 5.53%, consistently describing the evolution of dominant height. We observed that dynamic equations generating polymorphic curves effectively represent the sites, providing a connection and indication of the plantations volumetric production, with an agreement rate of 98.3%. Based on our findings, we recommend using the ADA and GADA approaches for estimating the dominant height of clonal teak plantations in the Eastern Brazilian Amazon.