More timber from fewer trees – determining what tree density optimises silver beech merchantable yield based upon a long-term thinning trial Tomás A. Easdale1*, Robert B. Allen2, Larry E. Burrows1, David Henley3 and Dudley A. Franklin4 1 Manaaki Whenua – Landcare Research, PO Box 69040, Lincoln 7640, New Zealand 2 8 Roblyn Place, Lincoln 7608, New Zealand 3 Scion, PO Box 29237, Christchurch 8440, New Zealand 4 7 Mt Thomas Rd, Rangiora 7471, New Zealand *Corresponding author: EasdaleT@landcareresearch.co.nz (Received for publication 24 August 2021; accepted in revised form 22 May 2022) Abstract Background: The tree stem density which optimises merchantable timber yield (volume per unit area) is unknown for most of New Zealand’s indigenous tree species. While moderate thinning of even-aged stands can promote yield, intense thinning may decrease yield by creating space that cannot be filled by residual trees, increasing tree mortality or reducing tree height. We quantified the effects of density on silver beech (Lophozonia menziesii (Hook.f.) Heenan & Smissen) tree growth, height and mortality, identified the density leading to optimal merchantable yield and assessed if this density varied with stand age. Methods: Tree stem diameter growth, height, and mortality responses to density were determined using tagged individuals monitored over time on a long-term thinning trial combined with flexible, multilevel, non-linear models. Empirical stand yield responses to density were determined and compared to yield–density relationships in simulated stands. The stand simulations projected beyond the monitored stand ages using the tree-level responses fitted to empirical data. Results: Low densities (≤400 stems ha-1) sustained fast tree growth for longer than high densities (≥700 stems ha-1) after thinning, but density did not consistently affect merchantable tree heights. The probability of tree mortality increased after intense thinning, but only temporarily, and never exceeding c. 0.01 year−1. A regression of yield–density relationships identified an empirical optimum of c. 570 stems ha–1 for stand ages of 48 and 58 years. At this density, merchantable yield at 58 years was seven-fold greater than that in unthinned stands. The simulations suggested moderately higher densities for optimal yield than our empirical optimum, a moderate increase in optimal densities with stand age, and that c. 90 % of potential cumulative yield was attained at 80 years. Conclusions: Because thinning increased tree growth, but had minimal effect on tree mortality, our results alleviate concerns about the stability and productivity of thinned stands. Densities that optimise yield are about two-fold greater than those previously recommended for silver beech and they remain relatively stable as stands age. This suggests that a single density will be adequate for a range of harvest ages, although harvest should take place before a stand age of 80 years. Such conclusions are relevant to managing regeneration within coupes harvested under existing legislation and to areas planted with silver beech. New Zealand Journal of Forestry Science Easdale et al. New Zealand Journal of Forestry Science (2022) 52:16 https://doi.org/10.33494/nzjfs522022x179x E-ISSN: 1179-5395 Published on-line: 07/06/2022 © The Author(s). 2022 Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (https://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 and ecosystem loss have revitalised the long-standing pursuit of efficient production from natural resources (Trewavas 2001; Keating et al. 2010). There are Introduction Human population growth and per capita consumption, land-use impacts on soils and concerns about biodiversity Keywords: Individual-based simulation; mixed-effects models; Nothofagaceae; resource-use efficiency; southern beeches; tree growth; tree mortality mailto:easdalet@landcareresearch.co.nz http://creativecommons.org/licenses/by/4.0/), Easdale et al. New Zealand Journal of Forestry Science (2022) 52:16 Page 2 compelling reasons for an eco-efficient transformation where more is produced from limited land, water, nutrients and energy (Keating et al. 2010). This raises a need for information, practices and technologies that can maintain or increase yield via efficient use of resources (Keating et al. 2010). A key contribution from plant breeding research has been the selection of crops that allocate a higher fraction of plant productivity to harvestable yield and that reduce resources wasted in plant competition (Weiner 2003). By examining the effects of thinning on tree density and yield, forestry science has long been making a similar contribution. Much has been learned but challenges remain (see Zeide 2001; Pretzsch 2005; Deng et al. 2012). Once a forest canopy is closed and trees make full use of available resources, the volume of wood produced per unit area (hereafter yield) by even-aged stands is nearly constant across a wide range of tree stem densities (hereafter densities), regardless of whether there are many small trees or a few large trees (Langsaeter 1941; Zeide 2001). This can be seen in thinned, even- aged stands (Nishizono 2010) or natural self-thinning stands (Xue & Hagihara 1998) and is consistent with the constant final yield of undisturbed monocultures (Weiner & Freckleton 2010). Thus, thinning dense even-aged stands is a means of redistributing resources to promote growth of residual trees and increase merchantable yield of timber (Zeide 2001). There are, however, associated trade-offs; although moderate thinning will increase merchantable yield, intense thinning will not. Excessive thinning can compromise yield by: (1) creating large canopy spaces that cannot be filled by the crowns of residual trees (Smith et al. 1997; Zeide 2001); (2) promoting lateral growth that produces short trees with less usable timber (Clutter et al. 1983; Smith et al. 1997); and, (3) increasing residual tree mortality (Harrington & Reukema 1983; Kariuki 2008) leading to further canopy gaps and inefficient use of space. While it is well established that merchantable yield has a hump-shaped relationship with density, it is generally unknown what density leads to maximum merchantable yield for a given species and site condition (Zeide 2001; Pretzsch 2005) as well as how the density giving optimal merchantable yield changes with stand age (Fig. 1; Pretzsch 2005; Zeide 2008). The cost of thinning interventions, and how long it takes for trees to reach harvesting size, underscores the importance of determining what density is optimal for yield. We investigate how density influences tree- and stand-level responses by silver beech (Lophozonia menziesii (Hook.f.) Heenan & Smissen) using a long-term thinning trial. More harvest volume is produced from silver beech than any other indigenous tree species yet this is a species that tolerates competition and grows slowly in natural forests (Richardson et al. 2011). The effects of thinning on residual trees are somewhat contentious. Although a shade-tolerant species with low apical dominance, such as silver beech, can become short-statured at low densities (Wardle 1984), the small diameters of trees in dense stands can lead observers to under-estimate tree age and over-estimate height growth at high densities (Smith et al. 1997). While strong competition in dense, even-aged stands leads to intense self-thinning driven tree mortality (e.g., Osawa & Allen 1993), intense thinning of such stands can raise mortality through physiological shock or by destabilising residual trees and increasing the risk of windfall (e.g., Harrington & Reukema 1983; Kariuki 2008). The first aim of our study was to determine how density influences tree- level stem diameter growth, height, and mortality as well as how this varies with stand age and tree size. Predictions are also contentious at a stand-level. Since tree competition can become more intense as tree size increases, we might expect that densities for optimal yield will decrease in older stands. However, long-term assessments of thinned Fagus sylvatica stands indicate that optimal merchantable yield shifts to higher densities during stand development (Pretzsch 2005). The second aim was to distinguish between these alternative stand- level responses by answering the following questions (see Fig. 1 for graphical representation): (1) what density optimises merchantable (harvestable ‘crop’) yield?; (2) what is the resulting yield?; (3) how does this yield compare with that of unthinned stands?; and, (4) how does the density giving optimal yield change with stand age? We address these questions with a combination of empirical analyses and simulations of yield and density relationships over time. Finally, we consider the management implications for New Zealand beech forests. Density M er ch an ta bl e yi el d ?1 ?2 ?3 ?4 Time1 Time2 FIGURE 1: Although the relationship between merchantable timber yield and density is of a known shape, key questions relate to the density that optimises merchantable yield for a given species, site condition and stand age (?1), ensuing yield (?2), yield gains relative to untended stands (?3), and temporal changes in these responses (?4). All four questions are addressed here for silver beech. Methods Study site Our study focuses on a silver beech forest in the Alton Valley (46°02’S 167°37’ E and 150-190 m elevation), Southland, New Zealand. The site has Orthic Brown soils (Dystrudept in US soil taxonomy; Hewitt 2010) and the terrain is mostly flat. The area receives, on average, 1250 mm of annual rainfall and the mean annual temperature is c. 9.5˚ C. The original old-growth forest dominated by silver beech was felled in 1951 but c. 40 seed trees per hectare were retained. Seed trees were poisoned 12 years after logging. Subsequent natural regeneration was variable but typically led to a high density of saplings (Franklin 1981; Easdale et al. 2009). Silver beech dominated the regenerating forest although mountain beech (Fuscospora cliffortioides (Hook.f.) Heenan & Smissen) was present locally and there were also low numbers of other indigenous angiosperm and podocarp trees (Easdale et al. 2009). Silver beech silvics Silver beech is a long-lived (up to 600 years) evergreen tree species that occurs from c. 37˚ 30’ to 46˚ 30’ S (Wardle 1984). Trees up to 42.7 m tall and 3.0 m diameter have been measured. Silver beech is restricted to montane and subalpine forests in the north of its range (up to 1400 m elevation), but in the south it is found near sea-level. It spans a wide range of annual rainfall from >8000 mm in western parts of the Southern Alps, where it can be the dominant species, down to c. 600 mm in south-eastern parts of the South Island, where it comes close to forming the driest New Zealand beech forests (Wardle 1984). This species seeds prolifically, and in 33% of years >2000 seeds m−2 are produced which commonly lead to an abundance of seedlings on the forest floor (Burrows and Allen 1991; Wardle 1984). Even seedlings of this shade-tolerant species which have been suppressed for decades are capable of responding to increased resource availability (Wardle 1984). Asymmetric competition for light is the major factor controlling diameter growth of small silver beech trees throughout New Zealand’s South Island, whereas physical environment (e.g., elevation) is relatively more important to large tree growth (Easdale et al. 2012), thus thinning of regenerating stands is expected to enhance growth of residual trees. Experimental design and sampling protocols A thinning trial of 15 contiguous 0.2 ha stands was established in the Alton Valley site in 1971. Stand treatments comprised different combinations of thinning interventions implemented in 1971 and 1980 to give a wide range of densities by 1980, with one unthinned control stand (Fig. 2; Franklin 1981). Other than this unthinned stand with c. 8000 stems ha–1 (see below), total stand-level densities after the second thinning in 1980 were 3000 (1 stand), 1500 (1), 730 (1), 390-415 (4), 325 (1), 285-295 (3) and 190-200 (3) stems ha–1 (Fig. 2). Trees were thinned mostly from below, aiming to retain dominant ‘crop’ trees with an even spacing while discarding forked trees where possible. Selected trees were ≥ 3 cm DBH in 1971 and ≥ 5 cm DBH in 1980 but smaller stems were retained at densities ≥ 1500 stems ha–1. Crop trees were pruned in 1975, 1977 and 1979 to a final height of 5.5 m. Easdale et al. New Zealand Journal of Forestry Science (2022) 52:16 Page 3 TABLE 1: Description of the study sites FIGURE 2: Layout of the Alton Valley thinning trial. Fourteen 0.2-ha stands were thinned in 1971 and 1980 to densities in stems ha–1 shown in the figure, with one left unthinned. 190 190 200 410 415 390 290 395 295 325 730 285 50 m Crop trees (n = 1431, including 54 mountain beech trees) were tagged in 1971 with tags drawn-out by copper wires and those surviving were remeasured in 1974, 1977, 1980, 1982, 1985, 1995, 1999 and 2009. The selected crop trees had initial diameters at breast height (DBH) of 8.2 ± 3.1 cm (mean ± SD) and total heights of 5–8 m in 1971, as well as a DBH of 14.1 ± 4.9 cm at the time of the second thinning in 1980. Only a subset of crop trees, with large DBH and a regular spacing, were initially tagged in the two stands with highest final densities and in the control plot, but all trees ≥10 cm DBH were measured in all stands in 1999 and 2009. Even though we used metal detectors to relocate fallen tags in 2009, 133 tags were not found and these trees could not be matched to prior measurements. A subset of crop trees had merchantable height (to the lowest major fork, bend, or the point where tapering decreased to c. 50-60% of DBH) and total height measured with ultrasonic hypsometers (Vertex III, Häglof, Sweden) in 1999 and 2009. The height measurements were made on five randomly selected trees within each of four DBH classes to allow for calibration of height–DBH curves across the range of tree sizes found in each stand. Data analysis To assess the extent to which thinning affected tree- level growth and whether low density compromised the height and survivorship of residual trees, we (step 1) assessed empirical responses of DBH, merchantable height, and mortality to 1980 densities using the repeated measurements of crop trees and available height measurements. We then (step 2) assessed empirical stand-level yield responses to 1980 thinning densities using the full measurements of stems ≥10 cm DBH in 1999 and 2009 (stand ages of 48 and 58-years; details below). Although the trial spanned densities of 190 to c. 8000 stems ha–1 with intermediate densities, 38 years of monitoring, and eight remeasurements, only eight densities were trialled and not all stands had attained merchantable yields by the last measurement. Thus, as a last step (3), we combined the tree-level models developed in step 1 into simulations of stand development to predict merchantable yield across a wider range of densities and stand ages than those measured. In assessing yield, we followed local industry standards (Ministry for Primary Industries 2013) that define merchantable yield as the combined volume per unit area of clear tree stems ≥ 30 cm DBH up to a height that excludes major forks or stem tapering. Tree-level responses We modelled the cumulative growth (DBH as a function of stand age) for silver beech crop trees that survived to 2009. This encompassed a total of 1153 trees with eight or nine repeated measurements and, except for one intensively thinned stand with 14 trees, included 32 to 157 modelled trees per stand (stems with unmatched tags or with <7 DBH measurements were excluded from growth analyses). Tree growth analyses relied on non-linear multilevel (‘mixed’) models to account for the covariance structure that results from repeated Easdale et al. New Zealand Journal of Forestry Science (2022) 52:16 Page 4 measurement of crop trees within stands (Gelman & Hill 2007). We modelled the effects of density on DBH in three steps. First, we fitted tree dbh (cm) as a Gompertz function of stand age (years): (1) with a = (a0 + aj + aij) and d = (d0 + dj + dij) where the parameters a, b and d respectively define the asymptote, location and slope of the curve (Sit & Poulin-Costello 1994). Multilevel growth curves were simultaneously fitted for individual trees i and for the average tree in each stand j by allowing a and d to vary both at stand-level (aj and dj) and at individual tree- level (aij and dij) via nested random effects and assuming normally distributed residuals εij. Here, the slope of the fitted curves represents DBH growth (Husch et al. 2003). Second, we extracted the stand-level parameters aj and dj fitted in the first step and, from a candidate set of models, identified two models that best explained their respective relationship with density after second thinning (1980). Third, we substituted parameters a and d into the original Gompertz equation with the models resulting from the second step, as shown below, and fitted this modified Gompertz model with random (tree- and stand-level) effects on the parameters controlling the asymptote and slope of the curve (i.e. DBH growth). This three-step procedure allowed us to explicitly model and calibrate the non-linear relationship between diameter curves and tree density within stands. We modelled merchantable height (Hmerch) using the 320 silver beech crop trees measured in 1999 and/or 2009. To this end, we first used multilevel models to identify the height–DBH function having the strongest fit among candidate functions given by Husch et al. (2003) and then tested the effects of density on fitted stand- level parameters (second step above). Inspection of crop tree mortality rates computed for each density and census interval showed that the probability of mortality increased exponentially towards lower densities (heavy thinning) and also increased at first (up to 19 years after the second thinning) and then decreased with time up to 29 years after the second thinning (Fig. S1, Supplemental information). These also suggested that timing since thinning may be more important than stand age in determining tree mortality. Thus, we modelled the probability of tree mortality (pa) as an exponential effect of density with a flexible lognormal function to account for the humped response to t_thin, time since second thinning: (2) where a, b, c, f and g are estimated parameters. Since the data consisted of multi-year census intervals (t) of different duration (i.e. periodic observations of survivorship/mortality), and periodic (multi-year) mortality (Pm) relates to annual mortality pa as: (3) we modelled annual probabilities of mortality by substituting pa in Equation 3 with Equation 2: (4) In this way, Equation 4 allowed us to estimate parameters for annual mortality (a, b, c, f and g) from periodic mortality data of variable duration. We accounted for all survivorship/mortality records (including any mountain beech selected as crop trees) so long as tagged stems were relocated in consecutive measurements and accounted for mortality events only where a stem was confirmed dead (i.e. stem measurements with tags ‘not found’ were excluded from analysis). Maximum likelihood methods were used to estimate the parameters most likely to have produced the data, given the models. For DBH growth and merchantable height, parameters were estimated with non-linear multilevel models and a Gaussian distribution using the “nlme” library in R (R Development Core Team, 2013). Mortality was modelled with non-linear models and a binomial distribution based on simulated annealing, a robust global optimisation algorithm (Goffe et al. 1994) implemented using the “likelihood” package in R (Murphy 2012). Alternative models for each response variable were compared using their associated Akaike Information Criterion (AIC) values, where lower AIC values indicate greater empirical support for a model. Differences in AIC values (∆AIC) < 2, between 4 and 7, and >10 respectively indicate negligible, moderate and strong empirical support between alternative models (Burnham & Anderson 2002). As the 1980 density in the unthinned control stand was unavailable, we: (i) fitted a regression for the relationship between tree density post-thinning in 1980 and density of live stems ≥10 cm DBH in 2009 (df =13, R2 = 0.92) for all but the unthinned stand and then: (ii) back-estimated tree density in 1980 from the fitted model and the estimated density of the unthinned stand in 2009. Stand-level responses Merchantable yield was first predicted for each tree from merchantable height and DBH using equations calibrated for Nothofagaceae (hereafter southern beech) found in New Zealand (Ellis 1979) and then added up by stand. Only trees ≥ 30 cm DBH were used in calculating merchantable yield (Ministry for Primary Industries 2013). To identify the tree density that optimised yield, we first fitted a flexible third-degree polynomial curve to describe the relationship between merchantable yield and density. We then identified the density which corresponded with the peak of the merchantable yield curve using the general-purpose “optim” optimizer function in R (R Development Core Team, 2013). Lastly, to assess how the density giving optimal yield varied over time, we developed a stochastic simulation model of stand development that combined the empirical growth and mortality responses derived for individual trees. This involved simulating the cumulative stem growth of individual trees at various stem densities and, for each simulated stand, estimating merchantable volumes as the combined volume of all live trees at a given stand age. Simulated growth curves varied as a function of stem density and stand age and built-in growth deviations around the ‘mean curve’, based upon our previously calibrated empirical responses. Incorporating growth deviations was important since trees only contribute to yield once they reach merchantable size (≥30-cm DBH) and small fluctuations in DBH distribution can cause large changes in yield (Smith et al. 1997). Merchantable heights were determined by a calibrated height-DBH allometry function and both height and DBH then informed the corresponding merchantable volume for each tree, according to a published allometry (Ellis 1979). The model also accounted for tree mortality, with trees stochastically surviving or dying based on previously estimated probabilities of mortality, as a function of stand age and time after thinning. For consistency with the field trial, we assumed that final thinning took place at 29 years. We: (1) ran simulations for post-thinning densities spanning 100 to 10000 stems ha−1, with multiple runs for each density, and concluded by: (2) identifying what post-thinning density optimised merchantable yield for each stand age. Our simulations were conditional upon the empirical tree-level results so that further specifics of the simulation approach are presented below, after the corresponding empirical results. Results Tree-level responses Tree growth Growth curves fitted to the empirical data showed that progressively larger DBHs resulted from gradually lower densities (Fig. 3a). Intense thinning led to longer-lasting growth responses, with responses to ≤400 stems ha–1 only starting to become discernible from each other after some 10 years from second thinning (39 years of age) (Fig. 3a). At 58 years, the mean DBH of crop trees in two of the lowest stocked stands was twice that of the unthinned control stand (37 cm DBH at 190 stems ha–1 vs. 17 cm DBH at c. 8000 stems ha–1). Tree growth was positively autocorrelated, with DBH at 20 years predicting DBH at 58 years (cross-stand mean r = 0.72). Starting with Equation 1, we progressively built in the effects of density on tree DBH growth. Parameter a was best described by a power function of density at 29 years, namely, the age at second thinning (∆AIC = −24 relative to a null model). Adding a term for log density at 20 years, the age at first thinning, did not improve predictions (∆AIC = −25 relative to a null model). Parameter d was best described as a logarithmic function of density at 29 years (∆AIC = −22 relative to a null model). Adding a term Easdale et al. New Zealand Journal of Forestry Science (2022) 52:16 Page 5 for log density at 20 years only moderately improved predictions (∆AIC = −26 relative to a null model) so we simply modelled d as a function of density at 29 years. Substituting parameter a and d in Equation 1 with power and logarithmic functions: and gave a modified Gompertz function for tree growth based on stand age and density after second thinning: (5) We parameterised Equation 5 with nested tree- (i) and stand-level (j) random effects on parameters a2 and d2 (i.e., and ) and obtained an unbiased model with the parameter values given in Table 1. The fitted model explained substantial variance when only accounting for the fixed parameters (R2 = 0.63) and captured nearly all of the variance with both the fixed and random parameters (R2 = 0.99). This meant that the random tree-level effects captured much of the DBH variability that could not be explained by density. Inspection of tree-level random parameters extracted from the “nlme” model showed that a2 and d2 were negatively correlated (r = −0.329) and that a2 had a negatively skewed distribution, which we overcame by a logarithmic transformation: where the negative sign turned a left skew into a right skew, the internal constant minimised the skew and the external constant centred the transformed values. Further assessment indicated that the standard deviation (sd) of tree-level parameters tra2 and d2 increased towards stands with lower stem density. After testing simple linear and non-linear functions (via “nls” in R), we modelled the standard deviation (sd) of tra2and d2 as a power function of ln density with fitted parameters given in Table 1. These models described much of the ‘random’ within-stand variance that remained unaccounted for by the ‘fixed’ parameters. Tree height Of the various height–DBH functions tested, a simple logarithmic function gave the best combination of balanced residuals and strength of support (∆AIC = −18 relative to a null model (i.e. fixed height estimates)). The resulting model was unbiased but had minimal explanatory value (R2 = 0.04) and gave a 0.7 m increase in merchantable height between a 30-cm and a 55-cm-DBH tree (about the largest recorded at last measurement). The fitted model and estimated parameter values are given in Table 1. A logarithmic model with variable slope parameters for stands had stronger support than the above (∆AIC = −34 relative to a null model), with up to 1.6 m differences in merchantable height among stands for a 30-cm-DBH tree (Fig. 3b). However, the fitted slopes (Fig. 3b) were unrelated to density at 20 or 29 years. Easdale et al. New Zealand Journal of Forestry Science (2022) 52:16 Page 6 10 20 30 40 Stand age (years) D B H (c m ) 20 30 40 50 60 7915 3000 1500 730 415 290 325 200 395 410 295 390 285 190 190 (stems ha−1) Density at 29 yr(a) 6 7 8 9 DBH (cm) M er ch an ta bl e he ig ht (m ) 30 40 50 60 415 325 190 285 410 295 390 200 190 1500 3000 7915 290 730 395 (stems ha−1) Density at 29 yr(b) FIGURE 3: Cumulative DBH–age growth curves (a); and merchantable height–DBH curves (b); fitted for the average crop trees across a range of densities at second thinning (29 years). Vertical dotted lines indicate the timing of the two thinning interventions in (a). The grey saturation of curves scales with the logarithm of density at 29 years. Tree mortality The probabilities that trees died were low. Only 26 of 1379 crop trees with original tags found were confirmed to have died within 29 years of monitoring after the second thinning. Mortality of ‘crop’ trees was lowest in dense stands and highest in stands with the lowest densities (Fig. S1, Supplemental information). The highest annual probability of mortality was 0.0106 per year for stands with 190 stems ha–1 between 15 and 19 years after second thinning. A model with an exponential effect for density and a flexible response to t_thin after second thinning (Fig. 4) had strong support relative to a constant probability of mortality (∆AIC = −30) and captured the empirical pattern of mortality effectively (Fig. S1, Supplemental information). Estimated parameter values are given in Table 1. Stand-level responses Empirical density and optimal yield A regression of merchantable yield calculated from field measurements (for all trees ≥30 cm DBH) showed that yield peaked at 62 m3 ha–1 with 571 stems ha–1 in 48-year-old stands (Fig. 5a). Subsequently, the optimal yield became more pronounced with an estimated peak yield of 114 m3 ha–1 at c. 567 stems ha–1 in the 58-year- old stands (Fig. 5a). At this age, merchantable yield was c. seven times greater than in the unthinned control stand. An assessment of stand data using all crop and non- crop trees ≥10 cm DBH (Table 2) showed some variation from the above results for crop trees only. Mean stem diameters were consistently larger at low density relative to high density but differences were partly dampened due to stem recruitment in the more heavily thinned stands (note an increase in stem densities from 1999 to 2009 at low densities). Mean top heights had differences of up to 4 m between stands and a weak positive association with density, with heights of 15.1- 18.0 m for densities of 400 stems ha−1 or less and heights Easdale et al. New Zealand Journal of Forestry Science (2022) 52:16 Page 7 TABLE 1: Set of fitted models and their associated parameter estimates, as calibrated from repeated measurement of tagged silver beech trees at the Alton Valley trial. These models and parameters where subsequently employed to simulate merchantable yield across a broad range of stand ages and densities. ra2-d2 is the correlation coefficient between tree-level random parameters a2 and d2. Variable Model Equation Parameter Estimated value Standard Error DBH (5) a1 138.0 23.7 a2 −0.2369 0.0282 b 1.680 0.00570 d1 0.01962 0.00792 d2 0.007785 0.00130 DBH variance (6) a 0.5917 0.430 b −0.5915 0.404 (7) c 0.003343 0.00177 d −0.4177 0.294 ra2-d2 −0.329 – Merchantable height (8) a 3.124 0.845 b 1.142 0.256 Probability of Mortality (2) a 1.327 1.51 b 0.2424 0.101 c 0.3850 1.29 f 14.98 2.37 g 0.5181 0.121 30 25 20 15 10 5 0 Density at 29 years (stems ha−1) Ye ar s af te r s ec on d th in ni ng 5e−04 0.001 0.0015 0.002 0.0025 0.003 200 500 1000 2000 5000 Density at 29 years (stems ha−1) Thinning FIGURE 4: Contour plot for probability of tree mortality as a function of years after the second thinning and density after the second thinning (29 years). Both contour lines and grey saturation of the background follow the ‘topography’ of mortality probabilities, with higher saturation for higher probabilities of mortality. of 17.1-19.4 m for densities over 700 stems ha−1. Basal areas generally corresponded with the merchantable yields reported above, with a maximum of 46.1 m2 ha−1 noted at 750 stems ha−1 at 58 years. Stand simulations and changes in optimal density with stand age With tree-level responses understood, we then assembled a stochastic simulation model. For a given density, mean DBH growth was given by Equation 5 (Table 1). Growth variability was incorporated by: (a) sampling parameters tra2 and d2 from a Gaussian distribution with sample size given by density, mean value of 0, and the standard deviation defined by eqs. 6 and 7 (Table 1): (b) back-transforming tra2 into a2; (c) shuffling values of d2 until they reached the correlation with a2 detected in the field data (r = −0.329) and (d) adding up the values obtained from the previous step to the mean values of a2 and d2 given in Equation 5. Each set of parameters described the cumulative growth of each tree in a simulated stand (Fig. S2, Supplemental information). Merchantable heights were defined from DBH by Equation 8. Simulated trees died if a random number between 0 and 1 drawn each year for each tree was lower than the fitted probability of mortality given by Equation 2 (Table 1) for a given density and t_thin, time since second thinning. Predictions of merchantable yield in simulated stands corresponded relatively well with the empirical analyses (Fig. 5). As in the field estimates, differences in such yield across a range of densities were small early on and became more pronounced with stand age (Fig. S3 and S4, Supplemental information). Simulations indicate that a density of 650 stems ha–1 gives maximum merchantable yield for 40 to 53 year-old stands; the density giving optimal yield shifts to 850 stems ha–1 for 54 to 64 year- old stands and stabilises at 950 stems ha–1 over 65 years (Fig. S3, Supplemental information). At a density giving optimal yield, simulations predict a mean merchantable yield of 126 m3 ha–1 at 58 years (Fig. 5b) and 190 m3 ha–1 at 80 years. The dome in the yield response curve was however, sufficiently flat to predict similar yields across a range of densities (Fig. S4). The temporal increase in merchantable yield is asymptotic with c. 90 % of the aggregate yield attained at 80 years. Merchantable yield is predicted to be four to four and a half times greater when compared to unthinned stands with c. 8000 stems ha–1 (Fig. 5b). Discussion Tree-level responses The implications of controlling density in even-aged, silver beech stands are clear and confirm initial assessments that the species responds well to thinning (Baker & Benecke 2001). Growth trajectories overlap early on as crop trees grow freely but they begin to diverge once neighbouring trees begin to compete. The higher the density of trees, the sooner competition is manifested in reduced growth (Weiner & Freckleton 2010) and this shade-tolerant species does not self-thin up to at least 58 years. Competition can affect growth at surprisingly low densities (Clutter et al. 1983) but the effects only become evident later on (>40 years at 190 stems ha–1). This supports a view that assessing Easdale et al. New Zealand Journal of Forestry Science (2022) 52:16 Page 8 ● ● ● ● ● ● ● ● ● ● ● ● ● 100 200 500 1000 2000 5000 10000 0 50 10 0 15 0 Density at 29 years (stems ha−1) M er ch an ta bl e yi el d ( m 3 ha −1 ) ●● ●● ●● ●● ●● ●● ●● ●● ●● ●● ●●●● ●● 48 yrs 58 yrs ● ●● R ec or de d Fi tte d Empirical results Simulation results(a) (b) ● ● ● ● ● ● ● ● ● ● ● ● ● 100 200 500 1000 2000 5000 10000 0 50 10 0 15 0 Density at 29 years (stems ha−1) M er ch an ta bl e yi el d ( m 3 ha −1 ) ●● ●● ●● ●● ●● ●● ●● ●● ●● ●● ●●●● ●● 48 yrs 58 yrs ● ●● R ec or de d S im ul at ed FIGURE 5: Empirical estimate of merchantable yield for the last two measurements (48 and 58 years) (a). Optimum yields with respect to density at second thinning (29 years) were fitted with third-order polynomial regressions. Dotted lines are 95% confidence intervals for the regressions. Simulated yield results (lines) for 48- and 58-year-old stands and corresponding empirical estimates (points) (b). Solid lines represent the mean of 100 simulation runs and dotted lines the 95% confidence intervals. Easdale et al. New Zealand Journal of Forestry Science (2022) 52:16 Page 9 TA B LE 2 : S um m ar y of s ta nd m et ri cs fo r al l l iv e st em s ≥1 0 cm D B H in t he A lt on t ri al a t 48 a nd 5 8 ye ar s (1 99 9 an d 20 09 r es pe ct iv el y) . M ea n to p he ig ht w as e st im at ed fr om m od el le d he ig ht s fo r th e 10 0 la rg es t st em s pe r he ct ar e. M ea n D B H a nd m ea n to p he ig ht a re p re se nt ed w it h th ei r as so ci at ed 9 5% c on fid en ce in te rv al s. P lo ts a re so rt ed fr om lo w es t t o hi gh es t p os t- th in ni ng d en si ty in 1 98 0. St em d en si ty ( st em s ha −1 ) M ea n D B H ( cm ) M ea n to p he ig ht ( m ) B as al a re a (m 2 h a− 1 ) Cr op tr ee s St em s ≥1 0 cm D B H St em s ≥1 0 cm D B H St em s ≥1 0 cm D B H St em s ≥1 0 cm D B H P lo t 19 71 19 80 19 99 20 09 19 99 20 09 19 99 20 09 19 99 20 09 4 30 00 19 0 17 0 38 0 35 .4 5 (3 3. 2 - 3 7. 9) 27 .8 2 (2 4. 6 - 3 1. 2) 15 .3 2 (1 5. 2 - 1 5. 5) 16 .5 7 (1 6. 5 - 1 6. 7) 17 .4 29 .1 5 15 00 19 0 17 0 60 0 31 .5 5 (2 8. 9 - 3 4. 1) 20 .3 1 (1 8. 6 - 2 2. 8) 14 .0 8 (1 3. 9 - 1 4. 2) 15 .2 ( 15 .1 - 15 .3 ) 14 .1 25 .8 8 75 0 20 0 20 0 52 0 28 .9 8 (2 6. 3 - 3 1. 2) 20 .7 6 (1 8. 7 - 2 3. 1) 15 .1 8 (1 5. 1 - 1 5. 3) 16 .3 6 (1 6. 3 - 1 6. 4) 14 .1 22 .6 2 79 15 28 5 26 5 43 0 27 .6 5 (2 5. 8 - 2 9. 3) 24 .7 3 (2 2. 5 - 2 7) 14 .1 ( 14 - 14 .3 ) 15 .1 ( 15 - 15 .3 ) 16 .8 24 .5 14 30 00 29 0 30 0 52 5 28 .1 6 (2 6. 4 - 3 0. 3) 23 .0 5 (2 0. 9 - 2 5. 4) 15 .8 5 (1 5. 7 - 1 6) 16 .3 6 (1 6. 2 - 1 6. 5) 20 .2 27 .6 15 15 00 29 5 30 0 59 5 26 .9 7 (2 5. 6 - 2 8. 3) 23 .2 5 (2 1. 4 - 2 5. 3) 14 .7 2 (1 4. 6 - 1 4. 8) 17 .9 ( 17 .8 - 18 ) 17 .8 30 .2 13 75 0 32 5 31 0 55 5 27 .3 5 (2 5. 6 - 2 9. 2) 23 .4 4 (2 1. 4 - 2 5. 5) 15 .1 4 (1 5 - 1 5. 3) 17 .9 5 (1 7. 8 - 1 8. 1) 19 .5 29 .3 12 15 00 39 0 37 0 50 5 27 .5 5 (2 6. 1 - 2 9. 2) 26 .5 6 (2 4. 4 - 2 8. 5) 15 .2 3 (1 5. 1 - 1 5. 4) 17 .5 3 (1 7. 4 - 1 7. 7) 23 .4 32 .4 11 30 00 39 5 37 0 45 0 28 .8 ( 27 .4 - 30 .4 ) 28 .5 9 (2 6. 7 - 3 0. 8) 16 .7 1 (1 6. 6 - 1 6. 9) 17 .9 4 (1 7. 8 - 1 8. 1) 25 .5 32 .6 1 30 00 41 0 41 0 50 0 27 .6 8 (2 6. 8 - 2 8. 7) 27 .8 7 (2 6. 1 - 2 9. 3) 15 .2 8 (1 5. 2 - 1 5. 3) 16 .0 8 (1 6 - 1 6. 1) 25 .3 32 .9 7 15 00 41 5 40 5 52 0 25 .5 7 (2 4. 4 - 2 6. 8) 25 .4 6 (2 4 - 2 7) 14 .1 ( 14 - 14 .2 ) 15 .3 4 (1 5. 3 - 1 5. 5) 21 .8 29 .3 3 30 00 73 0 72 5 93 0 23 .4 8 (2 2. 5 - 2 4. 7) 23 .6 4 (2 2. 5 - 2 4. 9) 16 .1 4 (1 6 - 1 6. 2) 17 .5 8 (1 7. 5 - 1 7. 7) 34 46 .1 9 15 00 15 00 83 0 88 0 22 .4 7 (2 1. 5 - 2 3. 4) 24 .0 3 (2 3 - 2 5. 1) 16 .8 8 (1 6. 8 - 1 7) 19 .4 1 (1 9. 3 - 1 9. 5) 35 .8 43 .8 10 30 00 30 00 12 85 13 60 17 .8 2 (1 7. 1 - 1 8. 6) 18 .7 5 (1 8 - 1 9. 6) 15 .7 ( 15 .6 - 15 .8 ) 17 .0 6 (1 7 - 1 7. 2) 35 .8 42 .4 6 ~ 80 00 ~ 80 00 19 25 20 55 14 .9 4 (1 4. 5 - 1 5. 4) 15 .8 4 (1 5. 4 - 1 6. 3) 14 .5 ( 14 .4 - 14 .6 ) 18 .3 3 (1 8. 2 - 1 8. 5) 36 .9 44 .2 Easdale et al. New Zealand Journal of Forestry Science (2022) 52:16 Page 10 thinning responses requires long-term monitoring. At 58 years, differences in mean DBH, between highest and lowest densities, were two-fold (17 cm to 37 cm) and, with growth at 190 stems ha–1 little constrained, the ratio might expand further in older stands. For some hardwoods, total height growth can be significantly less at lower densities than at moderate or higher densities (Clutter et al. 1983) despite a potential increase in mean height derived from thinning shaded and suppressed trees. We found only moderate statistical support for reduced total heights at low densities (result not shown) but no significant change in merchantable heights as a function of density. This relieves concern about costs from possible shifts in stem allometry for silver beech at low density (Wardle 1984) and supports a view that dominant tree heights are least affected by thinning (Lanner 1985; Skovsgaard & Vanclay 2007). The basis of thinning is to concentrate growth on selected trees either by avoiding untimely timber losses in stands that self-thin (Smith et al. 1997) or by transferring growth from many small slow-growing trees in stands that are slow to self-thin (e.g. Farnden & Herring 2002). Stands of silver beech tend to be dense and slow to self-thin, leading to stagnant growth both in regenerating stands (Wardle 1984) and in old-growth forests (Richardson et al. 2011; Easdale et al. 2012). The low mortality recorded in the unthinned stand was unforeseen. In part, this may reflect that the selected crop trees being monitored were the larger individuals, whereas the mortality occurred in the untagged smaller trees, although it is also not at odds with density- dependent mortality and simply signals that, up to 58 years, mortality of silver beech is low among the wave of ‘crop’ trees. This slow rate of self-thinning does limit timber yield and accounts for the potential need for silvicultural interventions in silver beech stands. Conversely, the implications of thinning leading to mortality of residual trees need not be of concern. Consistent with thinned eucalypts (Kariuki 2008), we found that intense thinning increased mortality in retained trees but the effects were transient and of minor impact. In silver beech the probability of mortality decreased for c. 16 years after a second thinning, but rates were generally low (≤ 0.007 year–1). The low mortality associated with heavy thinning corresponds with findings that patch cuts (<0.2 ha) in southern beech forests do not elevate silver beech tree mortality at the edge of cuts (Wiser et al. 2005). Given that many factors can influence tree growth, it is notable that stand age and density alone explained up to 63 % of the recorded variation in cumulative growth for silver beech at our study site. Yet, accounting for the remaining within-stand variability in growth was vital for realistic predictions of yield (Fig. 5). Simulations that assumed zero or constant growth variability across stand densities produced opposite and markedly biased predictions at high density (Fig. S5 vs S3, Supplemental information). We found larger growth variability at low density, which differs from a common observation that size inequality increases at high densities in even- aged stands (e.g. Nord-Larsen et al. 2006). This can be explained by a confounding effect of tree size. Plant density is only meaningful when related to mean plant size so that the performance of large trees at low density can resemble that of small trees at high density (Weiner et al. 2001). The differences in tree genetics, time of establishment and microsites can take a long time to be manifested in slow growing trees in dense stands, but rapidly become evident in larger trees at low density. Also worthy of attention was a strong correlation (r = 0.72) between DBH at 20 and 58 years. This signals auto-correlated growth, with early establishment or fast initial growth determining later growth because of competitive advantage, as previously documented by tree ring studies (e.g., Brienen et al. 2006) and serves to identify candidate crop trees at an early stage. Stand-level responses Although it seems indisputable that thinning can promote merchantable yield, it is usually unknown what density leads to maximum yield for a given species, site and stand age (Zeide 2001). Initial advice suggested thinning to 200–300 stems ha–1 (Franklin 1981; Wardle 1984) but both empirical evidence (Fig. 5a) and simulations (Fig. 5b and Fig. S3, ESM1) reveal that those densities are far too low and compromise yield due to unused space and, to lesser extent, increased tree mortality. Regressions with empirical data indicate that, up to 58 years, yield peaks at an optimum density of 570 stems ha–1 for silver beech. Simulation results suggest a higher optimum of 850 stems ha–1 for 54- to 64-year- old stands and 950 stems ha–1 for older stands. Both a quantitative study (Pretzsch 2005) and a theoretical review (Zeide 2008) have indicated that while low densities foster diameter growth early on, it is high densities that promote ingrowth of smaller stems into a merchantable DBH later on (see Fig. S2). So even though somewhat counterintuitive, the density for optimal yield is expected to be higher for older stands. The simulations confirm this prediction but indicate that the increase in density is relatively minor and the density for optimal yield is fairly stable across a breadth of stand ages (Fig. S4). The stand age for optimal yield in thinned stands appears promisingly shorter than the c. 135 years estimated for maximum yield in untended silver beech stands (Williams & Chavasse 1951). Since the simulation model was calibrated with the original selection of tagged crop trees in dense stands (and those could have outperformed untagged trees), the simulation optimum ought to be interpreted with some caution. Our simulations of stand development predicted an optimal merchantable yield of 190 m3 ha–1 at 80 years for 950 stems ha–1, which corresponds to clear logs at least 30-cm DBH and merchantable heights not exceeding 8 m. When all trees above 7-cm DBH are instead accounted for (a standard given by Smith et al. 1997) using the same volume function, the yield at 80 years is estimated to be 665 m3 ha–1 for 950 stems ha–1. This generates a mean annual increment of 8.31 m3 ha–1 year–1 which is comparable in magnitude to an earlier estimate of 5.6 to 6.3 m3 ha–1 year–1 at 70–80 years for thinned silver beech on fertile, low elevation sites in Easdale et al. New Zealand Journal of Forestry Science (2022) 52:16 Page 11 southern New Zealand (Valentine 1969; Wardle 1984). When compared with thinned second-growth southern beech forest in Patagonia, our estimated yield of 665 m3 ha–1 at 80 years for 950 stems ha–1 corresponds well with those from other southern beech forests at similar latitudes: 400–600 m3 ha–1 for a 70–80-year rotation of thinned Nothofagus pumilio at latitude 45.5°S (Nuñez & Vera 1992) and 690 m3 ha–1 for 70–75-year- old thinned Lophozonia alpina at c. 42.0°S (Grosse 1989). Not surprisingly, the above values sit well above a total yield of 160 m3 ha–1 for 60-year-old Nothofagus betuloides stands thinned to 1000 stems ha–1 in the southernmost southern beech forests at 54.8°S (Fig. 2 in Martínez-Pastur et al. 2010). Implications for managing New Zealand’s southern beech forest Current legislation largely restricts harvests in southern beech forest to patch cuts not greater than 0.5 ha (New Zealand Ministry of Agriculture and Forestry 2007). We know that recruitment in such cuts is variable but often prolific and usually well exceeds the density of large trees in mature forest (Wiser et al. 2007; Allen et al. 2012). Where recruitment is abundant, subsequent diameter growth tends to be stagnant in unthinned stands (Franklin 1995) and our results demonstrate how thinning can convert more of the regenerating stands into merchantable yield within a shorter time interval. In temperate forests, above-ground woody biomass production tends to be higher (up to c. three-fold higher) in forests on soils with high nutrient availability compared with soils with low nutrient availability (Vicca et al. 2012). Since the Alton Valley has moderately fertile Orthic Brown soils with low to moderate base saturation (Hewitt 2010), we might expect lower yields on soils that are less fertile. In addition, reductions in stand biomass productivity have been demonstrated with elevation for one southern beech species (Harcombe et al. 1998) and likely translate to similar reductions in silver beech. However, we expect that, once the absolute effects of site are considered, the relative effects of density on yield found here will apply to silver beech in other sites. Another consideration is managed species, with indications that responses to thinning likely differ between beech species. Assessment of thinning outcomes in a mixed red-hard beech (Fuscospora fusca (Hook.f.) Heenan & Smissen; Fuscospora truncata (Colenso) Heenan & Smissen) forest at Staircase Creek, Westland, indicated that optimum densities were lower for these species at c. 200 stems ha–1 and merchantable yield potential was substantive, with an estimated maximum 191 m3 ha−1 (stems 30-cm DBH) at 58 years (Easdale et al. 2010). Optimisation of yield is of course only one possible expectation from developing stands which could be managed also for other goals such as carbon sequestration (Evison et al. 2012) and biodiversity (Allen et al. 2012). Silvicultural guidelines often suggest that thinning should be sequential. This aims to gradually liberate space and resources as trees grow without limiting crown responses to later liberations (Clutter et al. 1983; Smith et al. 1997). Sequential thinning also aims to maintain some degree of natural competition that favours tall, clean, straight trunks (Smith et al. 1997) and promotes wind stability (Franklin 1995). Our results suggest a single thinning suffices and its timing, between two interventions, made no essential difference to silver beech growth. The slope of the DBH curve was largely explained by final density after a second thinning with negligible improvements when accounting for density after the first thinning. Noting a relatively flat peak in the yield response to density (Fig. S4), and guided more by empirical than simulation results, we suggest that 600−800 stems ha–1 should maximise yield at different harvest times for a merchantable size of 30 cm DBH. With total yield being greater in dense stands (Zeide 2001), lower merchantable thresholds would likely shift the optimum yield to higher densities, warranting a separate sensitivity analysis on the impacts of merchantable size. In general, the results support the view that the southern beech respond well after periods of suppression (Stewart et al. 1991) when tree stems have a DBH <50 cm (Wiser et al. 2005). This offers flexibility in the timing of interventions, where a single late thinning would ease the identification of crop trees (Franklin 1995), promote self-pruning (Smith et al. 1997) and reduce financial costs (Price 1989). Experience has shown that very late thinning of pole stands results in wind-throw and attack by Platypus beetles (Franklin 1995), so thinning should not be delayed excessively and ought to be completed by the time trees are 15 m tall and 15 cm DBH (Franklin 1995). Conclusions Our decadal-level study of tree- and stand-level responses to thinning showed the utility of remeasuring stands along a controlled and delineated density gradient. Because thinning increased tree growth, but had minimal effect on tree mortality, our results alleviate concerns about the stability and productivity of thinned silver beech stands. The density which optimised yield was about two-fold greater than that previously recommended for silver beech and also gave a seven- fold gain in merchantable yield at 58 years. This density remained relatively stable through stand development and suggested that a single density will be adequate for a wide range of harvest ages, although this should take place before stand age reaches 80 years. This rotation age is much shorter than that previously estimated for maximum silver beech yield in unthinned stands. Such conclusions are not only relevant to managing regeneration within coupes harvested under Part 3A of the Forests Act 1949 but also to any areas being planted with silver beech. Competing interests The authors declare that they have no competing interests. Acknowledgements This trial was established and maintained by the former New Zealand Forest Service and then supported by the Foundation for Research, Science and Technology (now Ministry of Business, Innovation and Employment) [C09X0308], the Ministry of Agriculture and Forestry (now Ministry for Primary Industries) [SUS 8905], and the authors. We are grateful to Lindsay and Dixon Limited for support and access to the study site, Gordon Baker for assistance with remeasurement, to Renske Terhürne for data assembly, to Guy Forrester, Jenny Hurst and Dean Anderson for statistical advice, to Cissy Pan for graphics assistance and to Christine Bezar and Leah Kearns for editorial support. Comments from Susan Wiser, Horacio Bown and anonymous referees helped to improve previous versions of this work. Digital records for the thinning trial are available from New Zealand’s National Vegetation Survey databank (https://nvs. landcareresearch.co.nz). Authors' contributions TAE conceptualised this study, remeasured and interpreted the experiment, curated data, analysed data, and wrote the manuscript. RBA conceptualised this study, remeasured and interpreted the experiment, wrote the manuscript, and obtained funding. DAF conceptualised this study, established the experiment, remeasured and interpreted the experiment, and reviewed the manuscript. LEB remeasured and interpreted the experiment, curated data, and reviewed the manuscript. DH remeasured and interpreted the experiment, curated data, and reviewed the manuscript. References Allen, R.B., Hurst, J.H., Wiser, S.K., & Easdale, T.A. (2012). Developing management systems for the production of beech timber. New Zealand Journal of Forestry, 57(2), 38–44. Baker, G., & Benecke, U. (2001). Growth response to silvicultural tending of red and silver beech regeneration. New Zealand Journal of Forestry, 46(3), 30–32. Brienen, R.J., Zuidema, P.A., & During, H.J. (2006). Autocorrelated growth of tropical forest trees: unravelling patterns and quantifying consequences. Forest Ecology and Management, 237(1-3), 179–190. https://doi.org/10.1016/j. foreco.2006.09.042 Burnham, K.P., & Anderson, D.R. (2002). Model Selection and Multimodel Inference: A Practical Information– Theoretical Approach. New York: Springer. Burrows, L.E., & Allen, R.B. (1991). Silver beech (Nothofagus menziesii (Hook. f.) Oerst.) seedfall patterns in the Takitimu Range, South Island, New Zealand. New Zealand Journal of Botany, 29(4), 361-365. https://doi.org/10.1080/002882 5X.1991.10415489 Clutter, J.L., Fortson, J.C., Pienaar, L.V., Brister, G.H., & Bailey, R.L. (1983). Timber Management: a Quantitative Approach. Malabar: Krieger. Deng, J.M., Ran, J.Z., Wang, Z.Q., Fan, Z., Wang, G., Ji, M., Liu, J., Wang, Y., Liu, J., & Brown, J.H. (2012). Models and tests of optimal density and maximal yield for crop plants. Proceedings of the National Academy of Science, 109(39), 15823–15828. https://doi. org/10.1073/pnas.1210955109 Easdale, T.A., Allen, R.B., Peltzer, D.A., & Hurst, J.J. (2012). Size-dependent growth responses to competition and environment in Nothofagus menziesii. Forest Ecology and Management, 270, 223–231. https:// doi.org/10.1016/j.foreco.2012.01.009 Easdale, T.A., Burrows, L.E., Henley, D., & Baker, G. (2009). Effects of thinning on silver beech growth. [Landcare Research Contract Report, LC0910/015], 24 p. Lincoln, New Zealand: Landcare Research. Easdale, T.A., Drew, K., Henley, D., & Baker, G. (2010). Effects of thinning on growth of red beech and hard beech. [Landcare Research Contract Report, LC0910/179], 24 p. Lincoln, New Zealand: Landcare Research. Ellis, J.C. (1979). Tree volume equations for the major indigenous species in New Zealand. [FRI Bulletin No. 237], Rotorua: New Zealand Forest Service. Evison, D., Easdale, T., Mason, E., & Sewell, A., (2012). Economics of managing New Zealand silver beech for timber and carbon. [NZARES No. 1170-2016- 93282], Nelson: New Zealand Agricultural and Resource Economics Society (Inc.). Farnden, C., & Herring, L. (2002). Severely repressed lodgepole pine responds to thinning and fertilization: 19-year results. Forestry Chronicle, 78(3), 404–414. https://doi.org/10.5558/ tfc78404-3 Franklin, D.A. (1981). A silvicultural regime for dense beech regeneration. [FRI What’s New in Forest Research 98], Rotorua: Forest Research Institute. Franklin, D.A. (1995). Management of beech. In: D. Hammond (Ed.), Forestry Handbook (pp. 95–97). Christchurch: N.Z. Institute of Forestry. Gelman, A., & Hill, J. (2007). Data Analysis Using Regression and Multilevel/Hierarchical Models. New York: CUP. https://doi.org/10.1017/CBO9780511790942 Goffe, W.L., Ferrier, G.D., & Rogers, J. (1994). Global optimization of statistical functions with simulated annealing. Journal of Econometrics, 60(1-2), 65–99. https://doi.org/10.1016/0304-4076(94)90038-8 Grosse, H. (1989). Renovales de raulí, roble, coigüe y tepa. Expectativas de rendimiento. Ciencia e Investigacion Forestal, 3(1), 37–72. https://doi. org/10.52904/0718-4646.1989.110 Easdale et al. New Zealand Journal of Forestry Science (2022) 52:16 Page 12 https://nvs.landcareresearch.co.nz https://nvs.landcareresearch.co.nz https://doi.org/10.1016/j.foreco.2006.09.042 https://doi.org/10.1016/j.foreco.2006.09.042 https://doi.org/10.1080/0028825X.1991.10415489 https://doi.org/10.1080/0028825X.1991.10415489 https://doi.org/10.1073/pnas.1210955109 https://doi.org/10.1073/pnas.1210955109 https://doi.org/10.1016/j.foreco.2012.01.009 https://doi.org/10.1016/j.foreco.2012.01.009 https://doi.org/10.5558/tfc78404-3 https://doi.org/10.5558/tfc78404-3 https://doi.org/10.1017/CBO9780511790942 https://doi.org/10.1016/0304-4076(94)90038-8 https://doi.org/10.52904/0718-4646.1989.110 https://doi.org/10.52904/0718-4646.1989.110 Harcombe, P., Allen, R.B., Wardle, J.A., & Platt, K.H. (1998). Spatial and temporal patterns in structure, biomass, growth, and mortality in a monospecific Nothofagus solandri var. cliffortioides forest in New Zealand. Journal of Sustainable Forestry, 6(3-4), 313–345. https://doi.org/10.1300/ J091v06n03_06 Harrington, C.A., & Reukema, D.L. (1983). Initial shock and long-term stand development following thinning in a Douglas-fir plantation. Forest Science, 29(1), 33–46. Hewitt, A.E. (2010). New Zealand Soil Classification. [Landcare Research Science Series 1, 3rd ed], Lincoln: Manaaki Whenua Press. Husch, B., Beers, T.W., & Kershaw, J.A. (2003). Forest Mensuration. NJ, USA: Wiley. Kariuki, M. (2008). Modelling the impacts of various thinning intensities on tree growth and survival in a mixed species eucalypt forest in central Gippsland, Victoria, Australia. Forest Ecology and Management, 256(12), 2007–2017. https://doi. org/10.1016/j.foreco.2008.07.035 Keating, B.A., Carberry, P.S., Bindraban, P.S., Asseng, S., Meinke, H., & Dixon, J. (2010). Eco-efficient agriculture: concepts, challenges, and opportunities. Crop Science, 50, S109–S119. https://doi.org/10.2135/cropsci2009.10.0594 Langsaeter, A. (1941). Om tynning i enaldret gran-og furuskog. Meddelelser fra det Norske Skogforsøksves, 8, 131–216. Lanner, R.M. (1985). On the insensitivity of height growth to spacing. Forest Ecology and Management, 13(3-4), 143–148. https://doi.org/10.1016/0378- 1127(85)90030-1 Martínez-Pastur, G.J., Lencinas, M.V., Peri, P.L., Cellini, J.M., & Moretto, A. (2010). Long-term forest management research in South Patagonia – Argentina: Lessons from the past, challenges from the present. Revista Chilena de Historia Natural, 83(1), 159–169. https://doi.org/10.4067/S0716- 078X2010000100009 Online report: Ministry for Primary Industries. (2018). Measuring Indigenous Trees and Logs - A Field Guide. Wellington: Ministry for Primary Industrieshttps:// www.mpi.govt.nz/dmsdocument/45-Indigenous- Forestry-Measuring-Indigenous-trees-and-logs-A- Field-Guide.Accessed 7 Aug 2018 Murphy, L. (2012). Likelihood: methods for maximum likelihood estimation. R package version 1.5. h t t p s : / / c r a n . r - p r o j e c t . o r g / w e b / p a c k a g e s / likelihood/ New Zealand Ministry of Agriculture and Forestry. (2007). Standards and Guidelines for the Sustainable Management of Indigenous Forests. Wellington: MAF. Nishizono, T. (2010). Effects of thinning level and site productivity on age-related changes in stand volume growth can be explained by a single rescaled growth curve. Forest Ecology and Management, 259(12), 2276–2291. https://doi.org/10.1016/j. foreco.2010.03.002 Nord-Larsen, T., Damgaard, C., & Weiner, J. (2006). Quantifying size-asymmetric growth among individual beech trees. Canadian Journal of Forest Research, 36(2), 418–425. https://doi. org/10.1139/x05-255 Nuñez, P., & Vera, O. (1992). Evaluación de intervenciones silvícolas en un renoval mixto de Lenga (Nothofagus pumilio) y Coihue (Nothofagus dombeyi), ubicado en la reserva forestal Coyhaique, XI Región. In Actas del Seminario de Manejo Forestal de la Lenga y Aspectos Ecológicos Relacionados (pp. 111–125). Esquel: CIEFAP. Osawa, A., & Allen, R.B. (1993). Allometric theory explains self-thinning relationships of mountain beech and red pine. Ecology, 74(4), 1020–1032. https://doi.org/10.2307/1940472 Pretzsch, H. (2005). Stand density and growth of Norway spruce (Picea abies (L.) Karst.) and European beech (Fagus sylvatica L.): evidence from long- term experimental plots. European Journal of Forest Research, 124(3), 193–205. https://doi. org/10.1007/s10342-005-0068-4 Price, C. (1989). The Theory and Application of Forest Economics. Oxford: Blackwell. R Development Core Team. (2013). R: A Language and Environment for Statistical Computing. Vienne: R Foundation for Statistical Computing. http:// www.R-project.org. Richardson, S.J., Hurst, J.M., Easdale, T.A., Wiser, S.K., Griffiths, A.D., & Allen, R.B. (2011). Diameter growth rates of beech (Nothofagus) trees around New Zealand. New Zealand Journal of Forestry, 56(1), 3–11. Sit, V., & Poulin-Costello, M. (1994). Catalog of curves for curve fitting. [Biometrics Information, Handbook No. 4], Victoria: BC Ministry of Forests. Skovsgaard, J.P., & Vanclay, J.K. (2007). Forest site productivity: a review of the evolution of dendrometric concepts for even-aged stands. Forestry, 81(1), 13–31. https://doi.org/10.1093/ forestry/cpm041 Smith, D.M., Larson, B.C., Kelty, M.J., & Ashton, P.M.S. (1997). The Practice of Silviculture: Applied Forest Ecology. New York: John Wiley. Stewart, G.H., Rose, A.B., & Veblen, T.T. (1991). Forest development in canopy gaps in old-growth beech (Nothofagus) forests, New Zealand. Journal of Vegetation Science, 2(5), 679–690. https://doi. org/10.2307/3236178 Easdale et al. New Zealand Journal of Forestry Science (2022) 52:16 Page 13 https://doi.org/10.1300/J091v06n03_06 https://doi.org/10.1300/J091v06n03_06 https://doi.org/10.1016/j.foreco.2008.07.035 https://doi.org/10.1016/j.foreco.2008.07.035 https://doi.org/10.2135/cropsci2009.10.0594 https://doi.org/10.1016/0378-1127(85)90030-1 https://doi.org/10.1016/0378-1127(85)90030-1 https://doi.org/10.4067/S0716-078X2010000100009 https://doi.org/10.4067/S0716-078X2010000100009 https://www.mpi.govt.nz/dmsdocument/45-Indigenous-Forestry-Measuring-Indigenous-trees-and-logs-A-Field-Guide https://www.mpi.govt.nz/dmsdocument/45-Indigenous-Forestry-Measuring-Indigenous-trees-and-logs-A-Field-Guide https://www.mpi.govt.nz/dmsdocument/45-Indigenous-Forestry-Measuring-Indigenous-trees-and-logs-A-Field-Guide https://www.mpi.govt.nz/dmsdocument/45-Indigenous-Forestry-Measuring-Indigenous-trees-and-logs-A-Field-Guide https://cran.r-project.org/web/packages/likelihood/ https://cran.r-project.org/web/packages/likelihood/ https://doi.org/10.1016/j.foreco.2010.03.002 https://doi.org/10.1016/j.foreco.2010.03.002 https://doi.org/10.1139/x05-255 https://doi.org/10.1139/x05-255 https://doi.org/10.2307/1940472 https://doi.org/10.1007/s10342-005-0068-4 https://doi.org/10.1007/s10342-005-0068-4 http://www.R-project.org http://www.R-project.org https://doi.org/10.1093/forestry/cpm041 https://doi.org/10.1093/forestry/cpm041 https://doi.org/10.2307/3236178 https://doi.org/10.2307/3236178 Trewavas, A.J. (2001). The population/biodiversity paradox. Agricultural efficiency to save wilderness. Plant Physiology, 125(1), 174–179. https://doi. org/10.1104/pp.125.1.174 Valentine, J.M. (1969). Some observations on the management of silver beech. New Zealand Journal of Forestry, 14(2), 219–226. Vicca, S., Luyssaert, S., Penuelas, J., Campioli, M., Chapin, F.S., Ciais, P., Heinemeyer, A., Hogberg, P., Kutsch, W.L., Law, B.E., Malhi, Y., Papale, D., Piao, S.L., Reichstein, M., Schulze, E.D., & Janssens, I.A. (2012). Fertile forests produce biomass more efficiently. Ecology Letters, 15(6), 520–526. https://doi. org/10.1111/j.1461-0248.2012.01775.x Wardle, J.A. (1984). The New Zealand Beeches. Wellington: New Zealand Forest Service. Weiner, J. (2003). Ecology – the science of agriculture in the 21st century. Journal of Agricultural Science, 141(3-4), 371–377. https://doi.org/10.1017/ S0021859603003605 Weiner, J., & Freckleton, R.P. (2010). Constant final yield. Annual Review of Ecology, Evolution and Systematics, 41, 173–192. https://doi.org/10.1146/annurev- ecolsys-102209-144642 Weiner, J., Stoll, P., Muller-Landau, H., & Jasentuliyana, A. (2001). The effects of density, spatial pattern, and competitive symmetry on size variation in simulated plant populations. The American Naturalist, 158(4), 438–450. https://doi. org/10.1086/321988 Williams, R.W.M., & Chavasse, C.G.R. (1951). The silviculture of silver beech in Southland. New Zealand Journal of Forestry, 6(3), 219–235. Wiser, S.K., Allen, R.B., Benecke, U., Baker, G., & Peltzer, D.A. (2005). Tree growth and mortality after small-group harvesting in New Zealand old- growth Nothofagus forests. Canadian Journal of Forest Research, 35(10), 2323–2331. https://doi. org/10.1139/x05-158 Wiser, S.K., Baker, G., & Benecke, U. (2007). Regeneration of red and silver beech: how important is the size of harvested area? New Zealand Journal of Forestry, 52(2), 31–35. Xue, L., & Hagihara, A. (1998). Growth analysis of the self-thinning stands of Pinus densiflora Sieb. et Zucc. Ecological Research, 13(2), 183–191. https:// doi.org/10.1046/j.1440-1703.1998.00256.x Zeide, B. (2001). Thinning and growth: a full turnaround. Journal of Forestry, 99(1), 20–25. Zeide, B. (2008). The science of forestry. Journal of Sustainable Forestry, 27(4), 345–473. https://doi. org/10.1080/10549810802339225 Easdale et al. New Zealand Journal of Forestry Science (2022) 52:16 Page 14 https://doi.org/10.1104/pp.125.1.174 https://doi.org/10.1104/pp.125.1.174 https://doi.org/10.1111/j.1461-0248.2012.01775.x https://doi.org/10.1111/j.1461-0248.2012.01775.x https://doi.org/10.1017/S0021859603003605 https://doi.org/10.1017/S0021859603003605 https://doi.org/10.1146/annurev-ecolsys-102209-144642 https://doi.org/10.1146/annurev-ecolsys-102209-144642 https://doi.org/10.1086/321988 https://doi.org/10.1086/321988 https://doi.org/10.1139/x05-158 https://doi.org/10.1139/x05-158 https://doi.org/10.1046/j.1440-1703.1998.00256.x https://doi.org/10.1046/j.1440-1703.1998.00256.x https://doi.org/10.1080/10549810802339225 https://doi.org/10.1080/10549810802339225 Supplementary Information This section presents a graphical examination of mortality rates (Fig. S1), a growth simulation example (Fig. S2), empirically calibrated simulation results (Fig. S3 and S4) and simulation outputs with alternative settings for growth variability (Fig. S5). Figures are presented in the same order as topics covered in the main text. FIGURE S1: Empirical and modelled probabilities of mortality of crop trees as a function of density at different time intervals after the second thinning (up to 29 years). Modelled probabilities (Equation 2 in Table 1, main text) were plotted for the midpoint of each interval (2.5, 10, 17 and 24 years after second thinning respectively). 0−5 years after second thinning ●●●●●● ●●●●●● ●●●● ●● ●● ●●●● ●● ●● 200 500 1000 2000 5000 10000 0. 00 0 0. 00 4 0. 00 8 Density at second thinning (stems ha−1) A nn ua l p ro ba bi lit y of m or ta lit y ●● Empirical Modelled 5−15 years after second thinning ●● ●● ●● ●● ●●●● ●●●● ●● ●●●● ●● ●● ●● 200 500 1000 2000 5000 10000 0. 00 0 0. 00 4 0. 00 8 Density at second thinning (stems ha−1) A nn ua l p ro ba bi lit y of m or ta lit y 15−19 years after second thinning ●● ●● ●● ●● ●● ●● ●●●● ●● ●● ●●●● ●● ●● 200 500 1000 2000 5000 10000 0. 00 0 0. 00 4 0. 00 8 Density at second thinning (stems ha−1) A nn ua l p ro ba bi lit y of m or ta lit y 19−29 years after second thinning ●● ●● ●● ●●●●●● ●●●● ●●●● ●●●● ●● ●● 200 500 1000 2000 5000 10000 0. 00 0 0. 00 4 0. 00 8 Density at second thinning (stems ha−1) A nn ua l p ro ba bi lit y of m or ta lit y Easdale et al. New Zealand Journal of Forestry Science (2022) 52:16 Page 15 0 20 40 60 80 0 20 40 60 80 100 stems ha−1 Stand age (years) D B H (c m ) 0 6 12 0 20 40 60 80 N 0 20 40 60 80 0 20 40 60 80 1000 stems ha−1 Stand age (years) 0 15 30 0 20 40 60 80 N 0 20 40 60 80 0 20 40 60 80 10000 stems ha−1 Stand age (years) 0 20 0 20 40 60 80 N FIGURE S2: Simulated diameter (DBH) growth curves for three different densities and their corresponding DBH distributions at 80 years. For comparability, only 100 trees are shown in each case. The dotted line presents the minimum merchantable DBH size. Note the greater variability of growth at low density. y S ta nd a ge (y ea rs ) 15 30 45 60 75 9 0 105 120 135 150 1 65 180 195 12 0 10 0 80 60 40 100 200 500 1000 2000 5000 10000 Density at 29 years (stems ha−1) Recorded yield Simulations 200 500 1000 2000 5000 0 50 10 0 15 0 20 0 Density at 29 years (stems ha−1) M er ch an ta bl e yi el d ( m 3 ha −1 ) 40 years 50 years 60 years 70 years 80 years 90 years 550 750 600 850 650 1000 650 1200 700 1200 700 1200 FIGURE S3: Simulation results, presented as a contour plot of merchantable yield (m3 ha–1) versus density and stand age (mean of 100 simulations). Both contour lines and grey saturation of background follow the ‘topography’ of merchantable yield, with darker grey denoting higher yields. The dotted line presents the optimal densities for different stand ages. FIGURE S4: Simulation results presented as decadal changes in merchantable yield (m3 ha–1) across a gradient of stem densities (mean of 100 simulations). Bold lines and associated figures indicate yields and density ranges spanning >98% of the peak yield at each age. Easdale et al. New Zealand Journal of Forestry Science (2022) 52:16 Page 16 FIGURE S5: Contour plot of merchantable yield for simulations that assume constant variability in tree growth across densities (a; mean of 100 simulations) and the corresponding contrast between recorded yield (points) and simulated yield (lines) at 48 years (b; open circles and thin lines) and at 58 years (b; filled circles and thick lines). Note that predictions greatly exceed observations at high density. Contour plot of merchantable timber yield for simulations that do not incorporate variability in tree growth (c) and corresponding contrast between recorded and simulated yield at 48 and 58 years (d). Note that predictions have a sharp boundary and equate to zero when predicted diameters are smaller than the merchantable threshold (30 cm). S ta nd a ge (y ea rs ) 40 60 8 0 10 0 1 20 140 1 60 1 80 20 0 2 20 240 12 0 10 0 80 60 40 100 200 500 1000 2000 5000 10000 Density at 29 years (stems ha−1) (a) ● ● ● ● ● ● ● ● ● ● ● ● ● 100 200 500 1000 2000 5000 10000 0 50 10 0 15 0 Density at 29 years (stems ha−1) M er ch an ta bl e yi el d ( m 3 ha −1 ) ●● ●● ●● ●● ●● ●● ●● ●● ●● ●● ●●●● ●● (b) 48 yrs 58 yrs ● ●● R ec or d. S im ul . S ta nd a ge (y ea rs ) 20 40 6 0 8 0 1 00 1 20 1 40 1 60 1 80 2 00 12 0 10 0 80 60 40 100 200 500 1000 2000 5000 10000 Density at 29 years (stems ha−1) (c) ● ● ● ● ● ● ● ● ● ● ● ● ● 100 200 500 1000 2000 5000 10000 0 50 10 0 15 0 Density at 29 years (stems ha−1) M er ch an ta bl e yi el d ( m 3 ha −1 ) ●● ●● ●● ●● ●● ●● ●● ●● ●● ●● ●●●● ●● (d) 48 yrs 58 yrs ● ●● R ec or d. S im ul . Easdale et al. New Zealand Journal of Forestry Science (2022) 52:16 Page 17