A preliminary growth and yield model for Eucalyptus globoidea Blakely plantations in New Zealand Serajis Salekin1,2,*, Euan G. Mason1, Justin Morgenroth1, Dean F. Meason2 1 New Zealand School of Forestry, University of Canterbury, Christchurch 8140, New Zealand 2 Scion, 49 Sala Street, Private Bag 3020, Rotorua 3010, New Zealand *Corresponding author: serajis.salekin@canterbury.ac.nz (Received for publication 14 May 2019; accepted in revised form 14 May 2020) Abstract 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 E-ISSN: 1179-5395 published on-line: 27/05/2020 © The Author(s). 2020 Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. Research Article Open Access Keywords: algebraic difference approach (ADA), basal area, Eucalyptus globoidea, growth and yield model, mean top height, mensurational-statistical model, stand density index. Salekin et al. New Zealand Journal of Forestry Science (2020) 50:2 Page 2 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 large- scale 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 (Castedo-Dorado et al. 2007; 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 mensurational- statistical 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 (Dmax), standard deviation of diameter (SDD), stand volume (V), self-thinning and height-diameter relationships (H-D) forward in time following measurements of stands. Then Dmax and SDD 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). Methods 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) (Dmax) of the trees were calculated from the individual tree measurements by following the procedure described by Goulding (2005). The standard deviation of DBH (SDD) was calculated for each plot measurement. Basal area (G) was calculated by summing the cross- sectional 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 (Dmax) and standard deviation of diameter (SDD). 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): H = 1.4 + (α + β/D)-γ (17) 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: D/(H - 1.4)0.4 = α × (D + β) (18) 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 Salekin et al. New Zealand Journal of Forestry Science (2020) 50:2 Page 3 TABLE 1: Summary of the data used for modelling* Variable Unit Statistical summary of variable Mean Min. Max. SD Plot size ha 0.06 0.04 0.10 0.02 Age (t) years 14 3 25 6 Individual tree height (h) m 12.90 0.10 39.80 9.05 Mean top height (MTH) m 18.98 3.50 28.80 7.05 Diameter at breast height at 1.4 m (DBH) cm 22.90 0.10 62.30 14.49 Max DBH (Dmax) cm 39.79 5.40 62.30 13.66 Standard deviation of DBH (SDD) cm 5.34 1.35 11.86 2.20 Basal area (G) m2ha-1 30.59 0.54 77.88 18.84 Volume (V) m3 ha-1 161.34 0.40 538.60 130.09 Stocking (N) stems ha-1 496.99 141.09 1375 317.33 Elevation (Elv) m. asl 211.70 80 300 100.41 Slope (°) 23.27 8 42 13.06 *A detailed individual PSP description including silvicultural treatments is provided in the Appendix. Salekin et al. New Zealand Journal of Forestry Science (2020) 50:2 Page 4 FIGURE 1: Permanent sample plot (PSP) locations and topography. Generic name Expression No. Po ly m or ph ic fo rm Schumacher 1 1 Schumacher 2 2 Gompertz 1 3 Gompertz 2 4 Weibull 1 5 Weibull 2 6 Hossfeld 7 Von Bertalanffy-Richards 1 8 Von Bertalanffy-Richards 2 9 Von Bertalanffy-Richards 3 10 A na m or ph ic fo rm Schumacher A1 11 Schumacher A2 12 Gompertz 13 Von Bertalanffy-Richards 14 Weibull 15 Hossfeld 16 Y2 = e ln(Y1)( t1 t2 )+α(1−t1t2 ) Y2 = e ln(Y1)( t1 t2 ) 𝛾𝛾 +α[1−(t1t2 )𝛾𝛾] 𝑌𝑌2 = 𝑒𝑒ln⁡(𝑌𝑌1)𝑒𝑒 −𝛽𝛽(𝑡𝑡2−𝑡𝑡1)𝑒𝑒𝛼𝛼[1−𝛽𝛽(𝑡𝑡2−𝑡𝑡1)] 𝑌𝑌2 = 𝑒𝑒ln⁡(𝑌𝑌1)𝑒𝑒 −𝛽𝛽(𝑡𝑡2−𝑡𝑡1)+𝛾𝛾(𝑡𝑡2 2−𝑡𝑡1 2)𝑒𝑒𝛼𝛼[1−𝑒𝑒−𝛽𝛽(𝑡𝑡2−𝑡𝑡1)+𝛾𝛾(𝑡𝑡2 2−𝑡𝑡1 2)] 𝑌𝑌2 = 𝑌𝑌1𝑒𝑒−𝛽𝛽(𝑡𝑡2 𝛾𝛾−𝑡𝑡1 𝛾𝛾) + 𝛼𝛼[1 − 𝛽𝛽(𝑡𝑡2 𝛾𝛾 − 𝑡𝑡1 𝛾𝛾)] 𝑌𝑌2 = 𝛼𝛼 − 𝛽𝛽( 𝛼𝛼 − 𝑌𝑌1 𝛽𝛽 ) (𝑡𝑡2𝑡𝑡1 )𝛾𝛾 𝑌𝑌2 = 1 1 𝑌𝑌1 (𝑡𝑡2𝑡𝑡1 )𝛽𝛽 + 𝛼𝛼[1 − (𝑡𝑡2𝑡𝑡1 )𝛽𝛽] 𝑌𝑌2 = 𝛼𝛼( 𝑌𝑌1 𝛼𝛼) ln⁡[1−𝑒𝑒(𝛽𝛽𝑡𝑡2)] ln⁡[1−𝑒𝑒(𝛽𝛽𝑡𝑡1)] 𝑌𝑌2 = 𝛼𝛼{1 − [1 − ( 𝑌𝑌1 𝛼𝛼) 1−𝜗𝜗] 𝑡𝑡2 𝑡𝑡1} 1 1−𝜗𝜗 𝑌𝑌2 = 𝛼𝛼{1 + [( 𝛼𝛼 𝑌𝑌1 ) 𝜗𝜗 − 1]𝑒𝑒[−𝛽𝛽(𝑡𝑡2−𝑡𝑡1)]} 1 𝜗𝜗 𝑌𝑌2 = 𝑌𝑌1𝑒𝑒 −𝛽𝛽( 1𝑡𝑡2 − 1𝑡𝑡1 ) 𝑌𝑌2 = 𝑌𝑌1𝑒𝑒 −𝛽𝛽[( 1𝑡𝑡2 ) 𝛾𝛾 −( 1𝑡𝑡2 )𝛾𝛾] 𝑌𝑌2 = 𝑌𝑌1 𝑒𝑒−𝛽𝛽𝑒𝑒−𝛾𝛾𝑡𝑡2 𝑒𝑒−𝛽𝛽𝑒𝑒−𝛾𝛾𝑡𝑡1 𝑌𝑌2 = 𝑌𝑌1[ 1 − 𝑒𝑒−𝛽𝛽𝑡𝑡2 1 − 𝑒𝑒−𝛽𝛽𝑡𝑡1] 𝛾𝛾 𝑌𝑌2 = 𝑌𝑌1 1 − 𝑒𝑒−𝛽𝛽𝑡𝑡2𝛾𝛾 1 − 𝑒𝑒−𝛽𝛽𝑡𝑡1𝛾𝛾 𝑌𝑌2 = 1 1 𝑌𝑌1 + 𝛽𝛽(1𝑡𝑡2 )𝛾𝛾 − (1𝑡𝑡1 )𝛾𝛾 TABLE 2. Different forms of differential equations. 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 (R2) (Equation 27). Adjusted R2 values were not considered for assessing differential equations as it is sensitive to grouped and repeated data (Warren 1971). (23) (24) (25) (26) (27) where N = Number of observations, O = Observed value, = mean of observed values, P = Predicted value, . K denotes the number of estimated parameters. The predictive ability of the models was evaluated using prediction errors or predictive residual error sum square (PRESS) statistics (Equation 28), Oi - P(i,-i) = e(i,-i) (i = 1,2,…..,n) (28) where Oi 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): (29) 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. R2”, 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). Results 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 Salekin et al. New Zealand Journal of Forestry Science (2020) 50:2 Page 5 Expression Reference No. (Soalleiro 1995) 19 (Jansen et al. 1996) 20 (Burkhart 1977) 21 (Candy 1989) 22 , and α, β, γ, and δ are model coefficients. V = α × G × MTH V = G × MTH(α+βt)e(γ+δt) V = G × (α + β MTH) V = e(α+βlogMTH)+γlogG TABLE 3. Volume yield equations. Here, V is volume/ha, G is basal area/ha, MTH is mean top height, t is age in years, and α, β, γ, and δ are model coefficients. RMSE = √∑ (Pi−Oi) 2N i=1 N MAE = ∑ |Pi−Oi| N i=1 N BIAS = ∑ (Pi− N i=1 Oi ) N R2 = ∑ P ′ i 2 ∑ O′ i 2 𝑅𝑅2 𝑎𝑎𝑎𝑎𝑎𝑎 = 1 − [ (1−𝑅𝑅2)(𝑁𝑁−1) 𝑁𝑁−𝐾𝐾−1 ] O O̅ P P′i = Pi − O̅ O′i = Oi − O̅ K Oi − Pi,−i = ei,−i (i = 1,2, … . . , n) RMSE = √∑ (Pi−Oi) 2N i=1 N MAE = ∑ |Pi−Oi| N i=1 N BIAS = ∑ (Pi− N i=1 Oi ) N R2 = ∑ P ′ i 2 ∑ O′ i 2 𝑅𝑅2 𝑎𝑎𝑎𝑎𝑎𝑎 = 1 − [ (1−𝑅𝑅2)(𝑁𝑁−1) 𝑁𝑁−𝐾𝐾−1 ] O O̅ P P′i = Pi − O̅ O′i = Oi − O̅ K Oi − Pi,−i = ei,−i (i = 1,2, … . . , n) O̅ P P′i = Pi − O̅ O′i = Oi − O̅ PRESS = ∑ Oi − (Pi,−i)2ni=1 = ∑ (ei,−i)2ni=1 O̅ P P′i = Pi − O̅ O′i = Oi − O̅ PRESS = ∑ Oi − (Pi,−i)2ni=1 = ∑ (ei,−i)2ni=1 O̅ P P′i = Pi − O̅ O′i = Oi − O̅ PRESS = ∑ Oi − (Pi,−i)2ni=1 = ∑ (ei,−i)2ni=1 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 (Dmax) model The Hossfeld polymorphic model (Equation 7) predicted the maximum diameter (Dmax) 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 goodness- of-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 Dmax enveloped all the measurements and followed a sigmoid shape (Figure 4C), which ensures Salekin et al. New Zealand Journal of Forestry Science (2020) 50:2 Page 6 Action RMSE MAE BIAS SE MPRESS MAPRESS Fitting 7.185 5.467 -1.777 1.116 - - Validation 3.852 2.512 0.066 1.112 0.009 0.946 TABLE 4. Mean top height (MTH) m, Von Bertalanffy-Richards polymorphic model fitting and validation statistics. FIGURE 2. Mean top height (MTH) model results: A) Residuals against prediction plot of first Von Bertalanffy-Richards polymorphic equation, light blue points represent model fitting, red points indicate validation residuals, and model fit is shown by the black line; B) Residuals frequency distribution, red dashed line shows the mean; and C) Model fit (blue lines) over measured MTH (thin black lines). biological rationality. Model parameters are provided in the Appendix (Table A2). Standard deviation of diameter (SDD) model Among all the candidate models tested (Table 2), the standard deviation of diameter (SDD) 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 Salekin et al. New Zealand Journal of Forestry Science (2020) 50:2 Page 7 distribution plot (Figure 5B). The prediction plot shows that the model included the full range of measured SDD (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 R2 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). Action RMSE MAE BIAS SE MPRESS MAPRESS Fitting 25.303 21.250 2.893 6.893 - - Validation 13.431 9.988 0.653 6.800 1.054 0.841 TABLE 5. Basal area/ha (G), first anamorphic Schumacher model fitting and validation statistics. FIGURE 3. Basal area (G) model results: A) Residuals against prediction plot of first Schumacher anamorphic equation, light blue points represent model fitting, the red points indicate validation residuals, and model fit is shown by the black line; B) Residuals frequency distribution, red dashed line shows the mean; and C) Model fit (blue lines) over measured G (thin black lines). 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 TABLE 2: Confusion matrix 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 Salekin et al. New Zealand Journal of Forestry Science (2020) 50:2 Page 8 Action RMSE MAE BIAS SE MPRESS MAPRESS Fitting 2.400 1.759 0.054 2.411 - - Validation 6.699 4.681 -0.061 2.388 0.059 0.932 TABLE 6. Maximum diameter cm (Dmax) model fitting and validation statistics. FIGURE 4. Maximum diameter (Dmax) model results: A) Residuals against prediction plot of Hossfeld polymorphic equation, light blue points represent model fitting, the red points indicate validation residuals, and model fit is shown by the black line; B) Residuals frequency distribution, red dashed line shows the mean; and C) Model fit (blue lines) over measured Dmax (thin black lines). 𝐻𝐻 = 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, Dmax and SDD 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 m2 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 Salekin et al. New Zealand Journal of Forestry Science (2020) 50:2 Page 9 Action RMSE MAE BIAS SE MPRESS MAPRESS Fitting 1.571 1.224 0.412 1.577 - Validation 1.959 1.513 0.337 1.569 0.407 0.596 TABLE 7. Standard deviation of DBH cm (SDD) model fitting and validation statistics. FIGURE 5. Standard deviation of DBH (SDD) model results: A) Residuals against prediction plot of Hossfeld polymorphic equation, light blue points represent model fitting, the red points indicate validation residuals, and model fit is shown by the black line; B) Residuals frequency distribution, red dashed line show the mean; and C) Model fit (blue lines) over measured SDD (thin black lines). 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 competition- mortality” 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 Salekin et al. New Zealand Journal of Forestry Science (2020) 50:2 Page 10 Action RMSE MAE BIAS SE R2adj MPRESS MAPRESS Fitting 39.122 27.983 -1.102 40.5 0.91 - - Validation 140.959 89.827 -0.582 39.413 0.95 -0.845 0.868 TABLE 8. Stand volume m3 (V), Jansen et al. (1996) model fitting and validation statistics. FIGURE 6. Stand volume (V) model results: A) Estimated stand volume from measured data; B) Residuals against prediction plot, light blue points represent model fitting, red points indicate validation residuals, and model fit is shown by the black line; and C) Residuals frequency distribution, red dashed line is shown the mean. Salekin et al. New Zealand Journal of Forestry Science (2020) 50:2 Page 11 FIGURE 7. A) Measured height-diameter (H-D), blue line shows the linear trend; B) Residuals against prediction plot, light blue points represent model fitting, red points indicate validation residuals, model fit is shown by the blue line; and C) Residuals frequency distribution, red dashed line is shown the mean. Action RMSE MAE BIAS SE R2adj MPRESS MAPRESS Fitting 3.08 2.41 -0.01 3.10 0.54 - - Validation 3.80 2.91 0.03 3.05 0.53 -0.001 0.530 TABLE 9. Height-diameter relationship (H-D) model fitting and validation statistics. FIGURE 8. A) Reineke’s SDI curve represented with self-thinning lines and B) SDI distribution plot. tree size (55% relative density). According to this theory, a uniform stand may have imminent competition- mortality 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 m2 ha-1 and volume ha-1 ranging from 52–252 m3. 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 long- term plot measurement data is indispensable to make comprehensive forecasting models needed to underpin management decisions. Additional files Additional file 1: Appendix 1. List of abbreviations ADA: Algebric difference approach; DBH: Diameter at breast height (1.4m); MTH: Mean top height; G: Basal area; SDD: Standard deviation of diameter; Dmax: Maximum diameter; V: Stand volume; SDI: Stand density index; RMSE: Root mean square error; SE: Standard error; PRESS: Predictive residual error sum square. Software A web-based E. globoidea growth and yield software is available from the following link: http://www. treesandstars.com/models/egloboidea.htm Competing interests The authors declare that they have no competing interests. Authors’ contributions SS and EGM conceived the study from the data provided by DFM. SS analysed the data and wrote the manuscript under supervision of EGM and JM. All the authors read and approved the final version of the manuscript. Funding Funding was received from Agricultural and Marketing Research and Development trust (AgMARDT) and Speciality Wood Product programme (SWP), New Zealand. Acknowledgements The authors are grateful to AgMARDT and SWP for providing financial support. They would like to thank Scion Ltd. for making the data available and Professor Horacio Bown for providing comments on the initial draft. They also appreciate the suggestions and comments provided by two anonymous reviewers and co-editor Dr John Moore on earlier version of this manuscript. References Apiolaza, L., McConnochie, R., & Millen, P. (2011). Introducing durable species to New Zealand drylands: genetics of early adaptation of Eucalyptus bosistoana. Paper presented at the Developing a Eucalypt Resource: Learning from Salekin et al. New Zealand Journal of Forestry Science (2020) 50:2 Page 12 Australia and Elsewhere symposium, Wood Technology Research Centre, Christchurch, New Zealand, 2 November 2011. https://luis.apiolaza. net/uploads/Apiolaza_et_al_2011_Eucalyptus_ bosistoana_in_NZ_drylands.pdf Arlot, S., & Celisse, A. (2010). A survey of cross-validation procedures for model selection. Statistics Surveys, 4, 40-79. https://doi.org/10.1214/09-SS054 Bailey, R.L., & Clutter, J.L. (1974). Base-age invariant polymorphic site curves. Forest Science, 20(2), 155- 159. Battaglia, M., Sands, P., White, D., & Mummery, D. (2004). CABALA: a linked carbon, water and nitrogen model of forest growth for silvicultural decision support. Forest Ecology and Management, 193(1), 251-282. https://doi.org/10.1016/j.foreco.2004.01.033 Bell, D.T., & Williams, J.E. (1997). Eucalypt ecophysiology. In Williams, J., & Woinarski, J. (Eds.). Eucalypt Ecology: Individuals to Ecosystems. Cambridge, UK: Cambridge University Press, pp. 168-196. Belli, K.L., & Ek, A.R. (1988). Growth and survival modeling for planted conifers in the great lakes region. Forest Science, 34(2), 458-473. Berrill, J., Nicholas, I.D., & Gifford, H.H. (2007). Preliminary growth and yield models for even- aged Acacia melanoxylon plantations in New Zealand. New Zealand Journal of Forestry Science, 37(1), 37-56. Berrill, J.P. (2004). Preliminary growth and yield models for even-aged Cupressus lusitanica and C. macrocarpa plantations in New Zealand. New Zealand Journal of Forestry Science, 34(3), 272-292. Berrill, J.P., & Hay, A.E. (2005). Indicative growth and yield models for even-aged Eucalyptus fastigata plantations in New Zealand. New Zealand Journal of Forestry Science, 35(2/3), 121-138. Berrill, J.P., & Hay, A.E. (2006). Indicative growth and yield models for stringybark eucalypt plantations in northern New Zealand. New Zealand Journal of Forestry, 51(1), 19-22. Blake, J., Somers, G., & Ruark, G. (1990). Perspectives on process modeling of forest growth responses to environmental stress. In: R. Dixon, R. Meldahl, G. Ruark, & W. Warren (Eds.), Process Modeling of Forest Growth Responses to Environmental Stress (pp. 9-17). Portland, OR, USA: Timber Press. Blanco, J.A., Zavala, M.A., Imbert, J.B., & Castillo, F.J. (2005). Sustainability of forest management practices: Evaluation through a simulation model of nutrient cycling. Forest Ecology and Management, 213(1), 209-228. https://doi.org/10.1016/j. foreco.2005.03.042 Borders, B.E. (1989). Systems of equations in forest stand modeling. Forest Science, 35(2), 548-556. Borders, B.E., Bailey, R.L., & Clutter, M.L. (1988). Forest Growth Models: Parameter Estimation Using Real Growth Series. USDA Forest Service General Technical Report. Minneapolis, MI, USA:NC-North Central Forest Experiment Station. Box, G.E. (1976). Science and statistics. Journal of the American Statistical Association, 71(356), 791-799. https://doi.org/10.1080/01621459.1976.104809 49 Castedo-Dorado, F., Diéguez-Aranda, U., Barrio-Anta, M., & Álvarez-Gonzàlez, J.G. (2007). Modelling stand basal area growth for radiata pine plantations in Northwestern Spain using the GADA. Annals of Forest Science, 64(6), 609-619. https://doi. org/10.1051/forest:2007039 Castedo-Dorado., F., Diéguez-Aranda, U., & Álvarez- González, J.G. (2007). A growth model for Pinus radiata D. Don stands in north-western Spain. Annals of Forest Science, 64(4), 453-465. https:// doi.org/10.1051/forest:2007023 Chou, C.K.S. (1991). Perspectives of disease threat in large-scale Pinus radiata monoculture - the New Zealand experience. European Journal of Forest Pathology, 21(2), 71-81. https://doi. org/10.1111/j.1439-0329.1991.tb00946.x Clutter, J.L. (1963). Compatible growth and yield models for Loblolly Pine. Forest Science, 9(3), 354-371. Clutter, J.L., Fortson, J.C., Pienaar, L.V., Brister, G.H., & Bailey, R.L. (1983). Timber management: a quantitative approach. New York, USA: John Wiley & Sons, Inc. Curtis, R.O. (1967). Height-diameter and height- diameter-age equations for second growth Douglas Fir. Forest Science, 13, 365-375. Drew, T.J., & Flewelling, J.W. (1977). Some recent Japanese theories of yield-density relationships and their application to Monterey pine plantations. Forest Science, 23(4), 517-534. Drew, T.J., & Flewelling, J.W. (1979). Stand density management: an alternative approach and its application to Douglas-fir plantations. Forest Science, 25(3), 518-532. Ek, A.R. (1974). Nonlinear models for stand table projection in northern hardwood stands. Canadian Journal of Forest Research, 4(1), 23-27. https://doi. org/10.1139/x74-004 Garcia, O. (1974). Height-diameter equations for Pinus radiata. [Nota Tecnica, No. 19] Instituto Forestal, Chile. Garcia, O. (1988). Growth modelling — A review. New Zealand Forestry, 33(3), 14-17. García, O. (1981). Simplified method-of-moments estimation for the Weibull distribution. New Zealand Journal of Forestry Science, 11(3), 304-306. Goulding, C.J. (2005). Measurement of trees. In M. Colley (Ed.), Forestry Handbook (4th ed., p. 316). Christchurch, New Zealand: The New Zealand Insitute of Forestry (Inc.). Hamner, B., & Frasco, M. (2018). Metrics: Evaluation Metrics for Machine Learning. R package version 0.1.4. https://CRAN.R-project.org/ package=Metrics Jansen, J.J., Sevenster, J., & Faber, J. (1996). Opbrengsttabellen voor belangrijke boomsoorten in Nederland. Wageningen, The Netherlands: IBN- DLO Kakitani, T. (2017). The global timberlization movement and the potential for durable eucalyts: downstream Salekin et al. New Zealand Journal of Forestry Science (2020) 50:2 Page 13 opportunities. Paper presented at the Durable Eucalypts on Drylands: Protecting and Enhancing Value symposium, Marlborough Research Centre, Blenheim, New Zealand 19 April 2017. Kitikidou, K., Papageorgiou, A., Milios, E., & Stampoulidis, A. (2016). Preliminary individual tree growth model, site index model, mortality model for Aleppo Pine (Pinus halepensis MILL.) in chalkidiki (Northern Greece). Journal Ciência Florestal, 26(4), 1247- 1257. https://doi.org/10.5902/1980509825136 Korzukhin, M.D., Ter-Mikaelian, M.T., & Wagner, R.G. (1996). Process versus empirical models: which approach for forest ecosystem management? Canadian Journal of Forest Research, 26(5), 879- 887. https://doi.org/10.1139/x26-096 Kuru, G.A., Whyte, A.G.D., & Woollons, R.C. (1992). Utility of reverse Weibull and extreme value density functions to refine diameter distribution growth estimates. Forest Ecology and Management, 48(1), 165-174. https://doi.org/10.1016/0378- 1127(92)90128-V Landsberg, J.J., Waring, R.H., & Coops, N. C. (2003). Performance of the forest productivity model 3-PG applied to a wide range of forest types. Forest Ecology and Management, 172(2), 199-214. https://doi.org/10.1016/S0378-1127(01)00804- 0 Lee, S.H. (1998). Modelling growth and yield of Douglas- fir using different interval lengths in the South Island of New Zealand. (PhD), University of Canterbury, Christchurch, New Zealand. Lin, H. (2017). Risk and impact of insect herbivores on the development of dryland Eucalyptus forestry in New Zealand. (PhD), University of Canterbury, Christchurch, New Zealand. Long, J.N. (1985). A practical approach to density management. The Forestry Chronicle, 61(1), 23-27. https://doi.org/10.5558/tfc61023-1 Lundgren, A.C., Gordon, A. & Hay, E. (1995). Taper equations for Eucalyptus pilularis, E. globoidea and E. muelleriana in New Zealand. Rotorua, New Zealand: Management of Eucalypts Cooperative. Mason, E.G. (2019). Influences of mean top height definition and sampling method on errors of estimates in New Zealand’s forest plantations. New Zealand Journal of Forestry Science, 49:1. https:// doi.org/10.33494/nzjfs492019x24x Meason, D., Bellè, P., Höck, B., Lumnitz, S., & Todoroki, C. (2016). Alternative species site mapping review and analysis. Retrieved 20 July 2016 from New Zealand Farm Forestry Association https://www.nzffa.org. nz/system/assets/5143/SWP-T004.pdf Menzies, H. (1995). Eucalypts show potential. Farm Forestry Review, 33-34. Methol, R.J. (2001). Comparisons of approaches to modelling tree taper, stand structure and stand dynamics in forest plantations. (PhD), University of Canterbury, Christchurch, New Zealand. Millen, P., van Ballekom, S., Altaner, C., Apiolaza, L., Mason, E., McConnochie, R.,Murray, T. (2018). Durable eucalypt forests–a multi-regional opportunity for investment in New Zealand drylands. New Zealand Journal of Forestry, 63(1), 11-23. Myers, R.H., & Myers, R.H. (1990). Classical and modern regression with applications (Vol. 2): Belmont, CA, USA: Duxbury Press. Näslund, M. (1936). Skogsförsöksanstaltens gallringsförsök i tallskog. Meddelanden från statens skogsförsöksanstalt (Reports of the Swedish Institute of Experimental Forestry) No. 29. Stockholm, Sweden: Centraltrvckeriet, Esselte h t t p s : / / p u b . e p s i l o n . s l u . s e / 1 0 1 5 9 / 1 / m e d d _ statens_skogsforskningsanst_029_01.pdf Ministry for Primary Industries. (2019). National Exotic Forest Description. Wellington, New Zealand. Nicholas, I., & Millen, P. (2012a). Eucalyptus globoidea. [Durable Eucalypt leaflet series]. Blenheim, New Zealand: New Zealand Dryland Forest Initiative. https://nzdfi.org.nz/wpcontent/ uploads/2014/08/Leaflet3E.globoideaDec2012. pdf Nicholas, I., & Millen, P. (2012b). Eucalyptus globoidea - a durable hardwood for planting in the Bay of Plenty region. Blenheim, New Zealand: Marlborough Research Centre. New Zealand Dryland Forest Initiative. (2013). Introductory Brochure. Blenheim, New Zealand: Marlborough Research Centre. https://nzdfi.org. nz/wp-content/uploads/2015/04/Eucalyptus- globoidea-a-durable-hardwood-for-planting-in- the-BoP-region.pdf New Zealand Forest Owners Association. (2017). New Zealand Plantation Forest Industry: Facts and Figures 2016/17. Wellington, New Zealand. https:// www.nzfoa.org.nz/images/stories/pdfs/Facts_ Figures_2016_%C6%92a_web_version_v3.pdf Page, D., & Sing, T. (2014). Durability of New Zealand grown timbers. New Zealand Journal of Forestry, 58(4), 26-30. Palahí, M., & Grau, J. (2003). Preliminary site index model and individual-tree growth and mortality models for black pine (Pinus nigra Arn.) in Catalonia (Spain). Journal of Forest Systems, 12(1), 137-148. Park, A., & Wilson, E.R. (2007). Beautiful plantations: can intensive silviculture help Canada to fulfil ecological and timber production objectives? The Forestry Chronicle, 83(6), 825-839. https://doi. org/10.5558/tfc83825-6 Pienaar, L.V., & Rheney, J.W. (1995). Modeling stand level growth and yield response to silvicultural treatments. Forest Science, 41(3), 629-638. Pilaar, C.H., & Dunlop, J.D. (1990). The permanent sample plot system of the New Zealand Ministry of Forestry. Bulletin des Recherches Agronomiques de Gembloux, 25(1), 5-17. Pretzsch, H., & Biber, P. (2005). A re-evaluation of Reineke’s rule and stand density index. Forest Science, 51(4), 304-320. R Development Core Team. (2017). R: A language and environment for statistical computing. Vienna, Austria: R Foundation for Statistical Computing; 2016. http://www. R-project. org. Salekin et al. New Zealand Journal of Forestry Science (2020) 50:2 Page 14 Reineke, L.H. (1933). Perfecting a stand-density index for even-aged forest. Journal of Agricultural Research, 4, 627-638. Sánchez-González, M., Tomé, M., & Montero, G. (2005). Modelling height and diameter growth of dominant cork oak trees in Spain. Annals of Forest Science, 62(7), 633-643. https://doi.org/10.1051/ forest:2005065 Saunders, M.R., & Puettmann, K.J. (2000). A preliminary white spruce density management diagram for the Lake States. [Staff Paper Series No 145]. St Paul, MI, USA: Department of Forest Resources, University of Minnesota. Tewari, V.P., & von Gadow, K. (2003). Preliminary growth models for Prosopis cineraria (L.) Druce plantations in the hot arid region of India. Forests, Trees and Livelihoods, 13(4), 361-373. https://doi.org/10.10 80/14728028.2003.9752471 Turner, J.A., West, G., Dungey, H., Wakelin, S., Maclaren, P., Adams, T., & Silcock, P.J. (2008). Managing New Zealand planted forests for carbon: A review of selected scenarios identification of knowledge gaps. Wellington, New Zealand: Ministry of Agriculture Forestry. Vanclay, J.K. (1994). Modeling forest growth and yield. Applications to mixed tropical forests. Wallingford, UK: CAB International. Vanclay, J.K. (2010). Robust relationships for simple plantation growth models based on sparse data. Forest Ecology and Management, 259(5), 1050-1054. https://doi.org/10.1016/j. foreco.2009.12.026 Warren, W.G. (1971). Correlation or regression: bias or precision. Journal of the Royal Statistical Society. Series C (Applied Statistics), 20(2), 148-164. https://doi.org/10.2307/2346463 Wickham, H. (2016). ggplot2: elegant graphics for data analysis. New York, USA: Springer-Verlag. https:// doi.org/10.1007/978-3-319-24277-4 Williams, J., & Woinarski, J. (Eds.). (1997). Eucalypt Ecology: Individuals to Ecosystems. Cambridge, UK: Cambridge University Press. Woollons, R.C., Snowdon, P., & Mitchell, N.D. (1997). Augmenting empirical stand projection equations with edaphic and climatic variables. Forest Ecology and Management, 98(3), 267-275. https://doi. org/10.1016/S0378-1127(97)00090-X Zeide, B. (1993). Analysis of growth equations. Forest Science, 39(3), 594-616. https://doi.org/10.1093/ forestscience/39.3.594 Zhao, W. (1999). Growth and yield of Pinus radiata in Canterbury, New Zealand. (PhD), University of Canterbury, Christchurch. (SD397.P6117.Z63) Zhao, W., Mason, E.G., & Brown, J. (2006). Modelling height-diameter relationships of Pinus radiata plantations in Canterbury, New Zealand. New Zealand Journal of Forestry, 51(1), 23-27. Salekin et al. New Zealand Journal of Forestry Science (2020) 50:2 Page 15 New Zealand Journal of Forestry Science Additional File 1: Appendix A preliminary growth and yield model for Eucalyptus globoidea (Blakely) plantations in New Zealand VAlg = π ∙ DBH2 ∙ H3 40000 ∙ (H − 1.4)2 ∙ (0.2134788 ∙ βc + 0.011344) 𝛽𝛽𝑐𝑐 βc = (1 − 1.4 H ) 0.280729 − 0.414429 (1 − 1.4 H ) 17.26155 Eucalyptus globoidea Age(years) Residual Stems ha-1 Age(years) Height (m) New Zealand Journal of Forestry Science . Preliminary model’s parameter estimates. 𝛼𝛼 β 𝛾𝛾 δ MTH p Sig. G p Sig. Dmax p Sig. SDD p Sig. V p Sig. α0 α1 β0 β1 H-D p Sig.