Modelling the E ﬀ ect of Microsite Inﬂuences on the Growth and Survival of Juvenile Eucalyptus globoidea (Blakely) and Eucalyptus bosistoana (F. Muell) in New Zealand

: The e ﬀ ect of microsite on juvenile forest plantation yield is rarely explored. This is because juvenile plantation growth is considered to be reasonably homogenous due to a lack of resource competition between trees prior to canopy closure. However, models of juvenile plantation height growth and survival that are sensitive to microsite variation could aid decisions relating to site preparation, plantation establishment and early silvicultural treatments. In this study, juvenile Eucalyptus bosistoana and E. globoidea height growth and survival proportion were modelled against topographic and environmental microsite characteristics as independent variables. The experiment included three di ﬀ erent sites situated in a sub-humid region of New Zealand. A total of 540 plots were planted with 18,540 trees in regular rows and columns. Micro-topographical variables signiﬁcantly inﬂuenced height growth and survival proportion of both E. bosistoana and E. globoidea, but species di ﬀ ered in their responses. More sheltered microsites yielded greater height growth and survival for both species. The height of both species was inﬂuenced by wind exposure, morphometric protection, and distance from the nearest ridge. E. bosistoana height was also inﬂuenced by topographic position and surface plan curvature. Survival was a ﬀ ected by surface proﬁle curvature for both species, while E. globoidea survival was also impacted by surface plan curvature and distance from the top ridge. This study identiﬁed microsite factors inﬂuencing juvenile height and survival of two Eucalyptus species.


Introduction
The term "site", used as a primary ecological unit, plays an important role as one of the principal factors in the survival and growth of trees at different scales [1]. It refers to a geographic location with a homogenous physical and biological environment [2,3]. In a forestry context, plantation forest sites, typically called stands, are bounded areas that receive similar silvicultural treatments [4,5]. Although plantation forests are homogenised through silviculture, their growth shows considerable spatial and temporal variability [6,7]. Skovsgaard and Vanclay [5] defined site productivity as the potential of a particular stand to produce aboveground biomass. Variation in site productivity has long been a subject of interest to growth and survival. This information is crucial for identifying optimal sites and positions within a site (e.g. hillslope) to successfully establish these species, especially at the large scale.
This study explores a comprehensive set of topographic, edaphic and climatic explanatory variable effects at the microsite level on the growth and survival of small plots of trees in a juvenile plantation. The main research objectives were to identify microsite level topographic and climatic variables that influenced the height growth and survival of juvenile E. globoidea and E. bosistoana, and to include these into a height growth and survival model, respectively.

Experimental Sites
The three experimental sites for this study were chosen based on (i) their locations, (ii) ease of access for carrying instruments and (iii) landowners' willingness to host the trials. The study sites were situated in a sub-humid zone close to Blenheim, New Zealand ( Figure 1). Sites A, B and C, so named to maintain landowner privacy, have areas 4.7, 3.7 and 2.2 hectares, respectively, and are planted with E. globoidea (Site A) and E. bosistoana (Sites B and C) ( Table 1).
The region in which the trial sites are located is sheltered by high country to the west, south and in some areas to the east, and it is one of the sunniest regions of New Zealand with 2487 mean annual hours of sunshine [56]. Warm, dry and settled weather predominates during summer, while winter days often begin with a frost, but are usually mild overall. Typical summer daytime maximum air temperatures range from 20 • C to 26 • C but occasionally rise above 30 • C. Typical winter daytime maximum air temperatures range from 10 • C to 15 • C [56]. South-westerly winds prevail in Blenheim, though dry Foehn winds from the northwest also occur [56]. The soils at these sites are formed from loess and classified as Pallic Argillic soils, which have thin subsoil bands of accumulated clay [57] ( Table 2).       (Table 1). Trees were planted in regular rows and columns within plots, with spacing equal to 2.4 m × 1.8 m in all sites. Weeds were controlled with herbicide in 1 m diameter spots around each tree immediately after establishment.
There were approximately 18,500 trees at the three sites with varying stand density ( Table 1). The height (h), diameter at breast height 1.4 m (DBH), and tree status (dead or alive) were measured for all trees. All tree measurements were undertaken during November 2015 and January 2016 and again in June-August 2017 (Table 1). Prior to these measurements, the height and tree status were measured at age 1.2 years for all trees. Individual tree height and survival data were averaged at each plot. The survival proportion (S) was calculated for each plot from the number of surviving trees.
Height data from all three sites were used to create juvenile height models. For survival data, only the A and C sites' survival proportions (S) were used to create juvenile survival models, as there was a thinning trial at site B prior to completion of field measurements for this study.

Topographic Data
Digital elevation models (DEMs) for all sites were produced by using a real-time kinetic geo-positioning system (RTK-GPS). The unit was carried on transect lines across the sites, with coordinates and elevation collected at five-metre intervals along the transects. The elevation points were interpolated into a DEM with five metre spatial resolution through topo-raster algorithm (ANUDEM), using a process detailed in Salekin, et al. [43].
Next, primary and secondary surface attributes were derived from the DEM for each plot (~15 m spatial resolution). The primary surface attributes included aspect and slope [59]. From these, the following secondary surfaces were derived: Total (CURV), profile (CVPRO) and plan curvature (CVPLA) [60,61]; topographic ruggedness index (TRI) [62]; topographic position index (TPI) [63]; topographic wetness index (TWI) [64,65]; wind exposure index (WEI) [66]; morphometric protection index (MPI) [67] and Euclidian distance between each plot's centre pixel and the nearest ridge line (DIST) (detailed descriptions of these indices are provided as supplementary material, S1). All surfaces were interpolated or derived using ArcMap 10.4 [68] or the System for Automated Geoscientific Analysis (SAGA) [69]. A summary of the topographic attributes is presented in Table 3. Table 3. Summary of the topographic attributes for study sites. TPI-topographic ruggedness index; TRI-topographic position index; TWI-topographic wetness index; WEI-wind exposure index; MPI-morphometric protection index; DIST-Euclidian distance between each plot's centre pixel and the nearest ridge line.

Climatic Data
Each site had a meteorological station, which was equipped with solar radiation, air temperature, soil moisture, wind and rain sensors. There were 20 additional air temperature sensors installed at sites A (n = 10) and B (n = 10), one meter above ground to measure the air temperature variation across and within the sites. All the sensors, including those within the meteorological stations, logged data at 30-minute intervals from March 2015 to April 2017.
Average and maximum daily air temperatures were calculated on a monthly basis for the duration of the measurement period. These calculations were based on data collected by the 10 sensors at each of sites A and B ( Figure 2). Air temperature differences between these sensors and the single air temperature sensor within each meteorological station were calculated ( Figure 2C,D). This was done by deducting the independent sensors' temperatures from the meteorological station's above-ground temperature reading. From the calculated air temperature difference and sensors' positions, a mixed-effect model was developed by linking primary topographic attributes. This model was used to simulate air temperature differences and estimate the air temperature within each microsite. Detail model development procedures and final models are described in supplementary material, S2.

Climatic Data
Each site had a meteorological station, which was equipped with solar radiation, air temperature, soil moisture, wind and rain sensors. There were 20 additional air temperature sensors installed at sites A (n = 10) and B (n = 10), one meter above ground to measure the air temperature variation across and within the sites. All the sensors, including those within the meteorological stations, logged data at 30-minute intervals from March 2015 to April 2017.
Average and maximum daily air temperatures were calculated on a monthly basis for the duration of the measurement period. These calculations were based on data collected by the 10 sensors at each of sites A and B ( Figure 2). Air temperature differences between these sensors and the single air temperature sensor within each meteorological station were calculated ( Figures 2C and D). This was done by deducting the independent sensors' temperatures from the meteorological station's above-ground temperature reading. From the calculated air temperature difference and sensors' positions, a mixed-effect model was developed by linking primary topographic attributes. This model was used to simulate air temperature differences and estimate the air temperature within each microsite. Detail model development procedures and final models are described in supplementary material, S2.

Modelling Strategy
In young plantations, prior to canopy closure, growth is modelled exponentially, with larger trees having greater leaf and root surface areas than smaller trees. For the measured height for Pinus radiata seedlings planted in plantations in New Zealand Mason and Whyte [32] expressed this growth function as,

Modelling Strategy
In young plantations, prior to canopy closure, growth is modelled exponentially, with larger trees having greater leaf and root surface areas than smaller trees. For the measured height for Pinus radiata seedlings planted in plantations in New Zealand Mason and Whyte [32] expressed this growth function as, Here, h 0 = mean height immediately after planting (in this case, 0.25 m), and h T = mean height at stand age T.
Equation (1) has been widely used for modelling juvenile crops [32,70]. Furthermore, Mason and Whyte [32] showed that the coefficients of Equation (1) could be extended as a linear function (Equations (2) and (3)) to independent variables and their interactions by inserting them into linear functions.
The mortality of trees in young plantations is not due to competition among them, but rather water stress or other site-specific factors. According to Mason and Whyte [32] juvenile mortality should be considered as a random process over time with as varying parameter and, therefore, should follow a Weibull probability distribution. The functional form should be anamorphic, as the percentage of deaths would be independent of stocking.
The survival function used by Belli and Ek [70] was one of exponential decay, which was converted to mortality by taking the same Weibull probability density function derivatives given by Mason [71]. Other modellers have used similar approaches [25,70,72]. In this case, the survival proportion function fitted a yield form described in Mason and Whyte [32] (Equation (4)).
where, S T = survival at stand age T, and α and β represent model coefficients.
It is expected that the coefficients should vary with independent explanatory variables, which can be extended linearly by following the same approach as for the height models (Equations (2) and (3)).

Model Testing and Validation
A mixed approach [73,74] was applied to evaluate the model, by performing a full set of residual analyses. Validation included a visual analysis of graphs of the residuals, the calculation of standard error (SE) (Equation (5)), root mean square error (RMSE) (Equation (6)), mean absolute error (MAE) (Equation (7)), and bias (Equation (8)).
where S = standard deviation of the mean, N = number of observations, O = observed value, P = predicted value. There are many established procedures to perform model validation [75]. Should independent datasets not be available, splitting data sets is a commonly accepted practice for model testing and validation, assuming the dataset is sufficiently large [73]. Dobbin, et al. [76] suggested a data splitting ratio of 75:25 (model fitting: validation), which was applied in this study.

Statistical Analysis
All statistical analyses were performed in the R statistical environment [77]. An assessment for multi-collinearity was performed for all explanatory variables by using the variation inflation factor (VIF) with the "vif.mer" function of the car package in R [78]. Elevation, slope, topographic ruggedness index (TRI), and total curvature were highly correlated with other potential independent variables. So by following Cook, et al. [79], the R procedure "anova" for comparing models was employed to determine whether or not these variables added any statistically significant extra information and they were excluded from the model building procedure. Following multi-collinearity analysis, model coefficients were fitted against the explanatory variables by using the "lm" function to predict α and β coefficients for Equations (1) and (4) fitted to data from individual plots.
Finally, height and survival models were fitted using the "nls" function with only significant independent variables. The height and survival models were validated against the validation datasets by using "rmse","mae" and "bias" functions from the metrics and qpcR packages [80,81].

Juvenile Height Models
Juvenile height models for both species (Equations (9)-(11)) had low and stable error statistics ( Table 4). The residuals of the models satisfied all required statistical assumptions (Supplementary material, S3), they were normally distributed and had minimum bias and heteroscedasticity. The E. globoidea model over-predicted height. With the exception of bias, all calculated statistics were lower for the fitting dataset than for the validation dataset (see Table 4).  (Table I). The E. bosistoana height model behaved differently at different sites. At site B, the model standard error was higher than that of site C (Table 4). At site B, RMSE, MAE and SE increased respectively to 0.603, 0.429 and 0.615 from the fit statistics, while BIAS reversed in turn to 0.024 from the fit statistics. In contrast to that, at site C all the fitting statistics features were reduced during validation (Table 4).

Key Variables for Microsite Height Growth
Juvenile E. globoidea height at Site A was significantly correlated with WEI, MPI and plot distance from the top ridge (DIST) ( Table 5 and Figure 3). Therefore, these variables were added to the final height growth model represented by Equation (9). The microsites exposed to high levels of wind had the least height growth, and tree height growth decreased with reduced morphometric protection (MPI). Trees close to the top ridge had the least height growth, while height growth increased proportionally with distance from the ridge until the age of 4.5 years. Beyond that age, trees at the mid-distance from the top ridge grew taller.  Figure 3). Therefore, these variables were added to the final height growth model represented by Equation (9). The microsites exposed to high levels of wind had the least height growth, and tree height growth decreased with reduced morphometric protection (MPI). Trees close to the top ridge had the least height growth, while height growth increased proportionally with distance from the ridge until the age of 4.5 years. Beyond that age, trees at the mid-distance from the top ridge grew taller.   Eucalyptus bosistoana height growth was influenced by different factors at different sites (Table 5). At site B, CVPLA, MPI, DIST, TPI, WEI and the interaction between WEI and DIST influenced tree height (Figure 4). In sites with local horizontal concave surfaces, trees were taller than the trees on horizontal flat or convex surfaces. TPI also showed a similar pattern, whereby trees were taller in valleys than on ridges. Until age 4.5 years, trees nearer the top ridge experienced faster height growth than trees further from the top ridge, but after that age, the converse was true ( Figure 4D). MPI and WEI were correlated positively and negatively, respectively, with height growth, suggesting that exposure to wind and other environmental stresses may have suppressed height growth. However, the lowest WEI with distant microsite had the highest height growth compare to low WEI and a position close to the ridge. On the other hand, high WEI with the farthest microsite, which means close to the valley floor, was the worst for tree height at the B site.
In the case of site C, E. bosistoana height was affected by WEI, WTI, TPI, MPI and DIST ( Figure 5). The MPI and WEI effects were similar to other results, with high MPI and low WEI resulting in increased tree height (Figure 5A,E). An increase of TPI affected the tree height up to age 2.5 years, after which the effect was reversed, such that trees in valleys had greater height growth, relative to trees on mid-slopes or ridges. The trees situated at mid-distance from the ridge top grew taller than those closest to, and furthest from, the ridge top. Interestingly, the surface wetness (TWI) minimally influenced tree height ( Figure 5B). Forests 2019, 10

Juvenile Survival Model
Analyses revealed that the smallest residual mean squares and the least biased residuals were produced by augmenting survival models (Equations (12) and (13)) with topographic attributes. The rate of mortality diminished with time in most plots, but mortality was greater in later years than in early years.
where, S EGT A and S EBT C are the survival proportion of E.globoidea and E. bosistoana at time T in site A and C; other terms have been defined earlier. Model parameters are provided in the supplementary material 5 (Table II). The residual distribution against predicted and independent datasets was normally distributed with minor distortions for all species and sites (Supplementary material, S3). Validation for both species was undertaken, and the survival proportion model had minimal differences in precision and bias between fitting and validation datasets ( Table 6). In the case of E. globoidea, the RMSE and MAE reduced during validation while they increased slightly with E. bosistoana model validation.

Key Variables Influencing Juvenile Microsite Survival
Eucalyptus globoidea survival was influenced by plan and profile curvature, WEI and distance from the ridge top (Table 7 and Figure 6). In concave and flat areas, mortality rate was steady whereas in convex areas mortality declined with time. This result was repeated for profile curvature, where on the raised surfaces trees survived in higher proportions than on hollow or flat surfaces. Survival decreased with increasing exposure to wind. Moreover, plots a long distance from the ridge top showed lower survival rates than those close to it. Eucalyptus bosistoana survival was influenced only by profile curvature ( Figure 6E). It showed that, in gullies, higher proportions of trees survived than on flat surfaces or ridges.

Juvenile Microsite Models
While earlier work has modelled juvenile trees on a broad scale (e.g., [30,32]), the juvenile microsite models described here have shown the utility of modelling juvenile crops at a finer scale. Individual juvenile trees have also been modelled by applying mathematical equations [25,82,83] and explaining different competing variables [82,84,85]. Kohama, et al. [86] and Weiskittel, et al. [21] studied juvenile and mature stand tree growth on a micro-scale, and Weiskittel, et al. [21] proposed a modelling framework but only for mature stand trees. However, juvenile and mature trees have different growth requirements and competition indices. The models presented in this study for juvenile trees have field applicability, which could be incorporated into a decision support system for silviculture at sites with similar characteristics, because previous research has shown that factors affecting juvenile crop growth and survival can have rotation-length implications [31,87,88].

Microsite Variables Affect Juvenile Tree Height Growth
This study showed that juvenile tree height growth and survival were affected by micrositerelated variables. Millner and Kemp [55] found similar results with other Eucalyptus species in Tuapaka, New Zealand. Topographic variables are major drivers of tree growth in many hilly regions [55,89], as they relate to both climatic and edaphic factors [90]. Height growth of both species in this study was greater in more sheltered microsites. Similar results were found by Brüchert, et al. [91] who showed that wind could influence the aerial architecture of the trees.
Generally, valley floors can be expected to have greater rooting depth [92] and less direct solar radiation as well as radiative heat [93], meaning that trees are better physiologically supported in terms of nutrients and moisture. However, this study found that mid-slopes, measured as distance from the ridge top, were best for E. globoidea. This may be due to mid-slopes provided optimum soil moisture availability to this species and it may be sensitive to the assumed elevated soil moisture and associated anaerobic soil conditions commonly found at valley floors [92,94].
Eucalyptus bosistoana grew taller in concave, depressed (valley) surfaces, and in locations farthest from ridges, which had relatively low wind exposure index. This can be explained in a similar way to E. globoidea but suggests that this species may be more water-demanding, or alternatively, tolerant of soil wetness than E. globoidea at young ages. Rohner, et al. [95] and Monserud, et al. [96] reported that the steep slopes resulted in shallow soil and less moisture availability due to lateral moisture flow. This is in line with the topographic position index effect, as it described each microsite with respect to slope.

Microsite Variability on Juvenile Tree Survival
Eucalyptus globoidea was apparently sensitive to assumed higher soil moisture levels and reduced growth, but the species could withstand drier conditions. Conversely, E. bosistoana showed less sensitivity to assumed higher moisture levels compared with E. globoidea. This suggested that E. globoidea may experience above-optimal levels of soil moisture (possibly due to anaerobic soil conditions) for tree survival, in both valleys and hollows. Conversely, E. bosistoana survived better in gullies, where there is presumably a chance to access higher moisture availability. Moreover, Ares and Marlats [89] found and concluded that in mountain regions of Argentina coniferous trees died on north-facing slopes, as this aspect receives more radiative heat than other aspects, which may increase water stress. Mason and Whyte [32] reported that frost negatively influenced juvenile tree survival, which could be an alternative or additional reason for increased mortality of E. globoidea in hollows.

Data Constraints
The initial height for the young Eucalyptus plantations was not recorded immediately after planting. For that reason, the initial height model was fitted by assuming Eucalyptus seedlings met New Zealand's Pinus radiata (D. Don) plantation standard, which was 0.25 m in height at time of planting [32]. The use of this standard height value might have influenced model stability at the early ages because the model extrapolated the height values for that age. Therefore, these models should be used cautiously over the period from planting to first measurement age. The genetics of each species may be a factor in the response to the environmental conditions. Gallart, et al. [23] showed that some P. radiata genotypes were more sensitive to microsite changes in soil physical properties than other, the same seed sources for both species were used across all three sites. Also, the design of these experiments was not orthogonal. Therefore, it is likely that there was some spatial autocorrelation between each replicate which may have influenced final results.
High resolution soil information describing physical and chemical characteristics were not available for these sites. Since the plantations were not established on sites with homogeneous soil conditions, this may be a confounding factor. High resolution soil data may have improved model precision and explanatory ability. Future studies should consider including fine-scale soil data, and plan accordingly for the associated time and financial costs.
Likewise, including fine-scale climatic variables into the models may have given greater explanatory power and understanding of causal processes [97,98]. However, in this study, only air temperature was considered but not selected for inclusion in the final model as it was not statistically significant. Solar radiation, wind speed, soil moisture and precipitation were not considered because data were not available at a sufficiently fine resolution as each site had only one meteorological station.

Summary and Conclusions
This study successfully demonstrated a statistically and biologically logical framework for modelling juvenile tree growth and survival at a microsite level. It also identified and explained height and survival variation of two Eucalyptus species at high spatial resolution on three sites. For both species, topographically sheltered surfaces yielded greater height growth and survival. Furthermore, there were different microsite factors that each species required for optimal growth and survival. The topographic features in this study indicated that soil moisture could be an important factor in explaining height growth and survival, with E. globoidea being negatively affected in valleys and hollows, where available moisture is presumably high. Conversely, E. bosistoana thrived in these environments. More research is required for these valuable Eucalyptus species over a larger environmental gradient of potential afforestation sites to understand the full suite of microsite factors and their interactions. However, the results of this study provide important information that will assist forest managers to identify sites and microsites where it is best to establish each species.
The study indicates that researchers wishing to evaluate effects of microsite variation on tree crop growth are likely to find high resolution digital elevation models more useful than randomly allocated soil pits and sparse soil analyses. Indices of land form variation derived from the digital elevation model were surprisingly well correlated with tree crop performance.