Developing a dynamic growth model for teak plantations in India

Tectona grandis (teak) is one of the most important tropical timber speciesoccurring naturally in India. Appropriate growth models, based on advanced modeling techniques,are not available but are necessary for the successful management of teak stands in the country.Long-term forest planning requires mathematical models, and the principles of Dynamical SystemTheory provide a solid foundation for these. The state-space approach makes it possible to accommodate disturbances and avarying environment. In this paper, an attempt has been made to develop a dynamic growthmodel based on the limited data, consisting of three annual measurements, collected from 22 teak sample plots in Karnataka, Southern India. A biologically consistent whole-stand growth model has been presented which uses thestate-space approach for modelling rates of change of three state-variables viz., dominant height,stems per hectare and stand basal area. Moreover, the model includes a stand volume equationas an output function to estimate this variable at any point in time. Transition functions werefitted separately and simultaneously. Moreover, a continuous autoregressive error structure isalso included in the modelling process. For fitting volume equation, generalized method of moments was used to get efficient parameter estimates under heteroscedastic conditions. A simple model containing few free parameters performed well and is particularlywell suited to situations where available data is scarce.


Background
About the species and its distribution Teak (Tectona grandis L. f.) is one of the most important tropical timber species and is suitable for multiple enduses. The potential for growing and managing teak in different ecological zones and under different situations is being increasingly recognized, leading to intensive domestication and cultivation of the species in countries/regions beyond its natural habitat (Perez and Kanninen, 2003).
Teak occurs naturally in parts of India, Myanmar, Laos and Thailand. It has been naturalized in Java, where it was introduced about 400-600 years ago (Kadambi, 1972;White, 1991). Early introductions of teak outside Asia were made in Nigeria in 1902, with the first plantations being of Indian origin and subsequently of Burmese origin (Horne, 1966). The first pure teak plantation in Tropical America was established in Trinidad in 1913. Teak planting in Honduras, Panama, and Costa Rica started between 1927and 1929(Ball et al., 2000. Teak is the world's most cultivated high-grade tropical heartwood, covering approximately 6.0 million hectares worldwide (Bhat and Hwan Ok Ma, 2004). Of this, about 94% are in Tropical Asia, with India (44%) and Indonesia (31%) contributing the bulk of the resource. Other countries like Thailand, Myanmar, Bangladesh and Sri Lanka contribute significantly with 17% in total. About 4.5% of the teak plantations are in Tropical Africa and the rest are in Tropical America, mostly in Costa Rica and Trinidad and Tobago (Pandey, 1998).
The most important teak forests in India are in Madhya Pradesh, Maharashtra, Karnataka, Tamil Nadu and Kerala, and also in Uttar Pradesh, Gujarat, Orissa and Rajasthan (Troup, 1921). Plantations have also been made in Haryana, West Bengal, Assam, Meghalaya and Dadra and Nagar Haveli (Chakrbarti and Gaharwar, 1995). However, many plantations are in a state of neglect.

Management of teak forests in Karnataka
The teak forests in India are managed based on the approved working plans (generally prescribed for ten years). The rotation period prescribed for teak in Karnataka is usually 120 years. The following mechanical (systematic; row thinning) and silvicultural (selective) thinning regimes are prescribed by the Karnataka State Forest Department.
An elite thinning at the end of 80th year (retaining about 150 trees ha −1 ).
However, this working plan prescription is not strictly followed in most cases due to various reasons. In some areas, trees are damaged bythe elephants and these damaged trees may be extracted every year.

Growth models
Reliable growth models are essential for effective longterm planning and decision making. This is especially important in intensively managed forest plantations, where it is necessary to evaluate alternative planting densities, thinning regimes, and rotation lengths (Garcia et al., 2011). Empirical models describe and generalize observed stand behaviour. They can be accurate and highly successful in fulfilling many of the forest manager's needs (Burkhart, 2008). However, they require extensive longterm experimental data, and extrapolation to unobserved situations can be uncertain. A theory-based model would describe hypothesized behaviour by reasoning from first principles. Theory-influenced semi-empirical models attempt to be compatible with observations, while at the same time use theoretical knowledge and hypotheses to improve performance under conditions for which data are scarce (Mohren and Burkhart, 1994).
In this paper, we present a biologically consistent whole-stand growth model for teak using data from the Karnataka State of Peninsular India. The main motivation was to produce a simple and robust base-line model for teak stands. The model uses the state-space approach for modelling rates of change of three state-variables: dominant height, number of trees per hectare and stand basal area. Moreover, the model includes a stand volume equation as an output function to estimate this variable at any point in time.
This research presents some new methodological contributions for developing growth models based on the state-space approach and differentiate from the past work using similar approach (e.g. García et al., 2011) in the sense that García et al. (2011) fitted each transition function separately, however, in this study the transition functions for mortality and basal area were fitted simultaneously using the full information maximum likelihood. Also, autoregressive procedure has been used in the present study to remove the problem of autocorrelation inherent in the dataset from repeated measurements. Further, volume equation was fitted using generalized method of moments to obtain efficient parameter estimates under heteroscedastic conditions even without estimating the heteroscedastic error variance.

Data
The data used in this study were collected on 27 pure and even-aged Teak research plots of different ages (11 to 38 years) and densities (498 to 2061 trees ha −1 ) established in 2010 in stands representing five forest divisions in Karnataka. The annual rainfall in these areas varies from 1600 mm to 4500 mm. The mean annual minimum and maximum temperatures vary from 11°C to 38°C. The plots are located throughout the teak growing areas in the state and were selected to represent the existing range of ages, stand densities and sites. The plots are rectangular and their size ranged from 220 to 436 m 2 , depending on stand density, in order to achieve a minimum of 30 trees per plot. The diameter at breast height (dbh, 1.37 m above the ground) was measured to the nearest 0.1 cm with digital callipers in all trees in the plots. Total height was measured to the nearest 0.1 m with a digital hypsometer. Age was calculated from the year of planting. The plots were re-measured annually during 2010-2013 and a total of three annual measurements were available for each plot.
The plot data includes: stand age (t, years), density (N, trees ha −1 ), quadratic mean diameter (D, cm), average height h; m À Á , dominant height (H, m), defined as the mean height of the 100 thickest trees per hectare, basal area (B, m 2 ha −1 ) and total stand volume (V, m 3 ). Total stand volume was calculated from a compatible tree volume equation (Tewariet al., 2013). The summary statistics of these stand variables are shown in Table 1. Although, the mechanical (systematic) and silvicultural (selective) thinnings are prescribed in the forest management working plans, however, thinning was carried out only in one stand. Moreover, an exhaustive examination of the data was carried out to reject sample plots on which illegal logging had been detected to make sure that the reduction in number of trees observed were due to natural mortality.

Model structure
Out of the 27 plots, 5 plots, on which illegal logging had been detected, were excluded from the analysis. The dynamic growth model developed from these data is based on the state-space approach and it is similar to those of García (2011), García et al. (2011) and García (2013). In this model it is assumed that the behaviour of any stand of teak evolving in time can be approximated by describing its current state with four state variables: dominant height (H), number of trees per hectare (N), basal area (B), and a measure of site occupancy (Ω), using transition functions to estimate the change of states as a function of the current state of the variables.
The transition functions are used to predict the growth by updating the state variables, ensuring two natural properties (García, 1994): (i) consistency, meaning no change for zero elapsed time; (ii) pathinvariance, where the result of projecting the state first from t 0 to t 1 , and then from t 1 to t 2 , must be the same as that of the one-step projection from t 0 to t 2 . Transition functions generated by integration of differential equations (or summation of difference equations when using discrete time) satisfy these conditions and allow computing the future state trajectory. These properties can also be achieved by using techniques for dynamic equation derivation known in forestry as the Algebraic Difference Approach (ADA; Bailey and Clutter, 1974) or its generalization (GADA; Cieszewski and Bailey, 2000).

Transition function for dominant height
The two-site-specific parameter equation, derived from the Korf base model and proposed by Tewari et al. (2014), was used as transition function for dominant height. The GADA model allows simultaneous concurrent polymorphism and multiple asymptotes, two characteristics of site equations that are sometimes desirable (Cieszewski, 2002). The mathematical expression of this model is where, H 1 and t 1 represent the predictor height (meters) and age (years), and H 2 is the predicted height at age t 2 .

Transition function for mortality
Natural mortality or survival is extremely variable and, particularly difficult to predict. The mortality transition function is based on assuming that the rate of change of N relative to dominant height increment depends on the current values of H and N: where, N is number of trees ha −1 , H is the dominant height and a i are parameters to be estimated. Grouping the terms of equation (2) and integrating both sides gives the following invariant: Equating the invariant for times 1 and 2 gives the following transition function for number of trees ha −1 : This predicts N 2 at a height H 2 given any initial values (H 1 ; N 1 ). Equation (4) was fitted using two different approaches. In the first case, no parameter restrictions were assumed while in the second case, parameter restrictions were obtained by assuming limiting values for the slope of the self-thinning lines based on the average square spacing S ¼ 100= ffiffiffiffi N p (Garcia, 2009). The rate of change of average square spacing relative to dominant height increment depends on the initial values of the same variables: Integrating equation (5) and taking logarithms, the following invariant is obtained: Any limiting value assumed for 2γ/α in equation (6) represents the slope of the self-thinning line assumed for the species. Comparing with equation (3), the following relationship is obtained: This relationship could be useful to estimate the parameters of the mortality transition function when the available data is limited, to assure a biologically consistent behaviour (García, 2009). Values of 2 and 3 were assumed for the -2γ/α ratio, i.e. for the slope of the selfthinning line.

Transition function for basal area
Instead of predicting basal area (B) directly, we choose to model growth of the product W = BH. The rate of change of W can be expressed as the difference between two components: the gross increment and the mortality. In pure and even-aged stands (such as those considered in this article), the gross increment can be written as b 1 ΩHN b 2 , where Ω is a relative occupancy factor that reduces growth in young or recently thinned stands that do not still fully occupy the site. The mortality as −k W N dN dH ¼ −kW d log N dH , where k represent the mean size of dying trees relative to the mean size of the survivors, considered as constant. Therefore, the general W growth model is: A simple closed-form solution of the differential equation (8) can be achieved in the special case where b 2 is equal to k: Intuitively, Ω represents resource (e.g. light, nutrient, moisture) interception. In the beginning, occupancy is low in young stands, and gradually increases approaching 1 when canopy closes. We assume that the rate of closure initially depends on H in the same way as the gross increment, and decreases to zero as full closure is approached.
Integration gives the invariant: With this occupancy model, integration of equation (9) gives the transition equation for basal area: with In order to understand the methodology and the derivation of the transition equation for basal area, the integration of the differential equation (9) has been provided as Appendix.

Total stand volume equation
Volumes per hectare can be estimated given the state variables using an output function. Two different total stand volume equations were fitted using as independent variables dominant height and basal area.
The models were fitted using the generalized method of moments (GMM) in the MODEL procedure of SAS/ ETS® (SAS Institute Inc., 2008). This method produces efficient parameter estimates under heteroscedastic conditions, without estimating the heteroscedastic error variance (Greene, 2012). The main drawback of this method is that error estimates for the predictions can not be generated without specifying the error structure (Parresol, 2001).

Parameter estimation and model selection
The model was fitted based on projections over pairs of consecutive measurements. Since the transition function for dominant height was fitted previously, the observed dominant heights in the database were substituted by the predicted values of this model in order to include the error of estimations in the whole-stand growth model. It was considered convenient to initialize the state variables at breast height, where H = 1.37 and B = 0 are known. The number N b of trees per hectare at this time was estimated by applying an 82.5% survival (since the survival rate observed after one year of planting in these stands was 80-85%) to the density from the spacing adopted at the time of planting which was 2 × 2 m except for one plot in which it was 1.5 × 1.5 m. Following Garcia (2011), the occupancy was initialized at breast height as: where, c i are parameters to be estimated; c 2 represents a breast-height density at which the stand would be fully closed. García et al. (2011) fitted each transition function separately. Here the transition functions for mortality and basal area were fitted simultaneously using the full information maximum likelihood (FIML) in order to obtain estimates that are consistent and efficient and contribute in increased precision of model predictions.
Since the convergence with this approach is very sensitive to the starting values of the parameters, both transition functions were first fitted separately to get the best set of initial values.
Individual tree forest growth and yield models usually employ a set of equations to describe stand development over time. It is not unusual in the forestry literature to treat the same variable as dependent in one equation and independent in another. When a variable is used as both an endogenous (dependent) variable on the lefthand side of one equation and on the right-hand side of another equation, this renders the system of equation simultaneous (Goldberger, 1964). In the present study, mortality is used as dependent variable in equation (4) and as an independent variable in equation (12), making the system of equations simultaneous. Simultaneous regression techniques lead to more efficient estimators. A gain in efficiency increases the precision of the resulting model predictions. The gain in efficiency is higher when the errors among different equations are highly correlated (Judge et al., 1988).
The error was modelled using a continuous autoregressive error structure (CAR(×)), which accounted for the time between measurements. Age at breast height was estimated using an iterative procedure based on the transition function of dominant height previously described. Different ages at breast height were used to estimate the dominant heights at the age of plots measurements and the age which minimized the sum of squares of the differences between observed and predicted dominant heights was selected as age at breast height. Those ages were only used to obtain correct consistent estimates of the parameters and their standard errors.
Estimation of the parameters was carried out with the MODEL procedure of SAS/ETS® (SAS Institute Inc., 2008). The CAR(×) error structure was programmed in this procedure which allows for dynamic updating of the residuals.
The comparison of the estimates for the different models was based on numerical and graphical analyses of the residuals. Two statistical criteria were examined: adjusted model efficiency (ME adj ) and root mean square error (RMSE). The expressions of these statistics are summarized as follows: where, y i ,ŷ i and y are the observed, estimated and mean values of the dependent variable, n is the total number of observations used to fit the model, and p is the number of model parameters.
In addition to the validation of individual components of the dynamic growth model, overall validation of the whole model was carried out. All possible combinations of inventories were used and the observed H, N and B values from the first inventory were used to estimate total stand volume at the age of the last inventory, including all the components of the whole stand model. In order to evaluate whether the model performs acceptably well when used for estimating total stand volume, a critical error, expressed as a percentage of the observed mean, was computed (Reynolds, 1984): where, n is the total number of observations in the data set, y i,ŷi and y are the observed, predicted and mean values of total stand volume, and τ is a standard normal deviate at the specified probability level (τ = 1.96 for α = 0.05). χ 2 n is obtained for α = 0.05 with n degrees of freedom. In addition to using this statistic, the plot of observed against predicted values of the total stand volume was also inspected.

Results and discussion
At first, the transition functions for mortality and basal area were fitted separately using nonlinear least squares and without assuming autocorrelation. However, a slight trend in residuals as a function of lag1 residuals within the same sample plot was apparent in all the models. Then, all the models were refitted using a first-order continuous-time autoregressive error structure. After applying this correction, the error trends in residuals disappeared (Figure 1).

Separate fitting
Three transition functions for mortality were analyzed and the results are presented in Table 2. The first one is based on equation (4) without assuming any restrictions for the parameters and the other two are based on assuming values of 2 and 3 for the slope of the selfthinning line according to equation (7). The three models explained more than 98% of the total variance, calculated without considering the initialized point at breast height, and the model with a value of 2 for self-thinning line resulted in the best values of the goodness-of-fit statistics.
Natural tree mortality is a complex process that is neither constant in time nor in space, so it is difficult to predict or explain the factors that control it (Van Laar and Açka, 2007). Data from permanent sample plots frequently show that a relatively large part of the plots have no occurrences of mortality even over periods of several years and sometimes an important reduction in trees per hectare occurs in short periods of time (e.g. Monserud and Sterba, 1999;Eid and Tuhus, 2001). Therefore, for a correct mortality modeling, an adequate combination of short and long intervals between measurements is essential. The time interval between consecutive measurements of the database used in this study is one year and the models fitted show a good performance for these short periods, within the age interval covered in the sample plots. However, the behavior of the models for long periods projected from the breast age showed high levels of mortality ( Figure 2).
The results of the fitting of the basal area transition function are shown in Table 3. It was not possible to obtain convergence or reasonable estimates with all the parameters. Therefore, as in Garcia et al. (2011), we fixed c 3 in equation (16) with the value 2.4 and c 2 in the same equation with the value 10000, implying full canopy closure at breast height with 10000 trees ha −1 . It had been shown that these values are not critical (García, 2011;García et al., 2011).
In a first step, the model (Equation 12) was fitted without assuming any other restriction for the remaining three parameters (k, b 1 and c 1 ), however, unrealistically small values for parameters k and c 1 (0.2121 and 0.0369, respectively) were obtained. In a second step, the transition function was fitted assuming a value of 0.4 for parameter k (García, 2011), however, the graphical representation indicated that the estimated value of c 1 was unrealistic again (0.0239). Finally, the model was fitted assuming k = 0.4 and two different values for c 1 , 0.05 and 0.01. Figure 3 shows occupancy curves for various values of c 1 generated with equation (13), starting from 2062 trees per hectare at  breast height. Also shown is the relative closure, R, which can be seen as representing the amount of foliage and roots relative to that of a closed stand. Closure and occupancy are related by Ω ¼ 1− 1−R ð Þ c 3 . Both selected values for c 1 , 0.05 and 0.1 seem reasonable, implying nearly closed canopies somewhere between 5 and 10 m top height. Or seen in another way; an equilibrium canopy depth of 5 to 10 m. The best goodness-of-fit statistics was obtained constraining c 1 to 0.05 and this value was used in the simultaneous fitting. The model explained more than 99% of the total variance, calculated without considering the initialized point at breast height.

Simultaneous fitting
The results of the simultaneous fitting of both transition functions (Equation 4 for mortality and Equation 12 for basal area) are shown in Table 4. Finally, two total stand volume equations were fitted separately. All the parameters were significant except for the intercept of equation (14). Both models performed equally good, although the best statistical criteria were obtained with equation (15) which explained more than 97% of the total volume variability with a root mean square error of 13.34 m 3 ha −1 . The final expression of the dynamic model is as follows: Evaluation Critical errors of 1.11%, 3.45% and 8.9% were obtained for total stand volume at time intervals of 1, 2 and 3 years. The critical errors obtained with the proposed model are lower compared to that obtained for projecting total stand volume in Scots pine (15.3%) in Galicia (Diéguez-Aranda et al., 2006) and the ones obtained for the same stand variable in radiata pine (10.9%, 11.9% and 17.3%) in Galicia, considering time intervals of 3, 6  and 9 years, respectively (Castedo-Dorado et al., 2007). However, the models used for Scots pine and radiata pine in Galicia included projection intervals longer than those used in this study. Considering the required accuracy in forest growth modelling, where a mean prediction error of ±10-20% is generally acceptable (Huang et al., 2003), it may be stated that the dynamic model provides satisfactory estimates within the range of ages and the projection intervals observed. A plot of observed against predicted values of total stand volume, estimated using the whole dynamic stand growth model, is shown in Figure 4. The linear model fitted to the scatter plot behaved well (R 2 = 0.9743), although there appears to be a slight tendency towards underestimation for low volumes and overestimation for high volumes.

Conclusions
Forest management planning relies heavily on mathematical models that involve time. Modern dynamical system theory provides a framework for a flexible representation of varying environments, as well as of responses to intensive silviculture and natural disturbances.
A biologically consistent whole-stand growth model has been presented for teak stands in Karnataka State of India. The transition functions for mortality and basal area were fitted simultaneously, using the full information maximum likelihood, which produced consistent and efficient estimates contributing to increased precision of model predictions. Also, autoregressive procedure was  The goodness-of-fit-statistics were calculated without considering the initialized point at breast height. *indicated the parameter value was fixed in the model, therefore, it was not estimated.