A stand-level growth and yield model for thinned and unthinned managed Norway spruce forests in Norway

ABSTRACT

A new stand-level growth and yield model, consisting of component equations for stand volume, basal area, survival, and dominant stand height, was developed from a dataset of long-term trials for managed thinned and unthinned even-aged Norway spruce (Picea abies (L.) Karst.) forests in Norway. The developed models predict considerably faster growth rates than the existing Norwegian models. Further, it was found that the existing Norwegian stand-level models do not match the data from the thinning trails. The significance of thinning response functions indicated that thinning increases basal area growth while reducing competition related mortality. No significant effects of thinning were found in the dominant stand height growth. Model examination by means of cross-validation indicated that the models were unbiased and performed well within the data range. An application of the developed stand-level model highlights the potential use for these models in comparing different management scenarios.

KEYWORDS:

• Volume
• basal area
• dominant height
• survival

Introduction

Since the 1950s, a large portion of forested lands in Norway have been planted with pure stands of Norway spruce (Picea abies (L.) Karst.) and Norway spruce is by far the most economically important tree species responsible for app. 80% of annual cut. At present, as much as 60% of harvested forests in Norway are regenerated by planting Norway spruce (Granhus et al. 2016). More recently the conversion of low-density, low-productivity broadleaved stands to higher producing Norway spruce stands has been considered in Norway as a means to increase carbon sequestration (Søgaard et al. 2019). Given the importance of Norway spruce for the Norwegian forest sector, it is central to have tools that can be used to make informed decisions about Norway spruce management.

During the 1960s to1980s there was a strong emphasis on the development of whole-stand models in Norway (e.g. Tveite 1977; Braastad 1982; Blingsmo 1984). The increased computing power in the 1990s and the availability of the first remeasurement of the permanent sample plots in the Norwegian National Forest Inventory (NNFI) induced a shift of interest towards individual-tree models (e.g. Eid and Tuhus 2001; Andreassen and Tomter 2003; Bollandsås 2007; Bollandsås and Næsset 2009). While individual-tree models provide more detailed information, as compared to whole-stand models, they also require more detailed information to be used in growth and yield simulation. Alternatively, whole stand models require relatively less detailed information for simulation and can lead to smaller stand-level prediction errors (Vanclay 1994; Weiskittel et al. 2011). As the majority of forest inventories are carried out through the use of fixed or variable-radius plots, practical forest management planning in Norway is mainly performed at the stand level, as opposed to the individual-tree level (Eid 2001). Thus, there is still the need for reliable whole-stand models in the Norwegian forest sector.

The currently available stand-level forest planning and scenario tools in Norway are AVVIRK-2000 (Eid and Hobbelstad 2000) and GAYA-JLP (Hoen and Eid 1990; Lappi 1992; Hoen and Gobakken 1997). The core of the growth model used in these planning and scenario tools consist of functions for quadratic mean diameter increment (Blingsmo 1984), height-development (Strand and Braastad 1967; Tveite 1976, 1977) and mortality (Braastad 1982). More recent developments of stand-level models include functions for stand-level volume (Næsset and Tveite 1999), mortality (Eid and Øyen 2003), height-development (Sharma et al. 2011), and prediction equations for quadratic mean diameter and trees per hectare (Eid 2001). Empirical models are generally only appropriate if they are applied within the range of data that was used for model development. The existing stand-level models in Norway are relatively old and their fitting data are even older. Since, in general, empirical models are only appropriate within the range of data that was used for model development, it is questionable if they represent the current growing conditions (climate) and management practices as needed.

Historically, low pulpwood values along with high costs of thinning operations have resulted in relatively low thinning intensity (no thinning or one late thinning) in spruce forests in Norway. However, recent increases in pulpwood prices and an improved demand for timber have resulted in increased interest in thinning of Norway spruce. The existing Norwegian stand-level models for even-aged spruce stands do not explicitly model differences in growth following thinning (no thinning response). It is widely accepted that there is a thinning response in Norway spruce and results from thinning trials of Norway spruce in Sweden and Finland have shown that thinning influences basal area growth and survival (Mäkinen and Isomäki 2004; Nilsson et al. 2010). Thus, not explicitly representing the thinning response in basal area growth and survival can lead to biased estimates of future density and volume. The lack of thinning effects in the current stand-level growth and yield models used in Norway suggests a gap in the tools necessary for determining optimal stand-level management strategies. Improved models that include a thinning response will provide forest managers with a better tool for planning stand level management of Norway spruce.

The objective of this study was to develop a stand-level model for forecasting growth and yield of thinned and unthinned Norway spruce stands. As such, data from long term thinning trials were used to develop a system of stand-level equations for predicting stand volume, basal area, survival, and dominant height. Comparisons of growth and yield projections among the new models developed here and the most recent stand-level models were made to show the improvements of the functions. Further, the model behavior is illustrated and used to examine the interactions of the models and effects of thinning on stand-development.

Materials and methods

Data

The data used in this study consists of a series of silvicultural trials with 94 unthinned and 135 thinned permanent sample plots, ranging from 0.4–0.6 hectares in size, established in even-aged planted Norway spruce stands. Establishment of the trials occurred between 1969 and 1976 in stands aged 20–31 years. Plot measurements were scheduled in five-year intervals, with a few deviations from the schedule, and to date each plot has received between 5 and 10 measurements. Within these data are 311, 442, and 246 growth observations in unthinned, once thinned stands, and twice thinned stands, respectively (Table 1). The trials are distributed across the natural range of Norway spruce in Norway, ranging from 58 to 67 degrees in latitude. The full range of the sample plot data are presented in Figure 1. A summary of these trials has been previously provided by Braastad and Tveite (2001).

Figure 1. Distribution of Norway spruce sample plots across Norway.

Table 1. Descriptive statistics of thinned and unthinned Norway spruce data.

Treatment	Growth periods		A	S	N	G	V	GTQ	HT	ALT	LAT
Unthinned	311	mean	41	16	2510	37.3	298			209	61
		SD	13	3	860	15.3	196			140	2.7
		min	18	11	777	6.1	21			20	58
		max	92	21	5021	79.8	893			530	67
1 Thinning	255	mean	50	16	1221	36.1	332	0.71	12.4	212	60
		SD	12	3	400	12.4	198	0.11	2.3	146	2.5
		min	25	11	600	6.3	26	0.47	8.5	20	58
		max	68	21	2640	67.4	917	0.92	21.7	450	64
2 Thinnings	187	mean	38	15	1796	25.4	149	0.79	11.4	260	61
1st thin		SD	8	2	312	9.2	82	0.08	1.8	132	2.3
		min	25	11	900	8	27	0.5	8.7	40	58
		max	61	20	2500	46.8	410	0.91	17.2	450	64
2nd thin	246	mean	55	15	1043	35.7	331	0.75	15.6
		SD	11	2	232	10.1	171	0.12	2
		min	31	11	460	11.8	72	0.35	11.8
		max	78	20	1700	60.3	812	0.92	21.1

Notes: A = stand age from planting (years); S = site index (m) at base age 40 years; N = number of trees (ha⁻¹); G = basal area (m² ha⁻¹); V = stand volume (m³ ha⁻¹); G_TQ = basal area thinning quotient (G after thinning divided by G before thinning); and HT = dominant stand height (m) at thinning, ALT = altitude (m); LAT = latitude.

The trails were established to examine the effects of thinning from below with focus on thinning intensity and timing of thinning on volume production and stand structure. Initially, the trials were planned to include treatments of different density reductions at different stages of stand development. Treatments included one to three thinnings from below with different combinations of target residual trees per hectare (TPH) of 2070, 1600, 1100, and 800 at stand dominant heights of 8, 12, 16, and 20 meters. Over the course of the experiment, the target residual TPH were often missed (thinning too light or too heavy) and/or the timing of thinning was performed too early or late. Further the experimental design lacked replication. The deviation from the experimental design restricts the examination of treatment groups but fortunately 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.

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). In-growth trees were included in plot measurements once reaching a dbh of 5 cm.

Model development

The whole-stand growth and yield model is comprised of individual equations for stand volume, basal area, survival, and dominant stand height. During the model development stage, multiple models for each stand component from the literature were compared to determine the best final model. In all cases, the Bayesian Information Criterion (BIC) was used in order to compare non-linear models with varying numbers of parameters (Claeskens and Hjort 2008). A list of variable names and definitions used in the following modeling framework is provided in Table 2. Further, a flow chart of the model application is provided in Figure 2.

Figure 2. Description of the model application in the new stand-level model. Numbers in brackets corresponds to the equations presented in the text.

Table 2. Definitions of variables and their abbreviations used in this work.

Abbreviation	Definition of Variable
A1,2	Stand Age (years from planting)
V1,2	Volume (m^3 ha^−1)
G1,2	Basal area (m^2 ha^−1)
H1,2	Dominant stand height (m), mean height of 100 thickest trees per hectare
N1,2	Number of trees per hectare
GA	Basal area after thinning (m^2 ha^−1)
GB	Basal area before thinning (m^2 ha^−1)
HT	Dominant stand height at thinning (m)
S	Site Index at base age 40 years (m)
NA	Trees per hectare after thinning
NB	Trees per hectare before thinning
QMD1,2	Quadratic mean diameter (cm)
CAI	Current annual volume increment (m^3 ha^−1 yr^−1)
MAI	Mean annual volume increment (m^3 ha^−1 yr^−1)

Note: All subscripts of 1 or 2 denote stand values at the measurement period at times 1 or 2, respectively.

Volume: For predicting stand-level volume, a non-linear version of the Schumacher (1939) volume equation was developed:(1) $V_2={\text{β}}_1{G}_2^{{\text{β}}_2}{H}_2^{{\text{β}}_3}exp\left(\frac{{\text{β}}_4}{A_2}\right)$(1) where ${\text{β}}_i$ are parameters to be estimated. Additional variables which described thinning response were initially included but none were significant in the final model and no patterns in the residuals were detected between thinned an unthinned stands. Predictions of future basal area (G₂) and future height (H₂) are needed for predicting volume with Equation (1).

Basal Area: Projection of stand basal area was modeled using a modified version of the equation of Hasenauer et al. (1997):(2) $G_2=G_1^{H_1/H_2}exp\left({\text{β}}_1{\left(\frac{N_2}{N_1}\right)}^{{\text{β}}_2}\left(1-\frac{H_1}{H_2}\right)TR\right)$(2)

The thinning response function (TR) is defined as:

$TR={\left({\displaystyle \frac{G_A}{G_B}}\right)}^{{\text{β}}_3{\displaystyle \frac{H_T}{H_2}}}$

This thinning response function allows for an increase in basal area growth following thinning for a period, but with time the basal area growth returns to that of the unthinned stand (Hasenauer et al. 1997). In addition to predictions of future dominant stand height (H₂), predictions of future trees per hectare (N₂) is also needed for predicting basal area with Equation (2).

Survival: Multiple stand-level survival models based on Zhao et al. (2007), with the inclusion of a thinning effect, were tested. The final best fitting model, similar to that of Gyawali and Burkhart (2014), was:(3) $N_2={\left[N_1^{{\text{β}}_1}+{\text{β}}_2\left(\frac{G_A}{G_B}\right){\left(\frac{S}{1000}\right)}^{{\text{β}}_3}\left(A_2^{{\text{β}}_4}-A_1^{{\text{β}}_4}\right)\right]}^{1/{\text{β}}_5}$(3)

Dominant Height: Stand height projection was modeled using a GADA formulation of the Chapman-Richards type height-age model (Cieszewski and Bailey 2000). Multiple model forms from Diéguez-Aranda et al. (2005) were tested, with the following model providing the best fit:(4) $H_2=H_1{\left(\frac{1-exp\left(-{\text{β}}_1{A}_2\right)}{1-exp\left(-{\text{β}}_2{A}_1\right)}\right)}^{{\displaystyle \frac{{\text{β}}_2+{\text{β}}_3}{X_0}}}$(4) where$X_0=0.5\left(log\left(H_1\right)+{\text{β}}_{2}L+\sqrt{{\left(log\left(H_1\right)+{\text{β}}_{2}L\right)}^2-4{\text{β}}_{3}L}\right)$ $L=log\left(1-exp\left(-{\text{β}}_1{A}_1\right)\right)$

Tree reduction: Thinning treatments are often specified as a reduction in either tree number or basal area. In the case that one of these variables is known following thinning and the other is not, for example under simulation studies, it will be necessary to predict the reduction in the other. To predict the number of trees remaining after a thinning when the reduction in basal area after thinning is known, the following model was developed:(5a) $\frac{N_A}{N_B}=exp\left({\text{β}}_1+{\text{β}}_2\left(\frac{G_A}{G_B}\right)\right)$(5a)

In the case that the number of trees following thinning is known and the basal area following thinning is unknown then Equation (5) can be solved to determine the reduction in basal area as:(5b) $\frac{G_A}{G_B}=\frac{log\left({\displaystyle \frac{N_A}{N_B}}\right)-{\text{β}}_1}{{\text{β}}_2}$(5b)

Model fitting and validation

Equations (1–5) were fit using non-linear seemingly unrelated regression using the “systemfit” package (Henningsen and Hamann 2007) in the open source statistical software R version 3.6.1 (R Core Team 2019). Seemingly unrelated regression accounts for the contemporaneous correlations among different regression equations and has been widely used in the development of growth and yield models (Furnival and Wilson 1971; Borders 1989; Robinson 2004; Coble 2009; Zhao et al. 2015). The system of equations was fit to the trial data (Table 1) as a growth series over each measurement period (5-year intervals). The growth series was constructed by considering all measurement intervals in succession. For example, for a plot with one initial measurement and two additional re-measurements, two growth periods were constructed.

The goodness of fit of all component equations of the whole-stand model was evaluated using four fit statistics: the mean residual (E), the mean of the absolute values of the residuals (MAE), root mean square error (RMSE), and the coefficient of determination (R²). These four criteria can be expressed mathematically as:(6) $E=\frac{{\sum }_{i=1}^n\left(Y_i-{\hat{Y}}_i\right)}{n}$(6) (7) $\mathrm{M}\mathrm{A}\mathrm{E}=\frac{{\sum }_{i=1}^n\left|Y_i-{\hat{Y}}_i\right|}{n}$(7) (8) $\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}=\sqrt{\frac{{\sum }_{i=1}^n{\left(Y_i-{\hat{Y}}_i\right)}^2}{n}}$(8) (9) $R^2=1-\frac{{\sum }_{i=1}^n{\left(Y_i-{\hat{Y}}_i\right)}^2}{{\sum }_{i=1}^n{\left(Y_i-{\overline{Y}}_i\right)}^2}$(9) where $Y_i$, ${\hat{Y}}_i$, and ${\overline{Y}}_i$ represent the observed, predicted, and mean values of each component, respectively, and n represents the number of observations.

Model validation was performed using two-fold cross-validation. The trial data were randomly split into two datasets roughly equal in number of observations, one validation and one fitting dataset. Equations (1–5) were fitted to the fitting dataset and then the validation dataset was used to calculated summary statistics (Equations 6–9) for the validation. A similar validation approach using Leave-one-out (LOO) cross-validation was also performed but provided similar fit statistics as the two-fold cross-validation. Since the LOO validation did not provide additional information those results are not shown.

Table 3. Stand data for four selected unthinned Norway spruce plots.

ID	S	A	N	G	H	QMD	V
1	11	46	2977	31.3	12.2	11.6	155.9
2	14	34	3090	26.2	12.2	10.4	125.9
3	17	25	2038	22.7	11.8	11.9	109.6
4	20	22	2710	24.6	10.7	10.8	115.9

Note: A = stand age from planting (years); S = site index (m) at base age 40 years; N = number of trees (ha⁻¹); G = basal area (m² ha⁻¹); H = stand dominant height (m), QMD = quadratic mean diameter (cm); and V = stand volume (m³ ha⁻¹).

Model comparison

To show the improvements of the new model, Equations (1–4), comparisons were made against the most recent existing functions for each of the respective components. Stand volume and mortality were compared with the functions of Næsset and Tveite (1999) and Eid and Øyen (2003), respectively. As there is no existing stand-level basal area growth equation for Norway spruce in Norway, the implied future quadratic mean diameter (QMD₂, cm) as calculated from predictions of G2 and N2 using Equations (2) and (3) where compared to the predictions of QMD₂ using the function of Blingsmo (1984). Each of the component equations were evaluated based on residual analysis and fit statistics from Equations (6–9).

Lastly, height-age relationships from Equation (4) were compared to those predicted from the site index function of Tveite (1977) and the height-age model of Sharma et al. (2011). As Equation (4) was based on age from planting, the site index curves produced from the functions of Tveite (1977) and Sharma et al. (2011) were adjusted to reflect the age from planting, instead of the age from which most trees reach 1.3 m in height, using age correction factors published in Søgaard et al. (2019).

Stand development in thinned and unthinned stands as implied by the models was examined by simulating forest growth under different stand productivities and three different thinning scenarios. One measurement from each of four plots from the fitting dataset, representing site indices of 11, 14, 17, and 20 m, where selected to demonstrate each component model (Table 3). Growth for each stand in Table 3 was simulated in 5-year time steps according to the work flow presented in Figure 2. The thinning scenarios examined were based on representative thinning treatments performed in the thinning trial data and included an unthinned scenario, a one-thinning scenario where 25% basal area is removed at a stand dominant height of 12 m, and a two-thinning scenario where 25% basal area is removed once at a stand dominant height of 12 m and then a second time at a stand dominant height of 16 m. All figures were generated using the ggplot2 package in R (Wickham 2009).

Table 4. Parameter estimates (standard errors) and fit statistics for Equations (1–5).

	Equation – component
	1 – V2	2 – G2	3 – N2	4 – H2	5 – NA/NB
${\text{β}}_1$	0.24961	4.77696	−1.0085	0.01605	−1.93267
	(0.01097)	(0.02191)	(0.2619)	(0.00142)	(0.02303)
${\text{β}}_2$	1.15036	0.30957	0.03675	0.61208	1.92953
	(0.01594)	(0.02651)	(0.00792)	(0.17427)	(0.02439)
${\text{β}}_3$	1.01153	−0.1479	3.76228	4.43722
	(0.01224)	(0.0175)	(0.7144)	(0.54254)
${\text{β}}_4$	2.320398		2.5541
	(0.59361)		(0.4084)
${\text{β}}_5$			−1.0097
			(0.2631)
	All Data
E	0.325	−0.012	−1.78	0.019	−0.00044
MAE	11.8	1.33	73.3	0.499	0.0436
RMSE	16.6	1.8	138.1	0.644	0.0533
R^2	0.987	0.963	0.965	0.976	0.846
	Cross Validation
E	0.312	−0.0945	−1.23	−0.026	−0.00043
MAE	11.6	1.32	73.5	0.472	0.0441
RMSE	16	1.83	158	0.627	0.0542
R^2	0.988	0.962	0.949	0.976	0.845

Results

Parameter estimates and fit statistics for Equations (1–5) are presented in Table 4. Estimated parameters for all equations were significant at the α = 0.01 level. When comparing the observed vs predicted values for each component in Equations (1–4) no obvious bias was seen (Figure 3) and further residual analysis (results not shown) indicated no obvious bias with stand age, site quality, degree of thinning, altitude, or latitude. Further, when estimated values of N₂ and H₂ were used in the prediction of G₂ and then estimates of H₂ and G₂ in the prediction of V₂, no bias was seen in the model predictions (Figure 3). Fit statistics for the total dataset and for the two-fold cross-validation indicated that the models fit well and explained a large portion of the variation, with R² values greater than 0.96 for Equations (1–4) and greater than 0.84 for Equation (5) (Table 4).

Figure 3. Observed versus predicted values for future Volume (V₂), basal area (G₂), trees per hectare (N₂), and dominant stand height (H₂) based on predictions from Equations 1–4 with parameter values from Table 2. Residuals are based on direct predictions from observations (circles) and, in the case of G₂ and V₂, include predictions from the necessary model components (triangles).

Dominant height projection

In the model development stage, the effect of thinning on height development was examined using indicator variables and thinning response variables. In no case were any parameters which including a thinning effect significant. Further, there was no obvious bias in dominant height projection when no thinning effects were included in the model. These results suggest that thinning from below does not significantly influence dominant height development in even-aged Norway spruce stands in Norway.

A comparison of height-age curves produced from Equation (4) and those produce by Tveite (1977) and Sharma et al. (2011) indicates differences in height development at older ages (Figure 4). The curves of Tveite (1977) tended towards lower heights at later ages, followed by the curves of Sharma et al. (2011), as compared to the curves from Equation (4).

Figure 4. Comparison of site-index curves between the new Equation (4) and the equations of Tveite (1977) and Sharma et al. (2011).

Survival

Parameters for site index and thinning response were significant in Equation (3) indicating that both stand productivity and degree of thinning effect future trees per hectare (N₂) (Table 4). The effect of these two variables is shown in Figure 5 which depicts the development of N₂ from a common starting value of 2500 trees per hectare under site indices 11, 14, 17, and 20 m and for different density removals as indicated by the ratio of basal area before and after thinning (G_TQ). The models indicate that the mortality rate is greater in stands of higher productivity. Further, the mortality rate is generally less in thinned stands.

Figure 5. Development in number of trees per hectare for different site indices (m, base age 40 years) and different thinning intensities based on a common starting value of 2500 trees per hectare. The basal area thinning quotient (G_TQ) is the ratio of basal area before and after thinning where a value of 1 indicates the unthinned stand condition.

Comparisons of the residuals between Equation (3) and those produced from the function of Eid and Øyen (2003) indicates the improvement of the new model (Figure 6). For both thinned and unthinned treatments, the function of Eid and Øyen (2003) over-predicted future N₂ which indicates an underprediction in mortality. While, the variance explained by the model of Eid and Øyen (2003) was high (R² = 0.94), the mean error was −86.8 trees per hectare, indicating the obvious bias seen in Figure 6.

Figure 6. Comparison of residuals for predicted future trees per hectare (N₂) between the new survival function, Equation (3), and the survival function of Eid and Øyen (2003).

Basal area and quadratic mean diameter

The significant parameter in the thinning response function (TR) in Equation (3) indicates that thinning positively influence basal area growth. Additional parameters to distinguish between one or two thinnings were also examined but none proved significant in the model. Further, no obvious bias was observed in the residuals between stands thinned once or twice (e.g. Figure 3).

As there is no stand-level basal area growth model for thinned or unthinned stands of Norway spruce stands in Norway, comparisons were made against the current diameter growth model of Blingsmo (1984). For this analysis the future quadratic mean diameter (QMD₂) was calculated from predictions of future basal area (G₂) and N₂ using Equations (2) and (3). Comparisons of the residuals for QMD₂ as calculated from the new models indicated an improvement over the function of Blingsmo (1984) (Figure 7). The diameter increment function of Blingsmo (1984) tended to under-predict QMD₂ whereas the implied QMD₂ as calculated from Equations (2) and (3) was unbiased.

Figure 7. Comparison of residuals for future quadratic mean diameter (QMD2, cm) between the function of Blingsmo (1984) and the implied QMD2 calculated from future basal area and trees per hectare as predicted from Equations (2) and (3), respectively.

Stand volume

Equation (1) predicts stand-volume as a function of basal area, stand dominant height, and stand age. Inclusions of parameters which explain differences between thinned and unthinned stands were not significant and residuals indicated no obvious bias in thinned and unthinned stands. Thus, Equation (1) can be equally applied in both stand types. Comparisons of residuals indicated that Equation (1) is unbiased, whereas, the volume function of Næsset and Tveite (1999) tended to over-predict volume when stand volumes were lower than 400 m³ ha⁻¹ and then under-predict volume thereafter (Figure 8).

Figure 8. Comparison of residuals for future volume (V₂, m³ ha⁻¹) as predicted from Næsset and Tveite (1999) and from Equation (1).

Stand development

Based on the simulated stand growth from initial stand states (Table 2), the results indicate the different responses from site productivity and thinning treatment. As seen in Figure 4, mortality is much higher on better-quality sites, gradually decreasing with decreasing site quality (Figure 9). Such results are expected as better-quality sites will grow into the self-thinning range much faster than lower-quality sites. As there were no significant effects of thinning in Equation (4), dominant height development was the same for all three scenarios.

Figure 9. Development in future trees per hectare and basal area for four Norway spruce stands (Table 3) representing different site indices (S-11, 14, 17, and 20 m) with three management scenarios including: Unthinned, one-thinning removing 25% of the basal area at a dominant stand height of 12 m, and two-thinnings removing 25% of the basal area at dominant stand heights of 12 and 16 m.

The effects of thinning on basal area development varied with site quality and thinning treatment. In the lowest site indices (S-11), one or two light thinnings (25%) did not result in a noticeable increase in basal area growth and by 120 years of age basal areas in the thinned scenarios remained well below the unthinned scenario. However, with increasing site quality the effect of thinning on basal area growth increased. The modeling results further indicate that thinned stands on good sites can potentially achieve a higher maximum basal area as compared to unthinned stands on good sites.

While maximum size-density relationships (MSDR) where not explicitly modeled for these data, the implied MSDRs can be examined by calculating the quadratic mean diameter (QMD, cm) from basal area and trees per hectare (N) and then relating N to QMD. In order to determine the MSDR line, the equation of Reineke (1933) was fit to using the fitting data described in Table 1 and was determined to be:$\log \left(N\right)=12.35-1.58\log \left(\mathrm{Q}\mathrm{M}\mathrm{D}\right)$where all variables are described as previously. Note the objective was not to determine the best fitting method for MSDRs but to produce a general MSDR line from which the implicit MSDRs in the stand models could be examined. Figure 10 indicates that the implied MSDRs from the basal area and survival models do coincide with the MSDRs observed in the fitting data.

Figure 10. Maximum size-density relationships for four Norway spruce stands (Table 3) representing different site indices (S-11, 14, 17, and 20 m) with three management scenarios including: Unthinned, one-thinning removing 25% of the basal area at a dominant stand height of 12 m, and two-thinnings removing 25% of the basal area at dominant stand heights of 12 and 16 m.

The lack of thinning response in basal area growth in site indices S-11 and S-14 resulted in less total standing volumes in the thinned scenarios as compared to the unthinned scenarios by 120 years of age (Figure 11). However, in S-17 the unthinned and once thinned scenarios produced similar standing volumes by age 120 and in S-20 the once-thinned stand produced more standing volume. When comparing the total potential volume removed, net volume (standing plus removed volume) was similar among all three thinning scenarios in S-11 and S-14. Net volume in site indices S-17 and S-20 among the thinned scenarios surpassed the unthinned scenario at later ages.

Figure 11. Development of standing volume and net volume (standing + removed volume) for four Norway spruce stands (Table 3) representing different site indices (S-11, 14, 17, and 20 m) with three management scenarios including: Unthinned, one-thinning removing 25% of the basal area at a dominant stand height of 12 m, and two-thinnings removing 25% of the basal area at dominant stand heights of 12 and 16 m.

Thinning also effects the current annual net volume increment (CAI, m³ ha⁻¹ yr⁻¹) with the effect dependent on site-productivity, the thinning treatment, and the amount of time since thinning. Figure 12 shows for these sample data (Table 2) that, among the thinning scenarios examined here, thinning initially reduced the CAI in all site indices. In the lowest site indices (S-11) the CAI in the one-thinning scenario was able to converge to the CAI of the unthinned scenario after a short time, but CAI in the two-thinning scenario remained below that of the unthinned stand. With increasing site productivity, the effect of thinning had an increasingly favorable response on volume growth and the models indicate that in S-14, −17, and −20 thinning can result in increases CAI in thinned stands above that of the unthinned stand. Such a result indicates the higher net volume as seen in Figure 11.

Figure 12. Comparison of current annual volume increment (CAI) and mean annual volume increment (MAI) stands of different site indices under three management scenarios: Unthinned, one-thinning removing 25% of the basal area at a dominant stand height of 12 m, and two-thinnings removing 25% of the basal area at dominant stand heights of 12 and 16 m.

The rotation age as defined by the culmination of mean annual net volume increment (MAI) also varied depending on site productivity and degree of thinning. For site indices S-11, −14, and −17, MAI similarly culminated in thinned and unthinned stands at 5.5, 8.1, and 12.2 m³ ha⁻¹ yr⁻¹ for those site indices, respectfully. However, in S-20 MAI culminated at 12.7, 14.7 and 13.7 m³ ha⁻¹ yr⁻¹ for the unthinned, once-thinned, and twice thinned scenarios, respectfully. Harvesting ages according to maximum MAI were similar among thinned and unthinned stands in S-11 and S-14 but thinning increased the age of maximum MAI in S-17 and S-20 by 4–5 years.

Discussion

Forest management planning in Norway is generally conducted based on information at the stand-level (Eid 2001), as opposed to the individual-tree level, and reliable stand-level models are needed for forecasting current stand conditions under different management scenarios to make informed decisions about treatment schedules. To fill this need, new stand-level models were developed based on an extensive dataset of thinned and unthinned long-term permanent sample plots in managed even-aged Norway spruce stands.

It is noteworthy that the new models predict considerably faster growth and mortality rates than the equations of Blingsmo (1984) and Eid and Øyen (2003), respectively. The growth function of Blingsmo (1984) tended to under predict future quadratic mean diameter while the survival function of Eid and Øyen (2003) tended to over predict future trees per hectare (e.g. Figures 6 and 7). Increased mortality rates are linked to increased growth rates as forest stands will grow into the range of distant-dependent mortality (self-thinning) faster. These changes in growth and survival can be primarily linked to differences in the underlying data sources used in model fitting. The data used in model fitting in this work (Table 1) consist of managed Norway spruce stands comprised of newer genetics, whereas the previous models are developed primarily from un-managed, naturally regenerated stands.

The new dominant height projection equation (Equation 4) differs from previously developed models (e.g. Tveite 1977; Sharma and Brunner 2016) for Norway spruce in Norway as the stand age used in Equation (4) is the age from planting instead of the age at which the height of most trees in the stand reaches 1.3 m. The use of age from planting may be more practical for forest owners and managers as the age when most trees reach a height of 1.3 m is subjective and without consistent inventories must me estimated. Thus, the new dominant height model will be useful for projecting height development in Norway spruce plantations where the exact stand age from planting is known.

When compared against the existing height development equations of Tveite (1977) and Sharma et al. (2011) Equation (4) indicated greater height growth at later stand ages across a range of site indices, with larger differences occurring as site index increased (Figure 4). One difference between these three models is the underlying data used in model fitting. Tveite (1977) combined data from naturally regenerated forests, which covered ages of a general rotation in Norway spruce, and young planted stands where the average age was about 35 years. Thus, most of those data represent stands of different genetics, climate, and management. While Sharma et al. (2011) used the same trial data used in this work (Table 1), with fewer measurements, they also incorporated data from the Norwegian National Forest Inventory where older aged stands are predominantly on sites of lower productivity (as indicated in their Figure 2). The inclusion of such data in the development of a height-age model, in addition to representing different genetics and management, could result in high-leverage on the asymptote resulting in the decreased height growth at later ages as seen in Figure 4.

An application of the model is presented to illustrate the model behavior with four selected plots (Figures 11 and 12). This limited exercise is not an exhaustive study of the model behavior, but an example of how the new models can be used to evaluate different management options in even-aged Norway spruce stands. The comparisons of predicted future yields for the four selected plots indicated the potential for greater total harvested volume at later ages when stands are thinned. However, standing volumes of thinned stands in lower site indices did not surpass those of the unthinned stands, indicating lower final yields (Figure 11). Further, modeling results indicate that thinning can potentially increase the rotation age, as defined by culmination of maximum MAI, depending on site quality and thinning treatment (Figure 12).

Results from the new-models indicate that they are unbiased and provide better predictions for growth and yield of thinned an unthinned even-aged Norway spruce stands than the most current stand-level equations available (e.g. Figures 6–8). Additionally, a comparison of predictions of G₂ with estimated values from Equations (3) and (4) and then predictions of V₂ with estimated values from Equations (2) and (4) with observed values of G₂ and V₂, respectively, indicated that the new stand-level model was unbiased when used as a system (Figure 3). The new equations indicate the need to explicitly account for degree of thinning in the stand-level model to account for changes in growth conditions following thinning. Overall the newly developed models shown quite different behaviors from the existing models and may result in deferent decisions about thinning scheduling and the optimal timing of final harvesting.

Equations (1–4) indicate that the primary differences in growth conditions between thinned and unthinned stands occur in basal area growth and survival. However, thinning does not influence dominant height growth. Further, an examination of the model behavior indicates that over a traditional rotation period, total net volume production is similar among thinned and unthinned scenarios (e.g. Figure 11). These results mirror results from analysis of thinning trials of Norway spruce in Sweden and Finland. For example, Mäkinen and Isomäki (2004) concluded that thinning from below in Norway spruce stands in Finland significantly reduced natural mortality and increased diameter growth, while having no significant effect on dominant height growth. The authors show that over a 30-year period there was no significant difference in total net volume production between thinned and unthinned stands. In Sweden, Nilsson et al. (2010) showed similar results in thinning trials of Norway spruce as those of Mäkinen and Isomäki (2004). While the data used in this work is not sufficient for a similar analysis of thinning experiments as performed by Mäkinen and Isomäki (2004) and Nilsson et al. (2010), the modeling results indicate similar reactions to thinning from below in Norway spruce stands in Norway as compared to both Sweden and Finland.

The modeling approach used in this work to incorporate thinning response into the basal area equation differs from the current documented basal area equations used in the decision support systems in Sweden (Heureka) and Finland (MOTTI) which rely on categorical variables to adjust growth over different periods following thinning. For example, the basal area growth model presented in Elfving (2010) incorporates two categorical variables which provide an increase in basal area growth over the period of 0–10 years following thinning and then a smaller increase over the period of 11–25 years following thinning, as compared to a stand with no thinning. A similar approach was used in the basal area growth model of Hynynen et al. (2002), where a categorical variable provides an increase in growth over the period of 0–5 years following thinning as compared to the unthinned stand.

Equation (2) provides a more detailed response in basal area growth to thinning by the inclusion of a continuous thinning response function, as opposed to a discrete categorical variable. This thinning response function accounts for both the intensity of thinning, through the ratio of basal area before and after thinning, and relative time since thinning, expressed as the ratio of stand dominant height at the time of thinning to the future stand dominant height. By expressing time since thinning in this manner, site effects are also accounted for as the ratio between stand dominant height in the future and at thinning for a given period after thinning depends on site productivity.

The data used for the development of the new stand-level model presented in this work come from long-term thinning trials in managed even-aged Norway spruce stands aged 20–70 years. These trials are representative of well stocked stands with adequate control of deciduous competitors. Further, the thinning treatments were implemented mostly from below, removing trees from the lower portion of the diameter distribution. Thus, the models may not properly represent stand development in un-managed forests or in different thinning applications. While the models appear to extrapolate well, following the self-thinning line, some care should be used if applied to stands outside of the range of data used for model fitting.

In conclusion, the new equations for volume, basal area, survival, and dominant stand height developed in this work collectively provide a stand-level model from which growth and yield of managed thinned and unthinned even-aged Norway spruce stands can be forecasted. Further, the stand-level model can be used to compare potential gains or losses in total or net harvested volume from different thinning treatments. The model will prove useful to forest landowners, managers, or analysts for making management decisions in managed even-aged Norway spruce plantations. Further, when coupled with additional auxiliary models, the new stand-level model can provide the core for a decision support system in managed Norway spruce forests.

Acknowledgements

This work was partially supported by the Research Council of Norway through project 301745 (Maximizing the Resilience and Carbon Sequestration in Managed Norway Spruce Forests) and the Forest Growth and Yield Research Group at the Norwegian Institute of Bioeconomy Research. The authors are indebted to the many people who have been involved in collecting and maintain the data series used in this work.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

This work was partially supported by the Research Council of Norway through project 301745 (Maximizing the Resilience and Carbon Sequestration in Managed Norway Spruce Forests) and the Forest Growth and Yield Research Group at the Norwegian Institute of Bioeconomy Research.