Sergio Orrego · Cristian Montes · Héctor I. Restrepo ·Bronson P. Bullock · Mauricio Zapata

Abstract Teak ( Tectona grandis L.f.) plantations are increasingly being established in tropical regions to meet a rising demand for its highly valued timber. Teak plantations have been established in the Atlantic Coastal Plain region of Colombia, a region climatically suitable for teak growth by having a monsoon climate with a unimodal precipitation pattern. Tree diameter at breast height (DBH, 1.3 m above ground) and mean top height, periodically measured over a 17-year period in 44 permanent sampling plots of size 0.06 and 0.10 ha, were used in this study. A stochastic diff erential equation (SDE), along with a Bertalanff y-Richards-type height growth model, was used to model and estimate top height growth of teak plantations in Colombia.Environmental noise and height measurement errors were explicitly considered as the main uncertainty sources of mean top height growth. The best model for estimating mean top height, based on statistical performance and biological rationale, had the asymptote def ined as a local parameter and the growth rate and shape specif ied as global parameters. This model outperformed its counterpart that had the growth rate specif ied as a local parameter and asymptote and shape as global parameters. The selected model also outperformed alternative approaches such as the mixed-Effects model, generalized algebraic diff erence approach, and the dummy variable method. Estimated trajectories for the mean top height of teak in Colombia are biologically sound based on the measured height series and previous studies in Latin America. Results suggest that most of the uncertainty associated with the mean top height growth of teak plantations in Colombia was largely explained by environmental noise. The best estimated model using the SDE approach can be useful for predicting height growth and evaluating site productivity of teak plantations in Colombia and in neighbouring countries with biophysical characteristics similar to those where teak has been planted in Colombia.

Keywords Mean top height · Stochastic diff erential equation · Forest productivity · Timber production ·Timberland investment

Introduction

Teak (Tectona grandisL.f.) is one of the most valuable tropical hardwoods for its natural durability and broad range of uses (e.g., outdoor furniture, boat decks, and indoor f looring). It is naturally distributed in Myanmar, India, Laos, and Thailand (Tewari et al 2014). However, natural teak forests are steadily declining in area, leading to conservation eff orts and the imposition of logging bans legally enacted almost three decades ago (Kollert and Kleine 2017). The logging bans and increasing demand for tropical hardwood (Kollert and Kleine 2017) are the main drivers of the observed increase of teak plantations.

The worldwide estimated area of teak plantations is 4.35-6.89 million ha (Kollert and Kleine 2017), with the greatest area (ca. 80%) in Asia. In tropical America, the f irst teak plantations were established in Trinidad in 1923, then in Honduras, Panama, Costa Rica (Tewari et al 2014), Brazil(Coutinho 2004), and Ecuador. Teak plantations have been present in the Atlantic Coastal Plain region of Colombia for 70 years, with f inancial analysis showing a high prof itability for this type of forest investment (Restrepo and Orrego 2015).

Height modeling is extensively used in forestry to evaluate site productivity of even-aged stands of a single forest species (Tesch 1981; Skovsgaard and Vanclay 2008; Burkhart and Tomé 2012), and height is considered a fundamental structural component of a forest growth model(García 1983; Fang et al 2001). To estimate teak growth,traditional approaches (Malende and Temu 1990; Bermejo et al 2004; Upadhyay et al 2005) and diff erential equations and advanced modeling techniques have been used, particularly in India (Tewari et al 2014; Tewari and Singh 2018).However, to the best of our knowledge, advanced modeling techniques and SDEs have not been used in Latin America to evaluate height growth and site productivity of teak plantations.

Height growth modeling of teak in tropical America has predominantly used traditional methods such as guide curve or algebraic diff erence approach, least squares or nonlinear mixed-Effects estimation methods, and Schumacher-,Korf-, Hossfeld- or Bertalanff y-Richards-type equations(Keogh 1982; Bermejo et al 2004; Torres et al 2012). For a provisional teak site classif ication, least squares methods were used to estimate a Schumacher-type height model with 144 observations of top height, dominant height and dominant-codominant height from 13 Latin American and Caribbean countries (Keogh 1982). In that study, the guide curve approach was used to estimate and identify 11 site index curves, with 25 years as the base age; results indicated that height growth rates were lower than at sites in India and Myanmar where teak grows naturally. To evaluate the observed underestimation of height growth rate using the Schumacher model, Keogh ( 1990) used more observations and a log-log transformation, the results were similar, leading to the recommendation of the use of repeated measurements.

Yearly inventory data and the guide curve method were used to estimate Bertalanff y-Richards-, Schumacher-, and Hossfeld-type height models for intensively managed teak plantations in Costa Rica (Bermejo et al 2004). The Hossfeld-type equation was considered the best and used to draw height curves for three site indices (23, 21, and 19 m at a base age of 10 years) representing best, average, and worst site qualities, respectively (Bermejo et al 2004). Torres et al( 2012) used Bertalanff y-Richards- and Korf-type equations to model teak height growth with repeated measurements from 44 PSPs in the Atlantic Coastal Plain region of Colombia. Nonlinear mixed-Effects models specif ied with a single random Effect and alternative variance structures (constant,power and exponential function of the mean) were f itted, and site quality classes were identif ied based on deviations from the average site. The Korf-type equation was the best model based on statistical performance. Five site index classes,using a base age of 12 years, were depicted, with site index of about 23 m for the highest quality site (Torres et al 2012).

The productivity of sites where teak is being planted in Colombia can be evaluated using an alternative to the traditional algebraic diff erence approach (Bailey and Clutter 1974; Clutter et al 1992), and its generalization (Cieszewski and Bailey 2000). In the present study, the SDE approach(García 1983; Seber and Wild 1989; Rennolls 1995) and the Bertalanff y-Richards (von Bertalanff y 1957; Richards

1959) functional form were used to model height growth for teak plantations in Colombia. The standard nonlinear least squares method, as recently proposed by García ( 2019), was utilized to estimate SDE parameters of the Bertalanff y-Richards-type model. The estimation was made using repeated measurements from 44 PSPs, the same data set used by Torres et al ( 2012), one of the few longitudinal data sets of teak height growth in Latin America.

The SDE modeling approach used in this study considers measurement error and a stochastic term. While the former accounts for the error associated with height measurements(Goelz and Burk 1996), the latter accounts for the uncertainty associated with height trajectories over time as a result of environmental-related noise. The explicit treatment of these two sources of uncertainty in SDE is not possible in the most common approaches to estimate height growth models such as the guide curve method (Clutter et al 1992),algebraic diff erence approach (Bailey and Clutter 1974;Cao 1993), and generalized algebraic diff erence approach(Cieszewski and Bailey 2000; Cieszewski 2001). Moreover,SDE allows for modeling growth as a system (e.g., a teak stand) considering changes in growth dynamics as a result of silvicultural interventions or natural process (Nord-Larsen 2006).

A general description of the study area is presented in the next section, along with the main features of the data, the methods used to calculate the mean top height, and the specif ication of the height growth as a SDE. The estimation of two alternative height growth models for teak plantations in Colombia are presented in the result section. In the models,either the asymptote or growth rate was def ined as a local parameter. The remaining parameters were def ined as global.A robustness analysis was conducted by comparing SDE results with the generalized algebraic diff erence approach,mixed-Effects models, and dummy variable method as alternative modeling techniques. We also discussed our results in the light of previous studies of teak height growth. Finally,some general concluding remarks are presented.

Materials and methods

Study area

The f irst teak plantations in Colombia were established in the Atlantic Coastal Plain region 70 years ago. The monsoon climate, similar to that where teak grows naturally in Asia, land availability, and private forest investments were the main drivers of the teak plantations in Colombia. The region has a typical unimodal precipitation pattern throughout the year, with a mean annual precipitation of 2500 mm and mean annual temperature of 27 °C. The landscape is predominantly plain with alluvial fertile soils. Most of these plantations were established on grasslands as an alternative to the livestock economic activity.

Data and mean top height calculation

Similar management practices were implemented in all the plots: planting density of 1600 seedlings/ha, weeding control by hand, pruning at years 5 and 9, and two thinnings at years 7-9 and 12-13 to maintain basal area at 26 m 2 /ha (Torres et al 2012). This study uses data from 44 PSPs of size 0.06 and 0.10 ha, established between 1982 and 2000, in which the DBH and height (H, 10% of the trees) of the trees were measured between 3 and 16 times.

All H-DBH measurement pairs were used to f it four H-DBH models by plot and inventory occasion, representing distinct functional forms (general Schumacher, Schumacher subtracting a breast height of 1.3 m, power, and polynomial),and commonly used in forestry for prediction of tree height:

The selected H-DBH model was used in the imputation of unmeasured tree heights. There exists empirical evidence about the importance of the def inition of top height and its Effect on the estimation of site index (Sharma et al 2002).The def inition of mean top height of the 100 largest trees per hectare is consistent with the f ield protocol followed for forest inventories of teak plantations in Colombia. It implies the use of 6 and 10 largest trees in plots of size 0.06 and 0.10 ha,respectively. The plot size and the spatial distribution of the largest trees on the ground, with largest trees potentially grouped in clusters, may induce a bias selection Effect in the estimation of the mean top height (Rennolls 1978; García and Batho 2005). This potential bias was overcome using the 9 and 15 largest trees for each plot and for each measurement, corresponding to those integers obtained after using the correction factor 1.6A− 0.6 as suggested by García and Batho ( 2005), whereAcorresponds to 6 and 10, respectively.

Height growth equation

Consider the model of height growth as a function of height(García 2019):

wheredH∕dtis the height growth,β0is the growth rate or time-scale,β1is the asymptote, andβ2is a power transformation which allows for f lexibility of the functional form. Equation 5 is equivalent to the well-known Bertalanff y-Richards equation (von Bertalanff y 1957; Richards 1959). The integration of Eq. 5 will lead to the representation of the height as a function of age and constants of the initial height (h0)and age (t0), as follows:

Forestry applications typically assume thatH(0)=0 ,but any other initial values might be used. Height growth is substantially inf luenced by environmental conditions (e.g.,climate, soil, and topography). Therefore, environmental anomalies or perturbations will lead to deviations of the height growth trajectory. These deviations are more likely to be relevant as a result of the climate variability and current human-induced climate change. If the environmental-related noise, also known as process error, is considered as an additive stochastic term, the Eq. 5 can be written as:

The f irst term of the right side of Eq. 7represents the expected trend of the process, and the second term,σp, corresponds to the uncertainty or volatility.Note that the uncertainty term does not depend onH.Moreover, heights for most of the trees are not measured but estimated, and even those measured are likely prone to error (Goelz and Burk 1996). Therefore, the measurement error is also considered (García 2019) as follows:

where theεijis assumed to be identically and independent normally distributed with mean 0 and varianceThe double subscript denotes the observationi=1,…,njin plotj=1,…,k. An appropriate parameterization of Eq. 5 using a Box-Cox transformation gives:

Parameters are estimated by using the strategy suggested by Furnival and Box-Cox (Box and Cox 1964), and formally developed by García ( 2019) and used in forestry.The standard function nls available in the software R (R Core Team 2020), was used to estimate the parameters of the SDE (Eq. 9). Models with asymptote (β1) and growth rate (β0) specif ied as local parameters (i.e., one parameter for each plot) were estimated, and the logarithm of the maximum likelihood function (log-lik) and the Akaike and Bayesian information criteria (AIC, BIC) were calculated.The estimated parameters of the estimated models were used to draw mean top height trajectories for each plot.

Robustness of the estimated SDE

The goodness of f it of the SDE was compared with that obtained by using the generalized algebraic difference approach (Cieszewski and Bailey 2000; Cieszewski 2001),mixed-Effects models (Torres et al 2012), and the dummy variable method (Wang et al, 2008).

Results

H-DBH equation and the mean top height

Fig. 1 Mean top height as a function of age for teak plantations in Colombia

The H-DBH model corresponding to Eq. ( 3) was selected as the best model, based on a comprehensive analysis of statistical performance, including verif ication of constant variance and normality assumptions. The use of a H-DBH equation at the stand level may result in lower estimated standard deviations in mean top height (Mason 2019), because of the lower height variability at the stand level compared to the plot level. However, the reduction in estimated standard deviations are lower compared to alternative methods, as that in which the average of the diameters of 100 largest trees per hectare are used to determine mean top diameter,and then calculate mean top height (Mason 2019). Both the initial value of the mean top height and its time trend suggest a high heterogeneity among plots (Fig. 1). Based on the initial values in Fig. 1, the initial height growth of teak is approximately in the range 3-6 m/yr. In general, the height growth rate of teak tends to be relatively high a decade later after the establishment of the plantation (Fig. 1 ).

SDE height growth models

Parameter estimates of the SDE, with either asymptote or growth rate specif ied as a local parameter, are shown in Tables 1 and 2, respectively. For the model in which the asymptote was assumed to be local, the estimated parameter ranged from 18.28 to 28.81 m, and the estimated global parameter (i.e., growth rate) was 0.106 (Table 1). The estimated process (̂σp= 0.0027) and measurement (̂σm=0 )errors suggest that all the uncertainty of the height estimation is explained by environmental noise.

The estimated model meets the assumptions of constant variance and approximately normally distributed errors(Fig. 2). The estimated mean top height trajectories are shown in Fig. 3 a. Finally, the results suggest an acceptable relationship between the estimated versus the observed mean top height (Fig. 3 b), with most of the points on the 45° straight line.

For the model in which growth rate was assumed to be local, the estimated parameter ranged from 0.075 to 0.358, and the estimated global parameter (i.e., asymptote) was 21.161 m (Table 2). The estimated measurement error (̂σm= 0.0418) exceeded the environmental noise(̂σp= 0.0040).

A plot of residuals versus estimated mean top height(Fig. 4 a) indicates that errors exhibit constant variance.However, the assumption of normal distribution for errors may not be reasonable (Fig. 4 b). The estimated mean top height trajectories for the model with growth rate as a local parameter are depicted in Fig. 5 a. The estimated asymptote is lower than many of the observed maximum values of height, suggesting that this model lacks biological interpretation. Finally, the relationship between estimated andobserved mean top height is shown in Fig. 5 b, indicating that most of the points deviate of the 45° straight line for the whole range of height values.

Table 1 Parameter estimates and goodness of f it of the stochastic differential equation model, Bertalanff y-Richards type, of the mean top height growth for teak plantations in Colombia. The asymptote was specif ied as a local parameter

The goodness of f it indicates that the SDE model with the asymptote as a local parameter (log-lik = − 331.901,AIC = 759.802, BIC = 955.967) was superior than the SDE model with the growth rate as a local parameter(log-lik = − 705.027, AIC = 1506.054, BIC = 1702.219)(Tables 1 and 2). The selected model captures the environmental variability of site productivity much better than the model with the growth rate as a local parameter (Fig. 3 a vs Fig. 5 a), which has a f ixed asymptote across all plots. In the selected model, the value of the shape parameter (Table 1)indicates that mean top height growth can be characterized by a curve where the inf lection point is attained at an age before the youngest represented in Fig. 3 a.

Table 2 Parameter estimates and goodness of f it of the stochastic differential equation model, Bertalanff y-Richards type, of the mean top height growth for teak plantations in Colombia. The growth rate was specif ied as a local parameter

Fig. 2 Diagnostics of the stochastic diff erential equation model, Bertalanff y-Richards type, of the mean top height growth for teak plantations in Colombia with local asymptote: a studentized residuals versus estimated mean top height (m), b quantile-quantile (QQ) plot

Robustness analysis

To compare our results obtained using the SDE approach,we estimated a nonlinear mixed-Effects model (Pinheiro and Bates 2000) using the package nmle of the software R (R Core Team 2020), considering the asymptote as a random Effect (Table 3). Although some results of the SDE and nonlinear mixed-Effects approaches are not directly comparable (e.g., log-lik), the AIC of the selected SDE model(AIC = 759.802) is lower than the mixed-Effects alternative(AIC = 812.682). Predictions for the random Effect in the mixed-Effects approach (Table 3) were systematically lower than the corresponding local parameters estimated using the SDE approach (Table 1).

Also, a model for mean top height using the generalized algebraic diff erence approach was estimated by considering asymptote and shape as local parameters (Table 4). The model was f it using the base-age invariant method as proposed by Bailey and Clutter ( 1974), estimated through the dummy variable approach. The dummy variable approach does not rely on increments; instead, it uses the integrated form and estimates one local parameter for each plot and global parameters for each equation. Unlike other authors(e.g., Borders et al 1984; Wang et al 2007), we did f it the integral form as presented by Cieszewski ( 2001). This method ensured an error-free independent variable (time)and avoided the problem encountered in the diff erence equation method.

Fig. 3 a Estimated mean top height trajectories and b estimated versus observed mean top height for teak plantations in Colombia using a stochastic diff erential equation model, Bertalanff y-Richards type, with the asymptote as a local parameter

Fig. 4 Diagnostics of the stochastic diff erential equation model, Bertalanff y-Richards type, of the mean top height growth for teak plantations in Colombia with local growth rate: a studentized residuals versus estimated mean top height (m), b quantile-quantile (QQ) plot

The AIC of the selected SDE model (AIC = 759.802)was lower than the generalized algebraic difference approach (AIC = 1110.907). Lower predictions for the asymptote in the generalized algebraic diff erence approach(Table 4) were systematically obtained compared to the results using the SDE approach (Table 1). Moreover, the selected SDE model was better than alternatives for mean top height using the dummy variable approach (results are shown in the supplementary material, Tables S1, S2,Figs. S1-S4).

Discussion

The SDE model, the Bertalanff y-Richards type, with the asymptote as a local parameter was the selected model for the mean top height growth for teak plantations in Colombia. The model is biologically and statistically better than its SDE counterpart with the growth rate as a local parameter,and specif ications such as the mixed-Effects model, generalized algebraic diff erence approach, and the dummy variable approach. The selected model suggests that the uncertainty associated with the height growth is mainly explained by environmental noise, whereas the measurement error is negligible, f inding which is consistent with the literature(García 1983).

Fig. 5 a Estimated mean top height trajectories and b estimated versus observed mean top height for teak plantations in Colombia using a stochastic diff erential equation model, Bertalanff y-Richards type, with growth rate as a local parameter

The selected model generated anamorphic trajectories of mean top height, which describe adequately the observed variability of height growth (Fig. 1) of those sites where teak has been planted in the Atlantic Coastal Plain region of Colombia, consistent with anamorphic site index curves obtained with a Korf type model (Torres et al 2012). However Torres et al ( 2012), used a nonlinear mixed-Effects model to study site index and site productivity of teak, which relies theoretically on the statistical assumption that plots were randomly selected from a population. This assumption could not be valid in top height modeling studies if sampling units were actually chosen by considering predominantly geographic coverage and not strictly a random selection. The selected model in this study led to anamorphic mean top height trajectories (Fig. 3 a), which indicate that the asymptote of the mean top height of teak is local. Local asymptotes seem to explain the variation between sites and some microsite characteristics, mainly nutrients, aff ecting the mean top height growth of teak in Colombia.

The best site seems to indicate a slowing of mean top height growth of teak plantations at year 20. A similar f inding was reported in the teak provisional site classif ication chart for the Caribbean, Central America, Venezuela and Colombia (Keogh 1982), and in a study of height growth of teak in Trinidad and Tobago made in 1957. However, teak height may continue steadily growing after year 20 in India(Keogh 1982).

This study constitutes an advance in the height growth modeling of teak in Colombia by using repeated measurements and a modern modeling technique as is the case of the SDE approach. Site productivity can be evaluated in a more reliable way by using the results from this study, which constitutes an empirical and substantial improvement compared to the provisional chart of site productivity for teak (Keogh 1982, 1990), and to the more recent study of height growth using mixed-Effects models, with the restrictive assumption of randomly selection of sampling units or plots. However,the slowing of mean top height was still observed. It is possible to speculate that more than the modeling strategy, the slowing of the mean top height trajectory in teak can be explained by a soil def iciency. If so, height growth of teak in Colombia may be heightened with an appropriate fertilization schedule.

Table 3 Parameter estimates and goodness of f it of the mixed-Effects model of the mean top height for teak plantations plots in Colombia. The asymptote parameter was specif ied as a random Effect

Conclusions

The SDE approach was proposed several decades ago for being totally consistent with the notion of a system (e.g.,a teak stand) evolving over time. The SDE allows for the specification of height growth as a function of height and environmental-related noise or perturbations determining height dynamics (García 1983, 1999). Because of intrinsic structural complexities associated with the SDEs and the methods for their analytical solution, SDEs did not become widely used in forestry. However, some of these issues have been overcome (García 2019). Here, the SDE approach and the Bertalanffy-Richards functional form were used to model height growth for teak plantations in Colombia, using repeated measurements from 44 PSPs,in which DBH outside bark for all trees, and H for a small number of trees (10%), were periodically measured in each plot. The specification of the asymptote as a local parameter resulted in the best model based on biological rationale and statistical performance, evaluated by comparing log-lik, AIC, and BIC values. The selected model outperformed its counterparts mixed-effects model, generalized algebraic difference approach, and the dummy variable method. In the selected model, the local asymptote is consistent with the variation of site productivity observed in the Atlantic Coastal Plain region of Colombia, and the uncertainty is associated entirely with the environmental noise and not with any measurement error.Site productivity can be evaluated in a more reliable way by using the results provided in this study, which constitutes an alternative approach compared to the previous evaluated studies. Moreover, results may be valuable for the establishment of new plantations in Colombia and for countries in Central America and the Caribbean. The selected model would also be valuable to evaluate the effect of climate variability and climate change on teak growth in the region.

Table 4 Parameter estimates and goodness of f it of the generalized algebraic diff erence approach of the mean top height growth for teak plantations plots in Colombia. The asymptote and shape parameters were specif ied as locals

Acknowledgements This research was done while the lead author was on sabbatic as a visiting researcher in the Warnell School of Forestry and Natural Resources, and the Plantations Management Research Cooperative (PMRC) at the University of Georgia, Athens, GA, USA.The lead author is grateful to the Director and Co-Director of PMRC and to graduate students for providing an exceptional academic and supportive environment. We thank the two anonymous reviewers for their suggestions and comments.

Publisher’s NoteSpringer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affi liations.