The relationship between volume increment and stand density in Norway spruce plantations

An understanding of the relationship between volume increment and stand density (basal area, stand density index, etc.) is of utmost importance for properly managing stand density to achieve speciﬁc management objectives. There are two main approaches to analyse growth–density relationships. The ﬁrst relates volume increment to stand density through a basic relationship, which can vary with site productivity, age, and potentially incorporates treatment effects. The second is to relate the volume increment and density of thinned experimental plots relative to that of an unthinned experimental plot on the same site. Using a dataset of 229 thinned and unthinned experimental plots of Norway spruce, a growth model is developed describing the relationship between gross or net volume increment and basal area. The models indicate that gross volume increases with increasing basal area up to 50 m 2 and thereafter becomes constant out to the maximum basal area. Alternatively, net volume increment was maximized at a basal area of 43 m 2 and decreased with further increases in basal area. However, the models indicated a wide range where net volume increment was essentially constant, varying by less than 1 m 3 ha − 1 year − 1 . An analysis of different thinning scenarios indicated that the relative relationship between volume increment and stand density was dynamic and changed over the course of a rotation.


Introduction
In production forests, mid-rotational thinning is implemented for the purpose of reducing natural density-related mortality and to release crop trees from competition, encouraging individual-tree growth and potentially producing higher valued timber products (Nyland, 2016). However, overly reducing the stand density can result in non-optimal total stand volume production. Physiologically, this loss in volume production can be explained as result of underutilizing available resources, such as light interception in the canopy and water and mineral absorption in the soil (Long et al., 2004). At the opposite extreme, excessive densities may also lead to non-optimal volume production if available site resources (e.g. light, water) become a limiting factor. Thus, there is a direct relationship between stand density and volume production. A quantification of this relationship is important for understanding the growth and economic implications of thinning treatments.
While long-term studies of thinning have been on-going since 1860, a general relationship between volume increment and stand density has not yet been determined. One of the issues contributing to the lack of understanding is the wide range of definitions of volume 'growth' and stand 'density' (Zeide, 2001). In terms of total stand volume, volume growth is defined as either net-growth, which includes standing and harvested tree volume, or gross-growth, which is net-growth plus volume of trees lost to mortality (Avery and Burkhart, 2015). Because gross volume growth includes mortality, differences in the relationships between gross or net volume growth and stand density are expected to occur after density-related mortality is initiated (Nyland, 2016). Further differences in volume growth are expected to occur between total and merchantable volume (Zeide, 2001). However, merchantability limits change with different product classes and largely depend on local markets. In this work, only total gross or net volume growth is considered.
Historically, the measure of stand density used to examine growth-density relationships has either been basal area, standing volume, or stand density index. As an estimator of gross volume growth, all three measures have been shown to produce the same growth-density relationship (Allen II and Burkhart, 2018). However, it is important to define how they are calculated for the growth period. For example, Pretzsch (2005) used the stand density index at the beginning of the growth period to examine growth-density relationships in Norway spruce and European beech stands in Germany. The use of basal area at the beginning of the growth period has mostly been used for the development of growth models to describe growth-density relationships and when the growth periods are relatively short (e.g. Zeide, 2004;Gizachew and Brunner, 2011;Allen II and Burkhart, 2018). Alternatively, Mäkinen and Isomäki (2004a) compared volume growth Forestry in thinning experiments in Norway spruce stands in Finland based on the average basal area over a 28-year period. This approach is often used in the analysis of designed experiments and/or when the growth periods are relatively long (e.g. Isomäki, 2004a, 2004b;Nilsson et al., 2010).
Analyses of growth-density relationships have been conducted using two approaches. The first approach is used when comparative thinning experiments, containing an unthinned control plot and one or more thinning treatments, are installed in the same stand. Within this approach, the density and volume growth of thinned plots can be calculated relatively (as a percentage) of the measures observed on the control plots. Often termed 'relative growth-density relationships', examining the relationship between growth and density in this manner allows for the quantification of the percentage increase or decrease in volume increment in the thinned plot as compared with the unthinned plot. The second approach relates absolute measures of density directly to absolute measures of volume growth, termed 'absolute growth-density relationships', and is often used when developing growth models and/or when data from comparative plots are unavailable.
Three hypothesized patterns concerning relative growthdensity relationships have been identified (Zeide, 2001). The first is the constant-optimum pattern of Wiedemann (1932), which suggests that the relative volume growth (RVG) of thinned stands is constant and optimal across a range of absolute basal areas and any increase or decrease outside of this range results in a decrease in RVG (Figure 1). The second is the optimal pattern of Assmann (1950), which suggests that there is a relative basal area (RBA), where RVG is maximized. In young stands, beginning at the basal area of the unthinned or lightly thinned control, RVG increases with a reduction in RBA up to an optimum and then decreases with further reductions in RBA. However, with increasing age the optimum RBA moves towards the unthinned control. Finally, the third pattern is the constant pattern of Mar:Moller (1954), which suggests that RVG is optimal beginning at the basal area (or volume) of the unthinned stand and is constant down to a 50 per cent reduction in basal area (or volume).
Alternatively, there is only one hypothesized pattern identified in the literature for absolute growth-density relationships; the constant-optimum pattern of Langsaeter (1941). Langsaeter (1941) suggested that volume increment increases with increasing standing volume up to a point where volume increment is maximized and constant across a range of standing volumes. At excessive densities volume increment should decline. The limits of standing volume where volume increment is constant and optimal are said to be dependent on species, site and age. This constant-optimum growth-density pattern is similar to that proposed by Wiedemann (1932) with the difference that Langsaeter considered volume growth in absolute units and uses standing volume as the measure of density.
Langsaeter's alternative approach to analysing the relationship between volume increment and stand density is quite interesting when compared with the works of Wiedemann (1932), Assmann (1950) and Mar:Moller (1954). If Langsaeter's hypothesized absolute growth-density pattern is correct, and due to the others defining maximum density as the density of the unthinned (or lightly thinned) control, then constant, optimum and constant-optimum growth-density patterns could all be observed in both absolute and relative growth-density relationships over the course of a rotation. For example, Figure 2 shows that given the same hypothetical constant-optimum absolute growth-density relationship, all three relationships could be observed in both absolute and relative densities depending on the density of the control plot and the absolute or relative densities of the thinned plots. Thus, the relative growthdensity relationship may change over the rotation period as the unthinned control moves through the absolute growth-density relationship. However, this potential linkage between absolute and relative growth-density relationships over a rotation is yet untested in the literature.
These alternative growth-density hypotheses have been tested in a variety of stand types around the world with various conclusions. Perhaps some ambiguity has arisen because there was no indication as to whether gross or net volume increment was considered in the previously mentioned absolute and relative growth-density hypotheses. In a review of the history and evolution of concepts in growth-density relationships, Zeide (2001) indicates that those hypotheses are for gross volume growth. However, more recent examinations of absolute growthdensity relationships have shown that gross volume increment increases with increasing stand density (e.g. Curtis and Marshall, 2009;Gizachew and Brunner, 2011;Allen II and Burkhart, 2018). Given that mortality will reduce net growth below gross growth it is possible that optimum or constant growth-density patterns may emerge in net growth at higher densities. Thus, an examination of growth-density relationships in both types of volume increment is warranted.
Analyses of thinning trials of Norway spruce stands in Fennoscandia have had varying results. Using data from Norway spruce trials in Finland, Mäkinen and Isomäki (2004a) determined that only early heavy thinning slightly, but significantly, reduced gross volume increment over a 28-year period. Based on data from Norway spruce thinning trials in Sweden, Nilsson et al. (2010) determined that medium and heavy thinning also slightly, but significantly, reduced gross volume increment over a 30year period. In contrast to the results from Finland, Nilsson et al. (2010) found that four light thinnings significantly increased net volume growth over the unthinned control and heavily thinned treatments in the same period. However, these results may not be directly comparable due with pre-thinning differences in the stands, different thinning treatments and different site productivities. Further, these two studies only examined relative growth-density relationships based on 28-30-year averages of volume increment and basal area which prevents the detection of differences in the growth-density relationships over time. A detailed examination of absolute and relative growth-density relationships for Norway spruce in Fennoscandia has yet to be conducted on the basis of long-term thinning trials.
The purpose of this study was to test theoretical relationships between gross or net increment and stand density in evenaged Norway spruce stands. The primary research goal was to quantify the absolute relationship between volume growth and stand density and to determine how different density removals affect the relative growth-density relationship between thinned and unthinned stands. An additional goal was to determine the effects of thinning on cumulative total gross and net volume The relationship between volume increment and stand density in Norway spruce plantations production over a general rotation period. To address the research goals the following hypotheses were tested: (1) absolute gross or net volume increment increases with increasing basal area; (ii) relative volume increment changes over the course of a rotation depending on the basal area of the unthinned stand and (iii) cumulative gross and net volume over a rotation is highest in the unthinned stand.

Methods
A dataset of long-term thinning experiments of Norway spruce was used to develop a modelling framework by which the alternative thinning hypotheses could be tested. This modelling framework includes individual functions for stand-level gross or net volume increment, dominant height and basal area. Additionally, a Monte-Carlo analysis was performed for nine different thinning scenarios to compare cumulative gross and net volume production at harvesting age from different thinning treatments.

Thinning trial data
The data used in this study come from a series of 11 thinning trials consisting of 135 long-term permanent plots established in even-aged Norway spruce plantations. Within those data are 442 observations (growth periods) in once thinned stands and 246 observations in twice thinned stands (Table 1). The intent of establishing the trials was to examine the effects of various levels of thinning intensity and timing of thinning, using dominant stand height as the thinning trigger, on volume production and stand structure. In all cases, thinning was performed from below. Initially, the trials were planned to include treatments of different density reductions at different stages of stand development. Over the course of the experiment, the target residual TPH were often missed (thinning to light or to heavy) and/or the timing of thinning was performed to early or late. The deviation from the original experimental plan was such that formal analyses of the trial data by treatment groups could not be performed. However, this deviation has resulted in a dataset which covers a wide range of thinning intensities and is highly useful for the development of models which explain growth response to thinning.
The thinning trials are distributed across the natural range of Norway spruce in Norway in the counties of Agder, Vestfold, Telemark, Viken, Innlandet and Trøndelag. A full description of these trials is provided by Braastad and Tveite (2001). Within this thinning trial series, there were an additional 65 plots which did not have information concerning the trees removed in the first thinning. Observations from those plots were not reported or used for analysis in this work. As the thinning trial data did  A = stand age from planting (years); SI = site index at base age 40 years; TPH = number of trees (ha −1 ); BPA = basal area (m 2 ha −1 ); V = stand volume (m 3 ha −1 ); TQ BA = basal area thinning quotient (BA after thinning divided by basal area before thinning, %); and H T = dominant stand height (m) at thinning.
not include unthinned control plots, an additional 94 unthinned control plots from other silvicultural trials, within the same geographic region, were added and consist of 311 additional growth observations (Table 1).

Tree measurements and stand calculations
Establishment of the trials occurred between 1969 and 1976 in 20-31-year-old stands (Table 1). Measurements were scheduled in 5-year intervals, with a few deviations from the schedule, and The relationship between volume increment and stand density in Norway spruce plantations to date each plot has received between 5 and 10 measurements. At plot establishment, every tree greater than 5 cm in diameter at breast height (dbh) was numbered and the species, dbh, and any damages were recorded. Total tree heights and crown heights were subsampled by randomly selecting one tree out of the first four trees for a height measurement, and then systematically sampling the height of every fourth tree thereafter. Non-sampled tree heights were estimated with a linear mixed-effects variant of the logarithmic height-diameter model of Lenhart (1968) localized with random parameters fit to each trial, plot and measurement period. Total stem volumes were calculated using the volume equations of Vestjordet (1967). Ingrowth trees were included in plot measurements once reaching a dbh of 5 cm. Quadratic mean diameter (QMD, cm), standing volume (V, m 3 ha −1 ) and basal area (BA, m 2 ha −1 ) were calculated for each plot. Net stand-level volume (NV, m 3 ha −1 ) was calculated as the sum of V and the cumulative volume of all trees removed in thinning. Gross stand-level volume (GV, m 3 ha −1 ) was calculated as NV plus the cumulative volume of any trees lost to mortality. Stand-level gross (GVPAI) and net (NVPAI) volume increment were calculated as the periodic annual increment over a measurement period.

Modelling and analyses
To address the research goals, three hypotheses were tested. The first was to test the relationship between absolute gross or net volume growth and absolute basal area. This was done by developing an equation explaining periodic annual volume increment as a function stand density, dominant stand height and dominant height increment. The second hypothesis was to examine if different relative growth-density relationships exist over the course of a rotation. As the equation for gross or net volume increment included both basal area and dominant height, additional equations for those variables were also developed. These equations were then used to simulate stand development over a rotation for nine different thinning scenarios (Table 2), including an unthinned control. The third and final hypothesis was to test differences in cumulative gross and net volume production over a rotation among the nine thinning scenarios. Based on the simulations, a Monte-Carlo analysis was performed to determine differences among the thinning scenarios.

Gross and net volume growth
Alternative growth-density hypotheses suggest that volume increment either increases with increasing stand density or reaches an upper limit at some density less than the maximum density which can be sustained on a given site. In order to test these hypotheses, the type II combination power and exponential function was chosen as it has the flexibility to reflect increasing, optimum, asymptotic patterns (Allen II and Burkhart, 2018). The functional form relating stand density to volume increment can be expressed as: where VPAI = gross or net volume periodic annual increment (m 3 ha −1 year −1 ), BA = stand basal area (m 2 ha −1 ), exp = the base of the natural logarithm and α 0 , α 1 = are the parameters to be estimated. From equation (1), the hypothesis that volume increment follows an increasing pattern with stand density can be tested on the basis of the significance of α 1 . If α 1 is significantly different from zero and positive, then an optimum pattern occurs and otherwise an increasing pattern occurs. In the case that α 1 is significantly positive, the basal area at which VPAI is maximized can be determined by calculating the first derivative of equation (1) and solving for the basal area where the first derivative is equal to zero. It follows that the basal area where PAI is maximized is equal to α 1 . Equation (1) expresses volume growth as purely a function of stand density. However, in addition to stand density, volume growth is also a function of average tree size, stand age and site productivity (Oliver and Larson, 1996;Zeide, 2004;Allen II and Burkhart, 2018). Therefore, the development of a model which can account for both the effects of stand development stage and soil/climate productivity is necessary to prohibit confounding those effects.
Predominant stand density measures, such as basal area, couple stem number per unit area with a measure of average tree size which is generally diameter-based. It has been shown that basal area can vary in equally dense stands and therefore the inclusion of a measure of average tree size assists in defining a unique growing condition (Zeide, 2005). From past studies of growth-density relationships, a positive parameter estimate on average tree size suggests that for a given basal area a larger average tree size results in greater stand volume growth (e.g. Zeide, 2004;Pretzsch, 2005;Allen II and Burkhart, 2018). However, this does not reflect that with increasing time, and thus increasing average tree size, stand volume growth declines (Avery and Burkhart, 2015) and an additional variable which reflects the gradual decline of volume growth with time is needed. Using the average tree height as a measure of average tree size in combination with its increment reflects this gradual decline in volume growth as height increment declines with time. Further, these two variables together describe both site productivity and local climate impacts over the growth period as well as age effects.
Two different definitions of tree height were examined including the average height of all trees and the average height of the 100 largest diameter trees, also known as dominant stand height. The dominant stand height definition provided slightly better fit statistics. Therefore, the dominant stand height (H, m), its increment (H PAI , m year −1 ) and their interaction were included in equation (1) as modifiers to the parameterα 0 . The final model has the form: In the model-development process, additional variables describing thinning were also included in equation (2). However, these variables were never significant in the model and therefore excluded. Additionally, stand age and site index were included as modifiers to parameter α 1 , which defines the basal area at which VPAI is maximized. These parameters were also not significant Forestry in the model and when removed no evidence of patterns in the residuals from the final equation were seen in stand age or site index.
Equation (2) can be used to determine the gross or net volume increment at a given basal area, dominant stand height and dominant stand height increment. It was additionally desired to understand how both types of volume increment develop over a rotation for both thinned and unthinned stands. This was accomplished by simulating gross and net volume increment over a traditional rotation period for different thinning scenarios. As equation (2) explains volume increment as a function of basal area and height, in order to perform this simulation additional equations for basal area and height development were needed.

Basal area development and response to thinning
In this analysis, it was desired to simulate basal area development of an unthinned stand based on a given initial stand density. However, no information on planting density was available for the unthinned control plots (Table 1). In order to simulate the basal area development in unthinned stands the model of Gizachew et al. (2012) was used. Their model is based on empirical results from Norway spruce spacing trials in Norway and was developed to explain basal area development of different initial planting spacings. To determine the validity of this equation, the first available density observation from each unthinned plot which had not yet entered into density-dependent mortality was used as a proxy for initial planting density. Based on residual analysis, the equation was found to be unbiased for the unthinned control plots and therefore acceptable to use for simulations.
To model the basal area of thinned stands, the basal area equation of Gizachew et al. (2012) was modified to account for the effects of thinning on radial growth. This included parameters which could explain differences in thinning removals and the timing of thinning. This model has the form: 11.715 , β 2 = log(1 − BA AT /β 0 ), β 1i = parameters to be estimated, SI = site index (m) at base age 40, TQ BA = basal area after thinning/basal area before thinning, BA AT = basal area after thinning (m 2 ha −1 ) and yst = years since thinning.
The asymptotic limit of basal area, parameter β 0 , was developed by Gizachew et al. (2012) from Norwegian forest inventory data. Thinning response is incorporated into the rate parameter, β 1 , as a linear function of basal area removal and site index. Parameter β 2 conditions the function such that the basal area at the thinning age (yst = 0) is equal to the basal area immediately after thinning. An additional variable was initially included to account for the effect of dominant stand height at the time of thinning on basal area development. This parameter was not significant in the model and therefore removed. Further, no obvious bias was seen in the residuals when examined over dominant stand height at the time of thinning.
For simulating volume growth in equation (2), an estimate of dominant height and dominant height increment is needed. Therefore, a dominant height projection model was developed using the Chapman-Richards type function: where H 1 = current dominant stand height at age A 1 , H 2 = future dominant stand height at age A 2 and θ 1 , θ 2 = are the parameters to be estimated. Previous results from thinning studies Norway spruce in other regions have indicated that thinning has no effect on dominate height development (Mäkinen and Isomäki, 2004a). In an initial screening of dominant height-age equations for these data, preliminary residual analysis indicated no patterns among the residuals for thinned and unthinned stands. Additionally, no variables describing thinning effects were significant when included in equation (4).
The relationship between volume increment and stand density in Norway spruce plantations

Model fitting and evaluation
Equations (2)-(4) were simultaneously fitted by nonlinear seemingly unrelated regression (NSUR), using the SAS/ETS MODEL Procedure (SAS Institute Inc., 2011). NSUR takes into account the contemporaneous correlations among nonlinear regression equations (Borders, 1989;Robinson, 2004). Model assumptions of normality and variance homogeneity of the residuals for all fitted equations were evaluated visually by examining normal Q-Q plots and residual plots, respectively. Further, residual plots in relation to age, site index and stand treatment (thinned vs. unthinned) were examined to determine if these stand variables were properly represented by each model.

Simulation and Monte-Carlo analysis
For the comparison of different thinning scenarios, a simulation was performed using equations (2)-(4). Different site productivities were evaluated based on site index classes of 11, 14, 17 and 20 m at a base age of 40 years, where age is defined as years since planting. A total of nine different thinning scenarios, including an unthinned control, were simulated for each site index class (Table 2). Simulated thinnings were performed based on the representative thinnings in the data used for model fitting (Table 1) and can be categorized as either medium or heavy thinning, based on percentage basal area removal, and early or late thinning based on dominant stand height. The lengths of the simulations varied by site index class and were based on general rotation ages of even-aged Norway spruce stands. These ages were 120, 96, 81 and 77 years for site indices 11, 14, 17 and 20 m, respectively, and come from analysis of the culmination of mean annual increment in Norway spruce stands (Søgaard et al., 2019).
The steps for simulating the nine different thinning scenarios presented in Table 2 are as follows: (1) For each site index class simulate dominant height development using equation (4) for every age over the rotation using the respective site index value for H 1 and 40 years for A 1 .
(2) Simulate basal area development for the unthinned scenario using the equation of Gizachew et al. (2012) for an initial planting density of 4000 trees ha −1 . (3) Simulate basal area development for the thinned scenarios for one or two thinnings: (a) First thinning (i) Using the simulated dominant heights and basal areas from the unthinned scenario, reduce the basal area in each thinning scenario by its respective percentage at its respective dominant height as detailed in Table 2. (4) Calculate gross or net volume increment based on simulated dominant heights, dominant height increments and basal areas for all nine scenarios.
Initially, it was desired to include a comparison of different initial stand densities. However, the result of the hypotheses examined, and the conclusions were the same and did not provide additional insight into the examination of growth-density hypothesis. Thus, one initial planting density of 4000 trees per hectare was chosen.
A Monte-Carlo analysis was performed to test significant differences in cumulative gross and net volume production among thinning treatments over the rotation. In this analysis, 1000 datasets were generated by randomly selecting 999 observations from the fitting dataset with replacement. For each dataset, equations (2)-(4) were fitted using NSUR creating 1000 unique sets of parameter estimates. Those parameters estimates were then used in the simulation of the nine thinning scenarios to create a distribution of cumulative gross and net volume at harvesting age for the four site indices. Differences among the thinning treatments were then determined by calculating 95 per cent confidence bands from the 2.5 and 97.5 percentiles from the distribution.

Visual assessment of relative growth-density relationships following simulation
To compare the implied relative growth-density relationships from the simulations, within each site index class the gross or net volume increment and basal area of the simulated thinned stands are related to those of the unthinned stand. The relative gross (RGVI) or relative net (RNVI) volume increments and RBA can be calculated as: where U and T denote the unthinned and thinned gross volume increment or basal area, respectively. Figures where then generated relating relative volume increment to RBA for each site index class at different stand ages.

Results
Parameter estimates and fit statistics for equations (2)-(4) are presented in Table 3. Analysis of residuals and normal Q-Q plots indicated that the model assumptions were met. Further, there was no evidence of bias among the residuals when plotted against age, site index, treatment types and climatic variables.

Absolute growth-density relationships
Hypotheses tests concerning the growth-density patterns for gross and net volume increment were made based on parameter α 1 in equation (2). Those parameter estimates were significantly greater than zero, indicating optimum patterns between gross Forestry or net volume increment and basal area. Further, the values of those parameters indicate that gross or net volume increment for Norway spruce stands are maximized at basal areas of about 62 and 43 m 2 ha −1 , respectively. To visually assess the behaviour of equation (2) when all other factors are held constant, predicted gross and net increments were plotted across the range of expected basal areas for Norway spruce and with a range of dominate stand heights and height increments (Figure 3). While equation (2) indicates that gross volume increment is maximized at 62 m 2 ha −1 basal area, Figure 3 shows that when height and height increment are held constant, GV PAI is almost constant at basal areas greater than 50 m 2 ha −1 . Similarly, NV PAI has only a slight variation over the basal area range of 30-50 m 2 ha −1 but declines at higher densities.
To highlight the importance of accounting for the effects of age and site, the results of the mean stand development from each thinning scenario in Table 2 were used to plot the relationship between net or gross volume increment and basal area (Figure 4 and Supplementary Figure 2). These figures differ from Figure 3 as both types of volume increment are determined from the simulated height and basal area development as predicted by equations (3) and (4). Thus, Figure 4 and Supplementary Figure 2 are the 2D representations of stand development in NV PAI or GV PAI and basal area, respectively.
A clear unimodal pattern can be seen in both NV PAI and GV PAI . If age were not accounted for this could be mistaken as density effect instead of the natural slowing of growth at higher ages which happens to coincide with stands growing into higher densities. Additionally, Figure 4 shows clear differences in the growth response patterns of thinned treatments among the different site index classes with the curves at lower SI deviating from the curves in the higher SI. Because thinning was implemented at a fixed height in the simulations, the age at which the treatments were applied increased with decreasing SI. Thus, the response patterns seen here also include an age effect when comparing among SI classes. Therefore, it is highly important to separate the age effect with the density effect to determine the direct relationship between volume growth and stand density to avoid confounding the age effect.
Both age and site effects were accounted for in equation (2) by including dominant height and height increment in the model. Additional variables which describe thinning were initially included but were not significant and no obvious bias was seen in the residuals among the thinned and unthinned plots. This suggests that for a given basal area, height, and height increment, any differences in volume increment between thinned and unthinned stands were not large enough to be detected by the modelling approach used in this work.

Relative growth-density relationships
By considering the unthinned stand as the reference stand, relative net (RNVI) and relative gross (RGVI) volume increment were calculated for all nine thinning scenarios at each simulated time step. Plotting those values over stand age shows when the RNVI or RGVI in the thinning treatments surpass that of the unthinned treatment, which has a relative value of one ( Figure 5 and Supplementary Figure 3).
The RNVI in all thinning treatments and for all site indices increased above the unthinned treatment at some point in the rotation ( Figure 5). In general, the age at which the surpassing occurred decreased with increasing site index resulting in the thinning treatments in the higher site index classes having a longer period of increased net growth over the rotation period. Alternatively, there was only one instance in SI-17 and two instances in SI-20 where RGVI was increased above one (Supplementary Figure 3). In all three instances, this was due to a slight increase in basal area of the thinned treatments above that of the unthinned treatment (Supplementary Figure 4).
The relationship between volume increment and stand density in Norway spruce plantations The change in RNVI over the rotation indicates differences between the unthinned and thinned treatments at different time periods which could indicate different response patterns. For clarity, four different time steps were chosen from the simulation data in SI-11 to compare and contrast the change in those relationships over time. Figure 6 shows that initially both NVPAI and RNVI increase with increasing stand density. With time the relationships develop into the constant pattern of Mar:Moller (1947), the optimum pattern of Assmann (1950) and the constantoptimum pattern of Wiedemann (1932). Thus, all three growthdensity patterns were observed in the simulation data and, as hypothesized in Figure 2, those patterns depended on the density of the unthinned control.

Simulation and Monte-Carlo analysis
The previous results indicate that thinning initially decreases both gross and net volume increment due to the immediate decrease in basal area. However, there are cases where both types of volume increment are increased in the thinning scenarios above that of the control. The question becomes how the differences in volume increment among the thinning treatments affect the total cumulative volume over the course of a rotation.
Results from the Monte-Carlo analysis show several cases where cumulative gross volume (CGV) is significantly different between the unthinned and thinned treatments (Figure 7). In site index classes 11, 14 and 17 all once and twice thinned treatments which received the heaviest thinning (50 per cent basal area removal) produced significantly less cumulative CGV than the unthinned treatment, regardless if the thinning was applied early or late. In the site index class 20, only the treatments which received two heavy thinnings had significantly less CGV as compared with the unthinned treatment.
Alternatively, there were no significant differences among all treatment scenarios within a given site index for cumulative net volume production. These results indicate that all eight thinning treatments examined in this work can produce as much net volume over the course of a rotation. Therefore, the recommendation that up to 50 per cent basal area can be removed without having a noticeable effect on total harvested volume is valid for these data under the rotation lengths examined here.

Absolute growth-density relationships
The results from this modelling study which examines gross and net volume increment from thinned and unthinned stands of Norway spruce in Norway suggests different absolute growthdensity patterns between the two types of volume increment. The models presented here indicate that when all other variables are held constant, gross volume increment increases with increasing basal area. However, Figure 3 shows that the pattern is curved and becomes essentially constant at higher densities (basal area > 50 m 2 ha −1 ). Alternatively, in net volume increment the results from the model suggests that volume increment Forestry Figure 4 Two-dimensional representation of net volume increment development for different thinning treatments based on equations (2)-(4) at different site indices (SI, m). The curves represent the integrated effects of both stand age and stand density. Thinning treatments are described with intensity as control (no thinning), medium (25 per cent basal area reduction), or heavy (50 per cent basal area reduction) and with timing as early (first thin at 12 m height) or late (first thin at 16 m height).
increases with increasing basal area up to a point where volume increment is maximized and then decreases with increasing basal area up to the maximum. The decline at higher densities is the result of greater competition and larger mortality as the basal area approaches the maximum. While the pattern is curved, Figure 3 indicates that there is a range in basal area where net volume increment is essentially constant and that this range increases as height increment, and therefore age, decreases. This result indicates that Langsaeter's hypothesized absolute growthdensity pattern is applicable for these Norway spruce data when net volume increment is considered.
In contrast to the work presented here, models developed from Norwegian National Forest Inventory (NNFI) data by Gizachew and Brunner (2011) indicated GV PAI strictly increasing up to a basal area of 60 m 2 ha −1 . Figure 3 shows that based on these data from thinning trials, GV PAI reaches a maximum and begins to flatten out after 50 m 2 ha −1 basal area. One major difference between these two studies is the underlying data source used in the development of the volume growth models. The thinning trials used in this work were established in evenaged Norway spruce plantations, whereas the NNFI data used in Gizachew and Brunner (2011) contain a range of forest types including both planted and naturally regenerated stands, the latter of which can be either even-or uneven-aged. The summary statistics for the two data sources show that the thinning trial data are representative of better stocked stands on higher site productivities as compared with the NNFI data.
Different response patterns in GV PAI from stands of different management types have been seen in other species. Based on thinning studies of loblolly pine in the south-eastern US, Allen II and Burkhart (2018) compared growth-density relationships between non-intensively and intensively managed plantations. In non-intensively managed plantations, growth models produced strictly increasing patterns between GV PAI and basal area and Allen II and Burkhart (2018) noted similarities between their relationship and the relationship of Gizachew and Brunner (2011). Alternatively, growth-density patterns from intensively managed plantations of loblolly pine showed more curvature, flattening out at higher densities, similar to the curves presented in Figure 3. Thus, differences in the response curves can be attributed to differences in site types and management.
The relationship between volume increment and stand density is quite complex and other factors outside of stand density which indicate a point in stand development and productivity are needed. Publications as early as Langsaeter (1941) indicated that the growth-density relationship is a function of stand age and site index. However, these variables are merely indicators of a point in stand development based on an average productivity. Further, there may be cases where stand age and site index are either unknown or roughly estimated which could lead to The relationship between volume increment and stand density in Norway spruce plantations Figure 5 Relative net volume increment (RNVI) development with stand age of thinned stands as compared with an unthinned (RNVI = 1). Thinning treatments are described with intensity as medium (25 per cent basal area reduction) or heavy (50 per cent basal area reduction) and with timing as early (first thin at 12 m height) or late (first thin at 16 m height). misleading results. The models presented in this work use height and height increment which in addition to explaining the effects of stand age also can explain the effects of current site productivity. Such models may be useful for forecasting volume increment from an existing forest inventory or in simulation studies when coupled with a height-age equation which incorporates climate variables. Further, in stands with non-homogeneous site productivities, such as are often found in Norway, even greater precision for forecasting volume increment could be found when connected to height-growth information obtained from remote sensing.

Relative growth-density relationships
As an alternative to absolute growth-density relationships, researchers have often employed relative growth-density relationships to examine differences among thinned and unthinned experimental plots. From the plethora of publications in the literature concerning relative growth-density relationships, it appears that much effort has been given to establish one unifying relationship or pattern. While there are conflicting claims and counter-claims concerning the basic relationships hypothesized by Langsaeter (1941), Mar:Moller (1954 and Assmann (1970), those patterns are often misspecified. Langsaeter's hypothesis is concerned with the absolute growth-density relationship, whereas the patterns of Mar:Moller and Assmann are concerned with relative growth-density relationships.
As hypothesized in Figure 2, Langsaeter's absolute growthdensity relationship could occur simultaneously with the relative growth-density relationships of Wiedemann (1932), Mar:Moller (1954 and Assmann (1970) due to the specification of the unthinned stand as the reference stand. The results presented in this study suggest that this hypothesis is plausible for net volume increment (e.g. Figure 6). While relative gross volume increment increases with increasing RBA, net volume increment of thinned stands is maximized at a basal area of 43 m 2 ha −1 and decreases at higher densities due to mortality. As thinning is generally applied relatively early in the rotation, the unthinned stand is often not at the maximum sustainable density on a given site, especially in plantations where extreme planting densities are not likely to occur. The modelling results show here that as the unthinned and thinned plots move towards the maximum density the relative growth-density relationship shifts from increasing to constant to optimal ( Figure 6). Thus, different patterns can be seen depending on the density of the reference stand and conflicting claims about the relationship may be an artefact of comparing different data sources. Figure 6 Absolute and relative relationships between net volume increment and basal area based on simulated thinning scenarios for SI-11 at different stand ages. Different relative growth-density patterns can also occur depending on the methodology used to examine those relationships. For example, in thinning studies of Norway spruce and Scots pine in Finland and Sweden, Isomäki (2004a, 2004b) and Nilsson et al. (2010) calculate relative volume increment and RBA as averages over a 30-year period. For comparison, the periodic annual volume increment and average basal area over a 30-year period following thinning were calculated from the simulations in this study and the relative relationships between thinned and unthinned scenarios were plotted. Similar to the results from Mäkinen and Isomäki (2004a), Figure 8 shows purely increasing relationships when these data are averaged over this 30-year period. This is not surprising given the results shown in Figure 5 and Supplementary Figure 3 as the unthinned control was still well below the maximum density during this period. While this methodology is valid for examining differences in thinned and unthinned experimental plots over a long growth period, it lacks the ability to detect different relative growth-density patterns which occur at different points in stand development. These results show that it is very important for studies to properly define both the volume increment analysed and the stand density so that appropriate comparisons among studies can be made.

Total volume production
While showing that the growth-density patterns are more dynamic than originally thought the questions becomes how much it matters and what effect does thinning have on total production at the harvest age. For example, in terms of net increment do the periods of increased growth counter the initial reduction in increment immediately following thinning. The results of the Monte-Carlo analysis indicated that there were no differences, according to over-lapping confidence intervals, among cumulative net volume production at the end of a rotation within a given site index between the nine thinning scenarios examined (Figure 7). Alternatively, one or two heavy thinnings reduced CGV production below the unthinned stand, depending on the site productivity. As the difference between gross and net volume is mortality, these results indicate that those treatments can produce as much total net volume while additionally reducing volume lost to mortality. Thus, the recommendation that up to 50 per cent basal area can be removed without having a noticeable effect on total net volume production over a rotation is plausible for the thinning scenarios and rotation lengths examined here.
In a previous analysis of these thinning trial data with measurements covering 25-30 years after thinning, Braastad and Tveite (2001) concluded that cumulative net volume increased with increased density with the largest volumes occurring in the unthinned plots. Figure 8 shows that this is due to a decrease in the average relative net volume increment with increasing thinning strength over the same period. The advantage of thinning on total net volume production is by reducing the time of reduced net volume growth at higher densities. Modelling results suggest that this can result in larger relative net volume growth later in  Table 2. the rotation (e.g. Figure 8). Therefore, the concluded similarities between cumulative net volumes between the thinning scenarios at the end of the rotation are highly dependent on the length of the rotation period. In this study, the rotation length was defined as the age where mean annual increment culminates. However, if the rotation length is shortened then these results may not hold.
Different effects of thinning on total gross and net volume production have been found in other Fennoscandia countries. Nilsson et al. (2010) found that only heavy thinning significantly reduced CGV production but show that light thinning shows an increase in cumulative net volume production. Alternatively, Mäkinen and Isomäki (2004a) found no significant differences in CGV between thinned and unthinned stands. Although, the results of these studies should be compared with caution as those analysis did not or could not account for differences in site productivity and consist of data with different thinning treatments and timing of thinning.

Conclusion
The relationship between volume growth and stand density is dynamic with differences occurring between different measures of volume growth. From the models produced in this work, when all other variables were held constant, gross volume increment increased with increasing basal area up to the maximum basal area on a given site. Alternatively, net volume increment followed an optimal pattern and was maximized around a basal area of 43 m 2 ha −1 . These differences indicate the importance of properly defining the measure of volume increment used in growthdensity studies and the need to ensure comparisons are made among the same measures.
While the relationship between net volume increment and basal area is curved, the model results indicated that there was a range of basal areas where net volume increment was essentially the same, varying less than 1 m 3 ha −1 year −1 and this range varied with current height and height increment. This result would confirm Langsaeter's hypothesis which says that volume growth is constant and optimal across a wide range of densities and that the range in density is dependent on age and site.
As hypothesized in this work, the confirmation of Langsaeter's hypothesis does not negate different response patterns in relative growth-density relationships but instead provides reasoning as to why different patterns may have been observed in past work. Due to the specification of maximum density for a given site being the density of an unthinned stand, all three hypothesized relative growth-density response patterns (i.e. increasing, constant, optimal) can occur depending on the absolute density of the reference stand. Figure 8 Relationship between gross or net volume periodic annual increment over the 30-year period following thinning and the midpoint (average) basal area over the same period. Volume increment and basal area of the thinning treatments have been made relative to the unthinned treatment.

Forestry
Results of simulated thinning scenarios indicated no significant differences in total net volume production at harvesting age between unthinned stands and various thinning strengths up to 50 per cent basal area removal. Alternatively, only two heavy thinnings (50 per cent basal area removal) significantly reduced total gross volume production, as compared with an unthinned stand, at the end of a rotation. This result suggests that for Norway spruce stands two heavy thinnings can produce as much total net volume at the end of a rotation, providing some earlier return on investment, while at the same time avoiding densityrelated mortality. However, these results should be taken with caution as the simulated results have been extrapolated outside of the range of stand ages used in model fitting. It is highly important that these thinning trials continue to be monitored throughout the rotation.
In conclusion, the results of this work rectify alternative historic hypotheses concerning the relationship between volume increment and stand density. The system of equations developed here indicates that different growth-density patterns can occur over a rotation and thus past controversy could have arisen from the desire to develop one unified pattern of response and the comparison of response patterns across different data sources. By use of the modelling framework presented in this work, differences in species reactions to thinning could be evaluated. Further, this modelling framework allows for the analysis of thinning trials more correctly by incorporating the effects of stand age, site productivity, thinning onset and thinning frequency to properly understand the growth-density relationship.

Supplementary data
Supplementary data are available at Forestry online.