A preliminary growth and yield model for Eucalyptus globoidea Blakely plantations in New Zealand

Background: New Zealand’s plantation forest industry is dominated by the exotic species radiata pine (Pinus radiata D.Don), which comprises approximately 90% of the net stocked area. However, there is interest in introducing new species to: (a) provide wood that is naturally decay-resistant as a substitute for wood treated with preservatives; (b) match species to the wide variety of environmental conditions in New Zealand; and (c) reduce reliance on P. radiata. Some Eucalyptus species are considered as potential alternatives to P. radiata, specifically those that can survive in resource-limited conditions and produce high quality wood. While Eucalyptus species are grown in plantations in many regions of the world, limited information is available on their growth in New Zealand. Eucalyptus globoidea Blakley is of particular interest and has been planted in trials throughout New Zealand. A complete set of preliminary growth and yield models for this species will satisfy the initial information requirements for diversifying New Zealand’s plantation forest industry. Methods: A set of growth and yield models was developed and validated, based on data from 29 E. globoidea permanent sample plots (PSPs) located mostly in North Island and a few in South Island of New Zealand. Trees were measured at different time intervals in these plots, with height and diameter at breast height (DBH) ranging from 0.1–39.8 m and 0.1–62.3 cm, respectively. An algebraic difference approach (ADA) was applied to model mean top height, basal area, maximum diameter, and standard deviation of DBH. Non-linear regression equations were used to project stand volume and height-diameter relationship, and Reineke’s stand density index (SDI) approach was employed to model mortality. Results: Mean top height, maximum diameter, and standard deviation of DBH were best fitted by Von Bertalanffy-Richards (SE=1.1 m), Hossfeld (SE=2.4 cm), and Schumacher polymorphic (SE=1.6 cm) difference equations, respectively. Basal area data were modelled with high precision (SE=6.9 m2 ha-1) by the Schumacher anamorphic difference equation. Reineke’s SDI approach was able to explain the self-thinning as a reduction in the number of stems per hectare. Stand-level volume per hectare and height-diameter relationship models were precise when including site-specific variables with standard errors of 40.5 m3 ha-1 and 3.1 m, respectively. Conclusion: This study presents a set of preliminary growth and yield models for E. globoidea to project plot-level growth attributes. The models were path invariant and satisfied basic traditional mensurational-statistical growth and yield model assumptions. These models will provide forest growers and managers with important fundamental information about the growth and yield of E. globoidea. New Zealand Journal of Forestry Science Salekin et al. New Zealand Journal of Forestry Science (2020) 50:2 https://doi.org/10.33494/nzjfs502020x55x


Introduction
A resilient forest economy would be diversified, with healthy forests of all ages producing a range of valuable products and services. The New Zealand forestry industry is almost entirely based on radiata pine (Pinus radiata D.Don) plantations (New Zealand Forest Owners Association 2017) due to its rapid growth rate across a broad range of sites (Turner et al. 2008) and established processing infrastructure and markets. However, there are opportunities to introduce new species and overcome some limitations of P. radiata (Millen et al. 2018). For example, non-naturally durable wood, less diversified forest ecosystems and slow growth in drought-prone environments limit radiata pine's potential. New Zealand's existing commercial forest sector could be complemented by introducing other species. This would help to reduce the reliance on largescale plantations containing pure stands of P. radiata which are at increased risk of pest and diseases attacks (Chou 1991) and produce a relatively narrow range of forest products.
Species of Eucalyptus have been considered as a commercial forestry alternative to P. radiata, especially those species that can grow well in dry conditions and produce high quality timber (Menzies 1995). However, despite strong advocacy for alternative species from various groups, only 1.2% of the total plantation forest area in New Zealand is comprised of various Eucalyptus species (MPI 2019). Growing Eucalyptus in New Zealand has, over the years, been challenging (Berrill & Hay 2005;Berrill & Hay 2006) due to specific site requirements such as sensitivity to soil moisture availability and frosty environments (Bell & Williams 1997;Williams & Woinarski 1997), pests and diseases that affect their health and productivity (Lin 2017), and a lack of markets for Eucalyptus wood products (Apiolaza et al. 2011). Recently, the situation has started to change, in part, because of the New Zealand Dryland Forest Initiative (NZDFI) and renewed consumer demand for Eucalyptus timber (Satchell & Turner unpublished data). The NZDFI has catalysed research into several naturally durable Eucalyptus species, chosen for their desirable properties (Nicholas & Millen 2012), for deployment on ex-pasture lands in relatively dry parts of the country (NZDFI 2013). Despite these advances, little is known about the growth dynamics of many of these Eucalyptus species in New Zealand.
Managed forests are dynamic biological systems that change in response to environments and silvicultural practices. Growth and yield models can support effective decision making by describing current and future forest dynamics (Blake et al. 1990;Blanco et al. 2005;Castedo-Dorado. et al. 2007;Clutter et al. 1983). Traditional time-based growth and yield models, called mensurational-statistical models, provide robust growth predictions but give little information about the mechanisms of forest dynamics (Korzukhin et al. 1996). Apart from being mathematically simple and biologically rational, Clutter et al. (1983) noted several important features: i) representation of growth and yield should be compatible; ii) the functions should be consistent; iii) the functions should be path-invariant; and iv) the functions should rise to asymptotes. Mensurational-statistical forest growth models are often based on large datasets, comprising repeated field measurements in permanent plots Pienaar & Rheney 1995) or information obtained from remotely-sensed data (Battaglia et al. 2004;Landsberg et al. 2003).
However, in scenarios where comprehensive data are not available, it may still be desirable to develop preliminary growth and yield models to forecast forest growth (Vanclay 2010), especially for new species (Berrill et al. 2007;Kitikidou et al. 2016;Palahí & Grau 2003). Such models are often imprecise, but can be useful (Box 1976) to obtain an initial forecast in order to make decisions about establishment, tending, and potential log marketing. Preliminary models are not only useful for characterising stand development but also provide insights into the yield potential of sites, a crucial factor for sound management of any forest stand (Tewari & Gadow 2003). Moreover, preliminary mensurationalstatistical models can be easily implemented and used by forest managers to generate initial estimates of growth and yield.
While preliminary models are available for Eucalyptus fastigata, E. nitens, and overall stringy-bark groups in New Zealand (Berrill & Hay 2005;Berrill & Hay 2006), no growth and yield models exist for the Eucalyptus species within the NZDFI's programme. Despite being planted in trials and plantations around New Zealand, managers and growers have limited knowledge of the expected growth and yield for these species. Development of species-specific, stand-level preliminary models will not only give them more information but also guide them about species choice for planting and future management.
The NZDFI selected a set of naturally durable Eucalyptus species (see, Page & Sing 2014) based on Australian timber durability standard (Class 1 and 2). E. globoidea is one of the top ranked species in that list and is commonly classified into the stringybark group. It was sparsely planted around New Zealand prior to the NZDFI programme. This species has naturally durable wood and is considered a highly durable timber (class 1 or 2) in the Australian standards (AS5606-2005) (Nicholas & Millen 2012a). E. globoidea is well adapted to dry parts of the New Zealand. Moreover, a strong consumer demand for naturally durable Eucalyptus wood has been identified (Kakitani 2017). Growth and yield model functions must adequately describe the system at any point in time by allocating local transitions (Garcia 1988), that is the rate of change of state as a function of the current state and of the current values of external control variables. Therefore, the main objective of this study was to develop a preliminary stand-level E. globoidea growth and yield model. This model projects estimates of mean top height (MTH), basal area/ha (G), maximum diameter at breast height (D max ), standard deviation of diameter (SD D ), stand volume (V), self-thinning and height-diameter relationships (H-D) forward in time following measurements of stands. Then D max and SD D can be used to fit a reverse Weibull function to describe the stand-level diameter distribution (García 1981;Kuru et al. 1992). Individual tree models were considered, but diameter distribution models have been shown to be superior to individual tree models when long-term projections are required by those planning harvests many years in the future (Methol 2001).

Data preparation and description
Tree-and plot-level E. globoidea data were available from a nationwide permanent sample plot system (Pilaar & Dunlop 1990). Data from twenty nine permanent sample plots (PSPs) established in plantations at ten different localities were available (Table 1 and Figure 1).
Trees were measured in the PSPs at 1 to 10-year intervals with an irregular frequency. Mean top height (average height of the 100 largest diameter stems per hectare) (MTH) and maximum diameter (the largest diameter measured at right angles to the stem over the stubs) (D max ) of the trees were calculated from the individual tree measurements by following the procedure described by Goulding (2005). The standard deviation of DBH (SD D ) was calculated for each plot measurement. Basal area (G) was calculated by summing the crosssectional area at breast height (1.4 m) of all trees in the plot, then dividing by plot size to provide a per hectare estimate. Stand volume (V) was calculated within each sample plot by estimating and summing up individual tree volumes calculated using a simplified E. globoidea stem taper volume equation presented in Lundgren (1995) (see the Appendix).
Plot-level summary data were organised by representing all possible measurement time intervals and used to fit differential equations. Simple time increment data from plot-level summaries were used to fit volume-per-hectare equations. The height-diameter relationship was modelled using individual tree measurement data from all plots.

Modelling and evaluation
The algebraic difference approach (ADA) (Bailey & Clutter 1974) was used to model mean top height (MTH), basal area/ha (G), maximum diameter (D max ) and standard deviation of diameter (SD D ). Sixteen well-known and frequently used polymorphic and anamorphic forms of differential equations (Bailey & Clutter 1974;Belli & Ek 1988;Ek 1974;Vanclay 1994;Zeide 1993) (Table 2) were fitted to the data using non-linear least-squares (Clutter 1963).
Volume per hectare yields (V) were modelled using various simple, established and commonly used functions (Table 3). Height-diameter (H-D) models were developed by fitting the Näslund (1936) equation with a range of different exponent terms (Zhao 1999): where H is tree height (m), D is diameter (cm) at breast height (1.4 m), and α and β are model parameters. The exponent term (γ) here is variable. This function is widely used and can be expressed in a linear form: The linearity of the height-diameter relationship is a unique property at a plot-level or at a stand-level (Curtis 1967;Garcia 1974) or at a stand-level (Zhao 1999) when a few plots are sampled from the same stand in a single site. However, fitting H-D relationships at a stand level results in underestimates of variability in MTH (Mason 2019), so sufficient heights (usually 12) were measured in each PSP to allow the fitting of plot-level H-D relationships. A better height-diameter relationship can  Weibull 1 5 Weibull 2 6 Hossfeld 7 Von Bertalanffy-Richards 1 8 Von Bertalanffy-Richards 2 9 Von Bertalanffy-Richards 3 10 Anamorphic form

Schumacher A1 11
Schumacher A2 12 Gompertz 13 Von Bertalanffy-Richards 14 Weibull 15 Hossfeld 16 be obtained by identifying and incorporating relevant factors accounting for differences within stands. These could include tree species, stand age, site characteristics, genetics, stocking, and silvicultural treatment (Zhao et al. 2006). This was achieved by separating and linearly expanding the regression coefficients described in Woollons et al. (1997). Generally, data on self-thinning or mortality follow a binomial or Poisson distribution. Therefore, development of models to describe this process needs larger datasets than for other models. Due to the small number of plots, a conceptual self-thinning/mortality model was produced by applying Reineke's stand density index (SDI) (Reineke 1933). This was done by estimating quadratic mean diameter at breast height (DBH) and basal area (G). The maximum SDI was assumed from the original range by selecting the highest stand density as there was no specific evidence of self-thinning in the dataset.
For validation there was no independent dataset available for this study, nor was the dataset large enough to be subdivided into fitting and validation datasets. Therefore, model validation was carried out by the 'leave-one-out' method of cross-validation (LOOCV), a method which is also called "Jackknife" (Arlot & Celisse 2010). Thus, the models were fitted times, leaving out each sample plot once, so that the number of model fits was equal to the number of plots (Sánchez-González et al. 2005), and residuals of predictions for the plots left out were compared with those of the overall model fit.
All the models except the self-thinning model were evaluated through the validation procedure. Validation included a visual analysis of graphs of the residuals, the calculation of root mean square error (RMSE) (Equation 23), mean absolute error (MAE) (Equation 24), bias (Equation 25), and adjusted coefficient of determination (R 2 ) (Equation 27). Adjusted R 2 values were not considered for assessing differential equations as it is sensitive to grouped and repeated data (Warren 1971). The predictive ability of the models was evaluated using prediction errors or predictive residual error sum square (PRESS) statistics (Equation 28), where O i is the observed value, P (i,-i) is the estimated value for observation i (where the latter is absent from the model fitting) and n is the number of observations. Each model has n PRESS residuals associated with it, and the PRESS (Prediction sum of square/P-square) statistic is defined as (Myers & Myers 1990): The bias and precision of models were analysed by computing means of the PRESS residuals. All statistical analysis was performed in the R statistical environment (R Development Core Team 2017). Different non-linear regressions were fitted using the "nls" function in the base package with appropriate significant variables. Evaluation metrics "adj. R 2 ", and "rmse", "mae", "bias" functions were used from the "Metrics" package (Hamner & Frasco 2018). Residuals were visually inspected for their normality and variance homogeneity. All graphical analyses were performed with the "ggplot2" package (Wickham 2016).

Mean top height (MTH) model
Among the tested differential equations (Table 2), the first Von Bertalanffy-Richards polymorphic model (Equation 8) exhibited the most precise fitting statistics based on goodness-of-fit. It minimised bias and standard error of prediction compared with the other models tested. However, the RMSE and MAE were higher in model fitting statistics, relative to validation, at which time they roughly halved to 3.8 m and 2.5 m, respectively ( Table 4). The model residuals were well distributed with minor heteroscedasticity at the beginning of the modelling period (Figure 2A and 2B). The model predictions covered the entire range of measured MTH values, except those for two stands ( Figure 2C). These two stands contained outliers, which were left out from the initial model building procedure. Model parameters are provided in the appendix (Table A2).

Basal area per hectare (G) model
Among tested models (Table 2), the first anamorphic Schumacher model with a single parameter (Equation 11) was found to be the best fit for basal area/ha projection. This model had the lowest error and greatest precision. Precision increased during validation with much less error indicating stable model performance (Table 5). The residual plot exhibited minor heteroscedasticity ( Figure  3A). The residual distribution was positively biased, which indicated a slight overprediction. Moreover, the model predicted basal area values covering the measured range, except for two stands (Figure 3). Model parameters are provided in the Appendix (Table A2).

Maximum diameter (D max ) model
The Hossfeld polymorphic model (Equation 7) predicted the maximum diameter (D max ) with greatest overall precision and least bias in comparison with other candidate model forms (Table 2). In this case, RMSE and MAE increased from fitting to validation statistic and bias went from positive to negative (Table 6), which indicated model under-prediction during validation. However, the standard error (SE) reduced slightly in validation indicating higher precision. The low MPRESS and MAPRESS values also indicated model goodnessof-fit (Table 6). Residuals were highly biased at the beginning and end of the modelling period ( Figure 4A), consistent with the limited availability of data. However, they were normally distributed ( Figure 4B). The function for predicting D max enveloped all the measurements and followed a sigmoid shape ( Figure 4C), which ensures  biological rationality. Model parameters are provided in the Appendix (Table A2).

Standard deviation of diameter (SD D ) model
Among all the candidate models tested (Table 2), the standard deviation of diameter (SD D ) was best predicted by the second Schumacher polymorphic model (Equation 2). The model had the lowest prediction errors. The RMSE (1.5 cm to 1.9 cm) and MAE (1.2 cm to 1.5 cm) increased slightly from fitting to validation. Minimal MPRESS and MAPRESS values also confirmed precision of the model (Table 7). Model parameters are provided in the Appendix (Table A2). Graphically, the model predicted values and residuals appeared to follow a normal distribution ( Figure 5A). The residuals plot shows overprediction and positive bias of the model with few outliers in the frequency distribution plot ( Figure 5B). The prediction plot shows that the model included the full range of measured SD D ( Figure 5C).

Volume per hectare (V) model
The best performing model of all those tested for volume per hectare yield (Table 3) was the four parameter Jansen et al. (1996) model (Equation 20). Fitting statistics showed minimal prediction error and high precision, though validation statistics were greater in both cases (Table 8), and these were confirmed by small MPRESS, small MAPRESS, and high adjusted R 2 value (Table 8). These results are also confirmed graphically (Figure 6), although there is a minor heteroscedastic tendency in the residuals ( Figure 6B). Model parameters are provided in the Appendix (Table A2).

Height-diameter (H-D) model
The stand-specific individual height-diameter (H-D) model gave precise predictions with an exponent of -2 (Equation 30). Stand-specific elevation and basal area (G) were found to influence the H-D relationship significantly (P<0.05) and adding them into the final model improved the fit. The goodness-of-fit values increased slightly from fitting to validation statistics, which indicated less precision in prediction (Table 9). Residuals were normally distributed, and the model fitted well (Figure 7B and C). Model parameters are provided in the Appendix (Table A3). (30)

Self-thinning model
The self-thinning model based on Reineke's SDI fitted the limited available data well. Stand density ranged from 150-1350 stems ha -1 , with most plots having a stand density between 400 and 650 stems ha -1 ( Figure 8B). The maximum carrying capacity was calculated as 1350, 25 cm diameter trees per hectare. Based on this value, a density management diagram was produced ( Figure  8A) with lines indicating understocking below 35% of maximum carrying capacity, full stocking between 35% to 55% of maximum carrying capacity, and over-stocking above 55% of maximum carrying capacity. Natural mortality started to occur when stocking approached the maximum carrying capacity ( Figure 8A).

Discussion
This study developed a preliminary set of stand-level growth and yield models for E. globoidea in New Zealand using sparsely available data. The state of a plot was adequately described by the following state variables: mean top height, basal area/ha, volume/ha, stocking, maximum diameter, standard deviation of DBH and  = 1.4 + (( 0 + 1 × ) + ( 0 + 1 × ) ) −2 the height-diameter relationship. The nature of the plot's growth is described by the rate of change of these variables over time by their corresponding transition function.
The final models were the best-fitted models, which generally had the highest accuracy among the tested set of equations from several differential forms. The best MTH, D max and SD D models took polymorphic forms, similar to earlier preliminary modelling studies in a range of species, including even-aged Cupressus lusitanica Mill. and C. macrocarpa Hartw. plantations (Berrill 2004), Acacia melanoxylon R.Br. (Berrill et al. 2007), Eucalyptus fastigata (Berrill & Hay 2005) in New Zealand and Pinus nigra Arn. in Catalonia, Spain (Palahí & Grau 2003). However, basal area/ha (G) was best fitted with an anamorphic form, which is unusual but can be found in similar types of data-limited situations. For example, Vanclay (2010) suggested one parameter anamorphic forms to deal with a similar small dataset.
Overall, model projections followed the growth pattern of this species. For example, the mean top height for a 15-year-old stand ranges from 15-22.5 m and basal area ranges from 14.5 to 62.5 m 2 ha -1 . Apart from a few outlier stands, our results followed similar trends to those reported by Nicholas and Millen (2012b). Similar to our study, their study and models developed in it were based on a very small number of measurement plots. Furthermore, Meason et al. (2016) reported that most of the small-scale plantation PSPs resided in New Zealand's North Island. This may affect the model's capability to perform over a wider range of environmental conditions. There were some errors in model prediction, which may be due to the irregular measurement intervals for the stands included in the study. Lee (1998) reported that long measurement intervals can produce apparently larger errors than short measurement intervals, but longer intervals were vital for avoiding biased projections over long intervals. The measurement  periods also differed among the PSPs, which may have caused bias and heteroscedasticity through the modelling period (Lee 1998). Furthermore, model precision could likely have been improved by reinforcing it with more biological, or silvicultural information, for example, thinning information or any kind of natural disturbance events (Park & Wilson 2007). In this study, such information was available for only a small number of plots (see Table A1 in the Appendix). Borders et al. (1988) reported autocorrelation in data while using similar types of datasets to those used in this study, especially in a data-limited situation. When data covered all age classes as well as sites, such autocorrelation can be tested independently by separating each time interval (Borders 1989), which was not possible in this case. However, this study aimed to use the available data to condition the shape of a previously-known process. Also, all these models are based on mensurational equations and could receive further reinforcement from a biological perspective, by adding physiology into the modelling framework.
The self-thinning model was based on the SDI concept of Reineke (1933) with the "imminent competitionmortality" theory of Drew and Flewelling (1977). Here, competition related mortality likely occurs within a zone defined by two lines: the maximum size-density relationship (100% relative density) and a second line paralleling the first at lower densities for the same mean  tree size (55% relative density). According to this theory, a uniform stand may have imminent competitionmortality within the zone, but the probability of that could be lowered by substantially lowering the density. Conversely, if a stand is allowed to grow for many years within the zone mortality will occur (Drew & Flewelling 1979). Therefore, the SDI approach can only give an indication but not any causal explanation, so mortality cannot be precisely predicted on this basis (Drew & Flewelling 1979).
The translation of specific management objectives into appropriate upper and lower levels of growing stock is the key and most critical step to design a density management regime (Long 1985). Moreover, stands in the dataset must exhibited self-thinning to fit these lines, which was not the case in this study. However, this study reported no live stems above the maximum line which indicated the SDI approach's applicability. Nevertheless, while the SDI approach can be easily estimated and applied (Long 1985), it requires further testing and refinement with more data, particularly from older, higher-stocked plots where self-thinning is evident. A more comprehensive dataset would enable the slope (power term in the SDI function) of the size-density relationship to be investigated. Both Pretzsch and Biber (2005), and Saunders and Puettmann (2000) reported that the SDI function's power exponent term changed with species and site, but in this study the default value (1.605) was used due to a lack of data.
These preliminary models offer an indication of how E. globoidea may grow in New Zealand. They can be useful tools for forest managers to make initial management decisions for mature E. globoidea. For example, when applied in combination these set of models predict that 15-year-old stands will have mean top heights ranging from 15-22.5 m, basal area/ha ranging from 10-58 m 2 ha -1 and volume ha -1 ranging from 52-252 m 3 . However, the set of models presented here did not cover all age classes, and so some interpolation or extrapolation may occur during projection. Silvicultural and natural disturbances were not accounted in this study; therefore, the models' performance can be altered once such effects are considered. The models are only strictly applicable to the site conditions specific to the plots used to develop them, hence need to be validated with data from plots in new sites. Due to small amounts of data, especially the number of PSPs, the models must be regarded as preliminary and must not be used beyond the data range of the fitting dataset.

Conclusions
The models developed in this study will provide valuable information and understanding about the growth patterns, stand density dynamics, and potential yield of E. globoidea stands in New Zealand. This information substantially increases the limited knowledge base about this species, thus will help growers, managers and investors to make appropriate planting and management decisions, as well as the potential economic return at harvest. The growth patterns will vary at individual sites; therefore, caution must be exercised how these results are applied. Moreover, more plot measurement data including site characteristics and silvicultural regimes may increase the precision of these models and reduce bias in future.
This study also showed sparse dataset can be useful to make indicative prediction models, which could give a preliminary information about specific species. However, acquisition and maintaining longterm plot measurement data is indispensable to make comprehensive forecasting models needed to underpin management decisions.

Software
A web-based E. globoidea growth and yield software is available from the following link: http://www. treesandstars.com/models/egloboidea.htm