Radial trends in black spruce wood density can show an age- and growth-related decline

Abstract

Context

Wood density variation affects structural timber performance and is correlated with several potentially confounding factors, such as cambial age, position in the stem and growth rate. To date, these relationships have not been comprehensively quantified in black spruce (Picea mariana (Mill.) B.S.P.).

Aims

The aim of this study was to describe the variation in annual ring density in black spruce as a function of cambial age, stem height and growth rate.

Methods

Radial density profiles from 107 black spruce trees were analysed using a two-stage modelling approach. First, the parameters of a nonlinear function were estimated separately for individual samples. Linear regression was then used to model the parameters obtained in the first stage as functions of internal and external tree descriptors.

Results

Annual ring density was high near the pith and declined rapidly in the first 15 annual rings before increasing to more stable values between rings 25 and 60. However, just below 25 % of the samples showed a gradual decline towards the bark, typically after ring 60.

Conclusion

Describing and quantifying radial density patterns, including the decline close to the bark, will help further our understanding of the links between tree growth and ring density over the life of the tree.

1 Introduction

Black spruce (Picea mariana (Mill.) B.S.P.) is one of the most important and valuable species for both lumber and pulpwood production in Eastern and Central Canada and an important commercial and reforestation species from the Atlantic Coast to Manitoba (Boyle et al. 1989; Zhang 1998). In the boreal forest of Canada, black spruce is harvested predominantly from natural stands. Consequently, an understanding of the influence of natural forest dynamics on black spruce wood attributes is important in order to formulate appropriate management strategies for maximizing the value of the resource at harvest.

As demand for wood and wood products has increased in recent decades, research efforts have focused more intensively on wood and fibre quality attributes in addition to volume production (Jozsa and Middleton 1994; Evans 1994; Koubaa et al. 2000; Franceschini et al. 2012). Of these attributes, wood density is one of the most important and useful wood quality indicators (Panshin and de Zeeuw 1980; Zhang et al. 1996) and is closely correlated with the mechanical properties of wood (Zhang et al. 1993; Koga and Zhang 2004) and with certain strength properties of paper (van Buijtenen 1982). Furthermore, with growing interest in forest carbon sequestration, wood density can also be used in allometric equations as an essential component of aboveground biomass estimations (Ketterings et al. 2001; De Vries et al. 2006).

Previous studies have quantified the effects of the different sources of variation, i.e. internal variables such as cambial age, ring width and genetic origin (Alteyrac et al. 2005; Jyske et al. 2008; Savva et al. 2010) and external variables such as position along the stem, site quality and climate (Wang et al. 2002; Bouriaud et al. 2004; Guilley et al. 2004), on annual ring density. Within-stem variations in ring density may consequently be related to various intrinsic and extrinsic sources, which have potentially confounding effects. In this context, a modelling approach based on the observed effects of important variables, such as cambial age, annual growth rate and position along the stem, on annual ring density can help further our understanding of these complex interrelationships.

In spruces, the radial pattern of variation in wood density typically shows a rapid decrease near the pith followed by a more gradual increase towards the bark (Panshin and de Zeeuw 1980), a trend that has been confirmed specifically for black spruce (Zhang et al. 2002; Alteyrac et al. 2005). Despite having been investigated extensively, the reported correlation between growth rate and wood density is not consistent among studies. For example, a negative correlation between annual ring width and wood density has been reported in both Norway spruce (Picea abies (L.) Karst.) and Sitka spruce (Picea sitchensis (Bong.) Carr.; Franceschini et al. 2010, 2013; Gardiner et al. 2011). Also in Sitka spruce, Simpson and Denne (1997) reported that the strength of this negative correlation increased from the stem base to the tree tip, for a given cambial age. Conversely, in balsam fir (Abies balsamea), Koga and Zhang (2004) found a significant negative correlation between growth rate and density below 3 m in the stem, but there was no significant effect above this height. Wimmer and Downes (2003), on the other hand, reported a positive relationship between ring width and ring density in Norway spruce, which was associated with higher late-season rainfall. However, a negative relationship was found when faster radial growth was a result of high early-season rainfall. Besides changes linked to climatic conditions, the relationship between ring width and wood density may also be altered by the application of fertilization (Mäkinen et al. 2002) and thinning treatments (Jaakkola et al. 2005).

In black spruce, some studies have reported that ring width is not significantly correlated with wood density in mature wood (Risi and Zeller 1960; Hall 1984). Conversely, negative correlations between wood density and growth-related traits (diameter at breast height (DBH), tree height and bole volume) were detected by Zhang and Morgenstern (1995). Zhang et al. (1996) also reported moderate negative correlations between wood density and growth rate in 15-year-old black spruce trees, but these could be reversed for given provenances, environmental conditions and geographical locations. More recently, Koubaa et al. (2000) observed the combined effects of cambial age and growth rate and concluded that the negative correlation between ring width and wood density in black spruce was not significant beyond a cambial age of 25 years.

The effect of height along the stem on the radial patterns of variation in wood density has received rather less attention. In Norway spruce, Mäkinen et al. (2007) found no significant effect of height on ring density, whilst Jyske et al. (2008) found a small positive height effect (3–6 %) for a given cambial age. In black spruce, Alteyrac et al. (2005) reported a large variation in density among stem heights, but the variation was smaller in mature wood than juvenile wood.

The aim of this study was to investigate the relationships between cambial age, growth rate and height along the stem with annual wood density variation in black spruce trees. More specifically, we present the development of a model which describes the general within-stem patterns of variation in average ring density in black spruce from naturally regenerated, unmanaged stands. The resulting model was used to produce simulations of annual ring density variation in black spruce trees of different ages and growth rates.

2 Materials and method

2.1 Sampling

Trees were sampled in an area 400 km northwest of Thunder Bay, Ontario, Canada. A total of 13 stands were sampled from natural forests located in the Lake St. Joseph Ecoregion, between 51°42′ N and 52°6′ N and between 90°3′ W and 94°31′ W, at an elevation of 400 m (Table 1). Mean annual temperature varied between −1.7 and 1.0 °C, mean annual precipitation ranged between 613 and 787 mm, mean summer rainfall between 244 and 299 mm, and the mean length of the growing season between 162 and 179 days.

Table 1 Mean (SD) stand-, tree- and ring-level characteristics for each site

Each stand was pre-selected according to the following criteria: (1) species composition (dominated by black spruce), (2) stand age (>50 years old) and (3) crown closure (>55 % crown closure). Since they were unmanaged stands, no previous silvicultural treatments had been performed. Three temporary sample plots (TSP) with a plot radius of 5.64 m (100 m²) were established at each site.

Three trees were selected for destructive sampling from each TSP: (1) the tree with the largest DBH (measured 1.3 m aboveground level), (2) the tree with the median DBH and (3) the tree with the DBH closest to the average of the quadratic mean DBH and the DBH of the smallest tree of the plot, hereafter described as dominant, co-dominant and intermediate trees (Smith et al. 1997), respectively. This ensured that a large range of radial growth rates was represented in the dataset. Mean tree- and ring-level characteristics for each dominance level are shown in Table 2. After felling, the total height of each sample tree was recorded before 5-cm-thick transverse discs were cut at predetermined heights. In the dominant and co-dominant trees, discs were taken at 0.5, 1.3 and 1.75 m, the live crown base (LCB, defined as the lowest live branch) and at two equidistant points between 1.75 m and LCB. Discs at successive heights from ground level were labelled M1, M2, M3, M4, M5 and M6, respectively. In the intermediate trees, only the M2 disc (1.3 m) was sampled as discs beyond this height were generally very small. Sampling positions were adjusted to avoid branch whorls or obvious stem damage, and discs <8 cm in diameter were omitted from the study.

Table 2 Mean (SD) tree- and ring-level characteristics for each dominance class

2.2 Wood density measurements

The discs were cut longitudinally through the pith and then into radial strips with a constant thickness of 2 mm in the longitudinal direction and a width of 25 mm tangentially. The samples were placed in a conditioning chamber set at 60 % relative humidity and 20 °C until they reached a constant mass, corresponding to an equilibrium moisture content of approximately 12 %. The samples were then scanned with an X-ray beam at a resolution of 40 μm in a QTRS-01X Tree Ring Analyzer (Quintek Measurement Systems Inc., Knoxville, TN, USA). No extractions were carried out prior to scanning because the extractives content of black spruce wood is known to be very low (Lohrasebi et al. 1999).

The X-ray densitometry values were calibrated using a gravimetrically determined target density for each specimen, calculated from the mass and volume of each sample. After scanning, the transition point between successive annual rings was determined as the point of equidistance between the maximum value of the previous ring and the minimum value of the current ring. Any incomplete or false rings and rings with compression wood or branch traces were excluded from the dataset. Only data where cambial age (CA) was less than or equal to 140 years were used for the development of the models because only one tree had a higher CA. Including damaged or small-diameter discs, samples from a total of ten trees were omitted from the analysis. The final dataset therefore consisted of density data from 35,127 annual rings in 450 radial strips from 107 trees sampled at 13 sites.

Due to very slow radial growth in several samples, it proved impossible to delineate the limits of earlywood and latewood with consistent accuracy in the densitometry data. To examine the factors that influenced density in narrow rings, thin radial–transverse cuts were prepared for examination under a microscope. On each sample, the 20 rings closest to the bark were examined to determine ring width, earlywood width and latewood width. The limits between wood types were determined visually and width measurements were made using the Motic Images Plus 2.0 digital microscopy software (Marconetto 2010). Eight pairs of high- and low-density ring series were compared in this way. A series of conditions were applied to ensure that each pair was directly comparable, i.e. the two discs were from the same site, from trees of the same dominance class and age (maximum difference of 10 years), and from the same sampling height. The absence or presence of a density decline near the bark was used to determine whether discs represented either a high- or low-density ring series. The ring series were selected randomly from all those that met the specified conditions. One pair of representative ring series is shown in Fig. 1.

Fig. 1

Radial–tangential microscopy images (×16) from a pair of ring series showing a decline and no decline, respectively, in wood density near the bark. The samples were taken from dominant trees from the same site (site 6) that were of a similar age (135 and 134 years old, respectively) and from adjacent sampling heights (M5 and M4, respectively)

2.3 Model development

Due to the hierarchical structure of the data, with observations on annual rings nested within discs, themselves grouped in trees clustered within sites, nonlinear mixed-effects modelling techniques were first applied (Guilley et al. 2004; Franceschini et al. 2010; Gardiner et al. 2011). However, in our initial modelling efforts, we found that the variation among individual discs followed many distinctive pith-to-bark patterns, which could not easily be described by a single equation. To circumvent this problem, a two-stage modelling strategy was used (Achim et al. 2006; McLane et al. 2011; Duchateau et al. 2013). First, the parameters of a nonlinear equation describing the pith-to-bark profiles of annual ring density were estimated for each disc separately. The parameter values obtained in the first stage were then modelled as functions of tree- and disc-level variables using linear regression.

2.3.1 Stage 1: Fitting nonlinear model forms to individual discs

Several nonlinear model forms used to model wood density in other species (Franceschini et al. 2010; Gardiner et al. 2011) were tested, but performed poorly due to a low rate of convergence among discs. This was attributed to a decline in wood density that was observed in the outerwood of a proportion of the samples that the models could not describe. Instead, we found that the radial pattern of ring density could be more appropriately described using a Michaelis–Menten equation with the addition of two exponential terms.

$$ \operatorname{RD}={b}_0\cdot {e}^{-{b}_1\cdot \operatorname{CA}}+\frac{b_2\cdot \operatorname{CA}}{b_3+\operatorname{CA}}\kern0.8mm -\kern0.8mm {e}^{-{b}_4\cdot \operatorname{CA}}+\varepsilon $$

(1)

where RD is the mean ring density of the ring (in kilograms per cubic metre), CA is cambial age (years), and b ₀, b ₁, b ₂, b ₃ and b ₄ are empirically determined parameters. Parameter b ₀ denotes the intercept when CA approaches a theoretical value of 0; b ₁ denotes the scale parameter of the early decrease in RD. The error term ε was assumed to be normally distributed, with N (0,σ ²). The first exponential term ${b}_0\cdot {e}^{-{b}_1\cdot \operatorname{CA}}$ has a value which decreases rapidly in the first few annual rings and then gradually approaches an asymptote close to zero. The second part of the equation is a Michaelis–Menten function where b ₂ reflects the pseudo-asymptotic RD of mature wood and b ₃ denotes the rate of progression between this asymptote and a minimum value reached in the juvenile period (Auty and Achim 2008; Auty et al. 2012). In the last part of the equation, the b ₄ parameter allows for a decline in RD that can occur at high cambial ages, which can override the asymptotic value of b ₂.

The optim function in the R statistical programming environment (R Core Team 2013) was used to search iteratively for the best combination of parameters in each disc, i.e. those which minimized the sum of squared residuals. Based on a visual inspection of each radial profile, the specific conditions for the inclusion of the exponential term containing the b ₄ parameter were: (1) when a declining trend in RD was observed close to the bark in the mature wood zone (i.e. CA > 25) and (2) when the decline continued for a period of at least 20 years. Convergence was achieved for 85 % of the radial samples using the optim package. For the remaining samples, the starting values were entered manually, so that we were able to estimate parameters for the entire dataset. Model fit was assessed through a combination of visual analysis of plots of residuals against fitted values and explanatory variables, fit indices (R ²) and selected model error statistics (Parresol 1999).

$$ \operatorname{ME}=\frac{{\displaystyle \sum {y}_i\kern0.8mm -\kern0.8mm {\widehat{y}}_i}}{n} $$

(2)

$$ \left|\operatorname{ME}\right|=\frac{{\displaystyle \sum \left|{y}_i\kern0.8mm -\kern0.8mm {\widehat{y}}_i\right|}}{n} $$

(3)

$$ \operatorname{RMSE}=\sqrt{\frac{{\displaystyle \sum }{y}_i\kern0.8mm -\kern0.8mm {\widehat{y}}_i}{n}} $$

(4)

$$ \operatorname{ME}\%=\frac{100}{n}{\displaystyle \sum \frac{y_i\kern0.8mm -\kern0.8mm {\widehat{y}}_i}{{\widehat{y}}_i}} $$

(5)

$$ \left|\operatorname{ME}\right|\%=\frac{100}{n}{\displaystyle \sum \frac{\left|{y}_i\kern0.8mm -\kern0.8mm {\widehat{y}}_i\right|}{{\widehat{y}}_i}} $$

(6)

which represent the mean error, mean absolute error, root mean square error, mean percentage error and mean absolute percentage error, respectively. In these equations, y _i is the observed value, ${\widehat{y}}_i$ the predicted value, and n is the number of observations.

2.3.2 Stage 2: Linear models for each parameter of the disc-level profiles

Once the parameters of Eq. 1 had been estimated for each disc, linear mixed-effects models were fitted to each vector of parameter values in turn. The general form of the models can be expressed as

$$ {Y}_i={\boldsymbol{X}}_i\boldsymbol{\beta} +{\boldsymbol{Z}}_i{\boldsymbol{b}}_i+{\varepsilon}_i $$

(7)

where Y _i represents parameters b ₀, b ₁, b ₂, b ₃ and b ₄ in Eq. 1; X _i is the vector of fixed effects for the parameter of interest; β the associated vector of fixed-effects parameters; Z _i is a vector of random effects groups; b _i is the vector of random effects variables; and ε _i the residual error term. The random effects b _i and within-group errors ε _i are assumed to be independent for different groups and to be independent of each other for the same group and normally distributed; that is, ε _i follows N (0,σ ²) and b _i follows N (0,Ψ) for each group, Ψ being the variance–covariance matrix of random effects (Pinheiro and Bates 2000).

In order to find the most representative set of predictor variables and avoid potential issues with multicollinearity, the potential explanatory variables were divided into two categories.

• 1. Tree-level descriptors: tree age at breast height (Age, years); total height of the tree (H); and mean ring width in the mature wood (i.e. from cambial age 25 to 60) at 1.3 m (RW_1.3, in millimetres)

• 2. Disc-level descriptors: disc position along the stem (h, metres from stem base); average ring width of the first three rings near the pith (RW_p, in millimetres); mean ring width in mature wood (RW_h, in millimetres); and average growth rate (diameter or ring area) for different ranges of cambial age close to the bark (RW_b, in millimeters, or RA_b, in square millimetres). In addition, we fitted a linear regression of annual ring width against cambial age for the rings close to the bark. The slope coefficient of the regression—described by Bigler and Bugmann (2003, 2004) as ‘growth trend’ (GT, in millimeters per year)—was also screened as an input variable in the models since it is also closely related to tree vigour. Definitions of the abbreviations used for the explanatory variables in this study are given in Table 3.

  Table 3 Definitions of the abbreviations for the candidate explanatory variables used in this study

Considering the logical link between the selected variables and the behaviour of each parameter in Eq. 1, we hypothesised that RW_p was likely to be included in the models for parameters b ₀ and b ₁. Similarly, we hypothesised that RW_1.3 was related to parameter b ₂. In addition, parameter b ₄ was thought to be related to the age of the sample and the growth rate either close to the bark or in the mature wood zone. Thus, Age, RW_b or RA_b, and RW_h were candidate variables that described the competitive status of the tree at each height in the stem. Tree vigour-related variable GT was also considered.

After we fitted the mixed-effects models for parameters b ₀, b ₁, b ₂, and b ₃, a correlation analysis revealed that the parameters obtained at the disc level were not independent of each other. For those parameters, therefore, we used seemingly unrelated regression (SUR) for parameter estimation. SUR accounts for correlated error terms among different equations by estimating the parameters of the equations simultaneously (Lei et al. 2005). Additionally, since parameters b ₁ and b ₃ did not follow a Gaussian distribution, models were fitted to their natural log-transformed values.

In samples where there was no observed decline in RD near to the bark, the second exponential term in Eq. 1 was omitted from the model. Logistic regression was therefore used to predict the occurrence of the b ₄ parameter as a function of Age, RW_b, RA_b, RW_h or GT. The predicted probabilities take values between 0 and 1 and were categorized into binary outcomes using a cut-point analysis (Hein and Weiskittel 2010). The sensitivity and specificity of all potential threshold values (from 0 to 1, with an interval of 0.001) were calculated and the optimum was defined as the value giving the highest sum of sensitivity (true positives) and specificity (true negatives). Then, a linear mixed-effects model (Eq. 7) was fitted to the b ₄ values for the instances where a decline was observed. The estimated parameters were then combined with results from the logistic model so that the second exponential term of Eq. 1 was only present when the predicted value was above the cut-point threshold.

The coefficients of the individual parameter models were then substituted back into Eq. 1 for the prediction of wood density. The performance of the full model was then evaluated using the error statistics given by Eqs. 2–6, and from fit indices calculated from the fixed effects of the model. Simulated wood density profiles were then constructed using predictions from the full model. The simulations were based on the mean values of the explanatory variables in the dominant, co-dominant and intermediate diameter classes. Predictions were then plotted to give a visual representation of the correlations between the explanatory variables and wood density profiles. To show the correlations of within-tree characteristics with wood density, different percentiles of the measured data were used for selected explanatory variables, whilst the remaining tree-level variables were fixed to their mean values within each dominance class. Disc-level variables were fixed to selected values (e.g. mean quartile values) within each disc level and dominance class combination.

A potential issue with the use of an exponential function in Eq. 1 to model the decline near the bark is that predictions cannot be extrapolated beyond the range of the data since they would approach zero very quickly. Therefore, we attempted to estimate a limit beyond which the predictions of the decline are no longer valid. For all discs in which the decline was observed, we calculated a relative rate of decline for which our reference was the average RD of rings 25–60 (RD₆₀). For each disc where the decline occurred, the difference between RD₆₀ and the minimum RD (RD_min) at the bark was calculated. The 95th percentile of the observed ratios of (RD₆₀ − RD_min)/RD₆₀ was 27.2 %, which we used in the simulations as a limit of applicability of our model.

3 Results

When averaging the pith-to-bark profiles for all discs, RD was high near the pith and decreased rapidly in the first 12 growth rings, i.e. declined from 591.0 kg m⁻³ to a minimum of 473.8 kg m⁻³ between rings 10 and 15. This was followed by a slow increase until a more constant value was reached between rings 25 and 60. Based on the criteria described above, the dataset was separated into two groups. It was observed that 106 of the 450 radial samples showed a gradual decline in RD near the bark and therefore required the inclusion of parameter b ₄. In the remaining samples, there was no density decline towards the bark (Fig. 2). In the samples with a decline in ring density, annual ring width ranged from 0.1 to 3.8 mm, with a mean value of 0.8 mm, whilst RD ranged from 250.2 to 784.9 kg m⁻³ (mean = 504.4 kg m⁻³). In the samples with no decline, ring width ranged from 0.1 to 4.8 mm (mean = 1.0 mm) and RD ranged from 262.5 to 839.5 kg m⁻³, with a mean value of 499.8 kg m⁻³.

Fig. 2

Ring density vs. cambial age for the 107 black spruce sample trees. Lines represent smoothing splines fitted to the data at each sampling height. a Discs (106) with a density decline close to the bark. b Discs (344) with no density decline close to the bark

A summary of the parameter values obtained in stage 1 of the model fitting process (Eq. 1) can be found in Table 4. For the individual parameter models in stage 2, various combinations of just six explanatory variables were found to be significantly correlated with the variation of the parameters estimated using Eq. 7 (Tables 5 and 6). The occurrence of parameter b ₄ was found to be significantly related to Age and RW_h, but was not related to GT, RW_b or RA_b. A predicted value threshold of 0.104 was used as the cut-point in the logistic model, above which a decline in RD near the bark was considered to be present (sensitivity = 0.63, specificity = 0.77). Then, in the model fitted to the observed occurrences of b ₄, there were significant effects of Age and h, but no relationship with growth rate was detected (Table 6).

Table 4 Descriptive statistics for parameters b ₀, b ₁, b ₂, b ₃ and b ₄ obtained by fitting Eq. 1 to samples from individual discs

Table 5 Coefficients and associated SEs for each explanatory variable for the b ₀, b ₁, b ₂ and b ₃ parameter models and adjusted R ² values

Table 6 Coefficients and associated SEs for each explanatory variable for the b ₄ parameter model and adjusted R ² values calculated from the fixed effects

The disc-level models (stage 1) had moderate fit (mean R ² = 0.37) and were generally unbiased, despite the large variability in the radial density profiles (Fig. 3a and Table 7). When predictions were calculated from the full model, the R ² decreased to 0.17, although the mean error and mean percentage error (0.2 kg m⁻³) remained very low at 0.1 kg m⁻³ and 0.2 %, respectively (Fig. 3b and Table 7).

Fig. 3

a Model residuals of the stage 1 predictions vs. cambial age. b Residuals of the stage 2 predictions against cambial age. Percentiles of the residuals are represented by the different shades of grey. The grey dashed line represents the median value of the residuals for a given cambial age

Table 7 Fit indices and error statistics calculated from the fixed effects of the disc-level models and full model given by Eq. 1

In general, lower minimum values of RD around the 15th annual ring and higher maximum values near the bark were attained in the lower stem (Fig. 4). For a given cambial age, the RD of a slower-grown (i.e. lower value of RW_1.3) intermediate tree was generally higher than for the other dominance classes, although the decline in wood density near the bark was more likely to occur (Fig. 5). However, even the RD profiles of dominant trees could also show a declining trend beyond a certain age (not shown).

Fig. 4

Simulated radial density profiles with cambial age at different sampling heights (h) for a typical dominant tree. Tree-level variables (i.e. H, Age and RW_1.3) were set to their mean values for dominant trees, whilst disc-level variables (i.e. h, RW_p, RW_h) were set at their mean values for each disc height in a dominant tree. Lines represent the mean predictions at successive disc heights (M1–M6) from the stem base

Fig. 5

Simulated radial density profiles with cambial age at breast height (h = 1.3 m) for trees in each dominance class. The horizontal line denotes the proposed limit of application of our model. Tree-level variables (i.e. H, Age and RW_1.3) were set to their mean values within each dominance class, whilst disc-level variables (i.e. RW_p, RW_h) were set at their mean values at breast height within each dominance class

The more detailed analysis of ring series close to the bark did not reveal any relationship between the occurrence of a decline and annual ring width. However, the proportion of latewood was, on average, lower in the discs with a density decline (Fig. 6).

Fig. 6

Graph of average latewood proportion against the year of annual ring formation for the two groups of ring series that were analysed under a microscope

4 Discussion

The wood density data in our study showed the typical radial pattern found in spruces that has also been observed in other studies (Alteyrac et al. 2005; Gardiner et al. 2011). Average ring density was high near the pith and decreased rapidly outwards, followed by a slow increase until values stabilised between growth rings 25 and 60. This pattern corresponds to the type II pith-to-bark profile described by Panshin and De Zeeuw (1980).

Investigating the complex interrelationships between cambial age, growth rate and height along the stem is recognised as one of the most difficult aspects of wood density studies (e.g. Alteyrac et al. 2005; Gardiner et al. 2011). The highly variable radial patterns among samples, particularly in the juvenile wood zone, proved difficult to incorporate into a single analysis. However, modelling individual profiles in a first stage and then constructing linear models for those parameters in a second proved to be a robust solution. The choice of an equation containing three separate components allowed enough flexibility to represent the large variation in density profiles. The first two terms in the equation are capable of describing the type II density pattern previously described, whilst the final exponential term was included when it was necessary to model a decline close to the bark.

The relatively poor fit of the full model in this study could be explained by the fact that, despite its classification as a type II species, wood density profiles in black spruce are, in reality, extremely variable. This variability will inevitably increase the absolute magnitude of the model errors by increasing the deviation between the observed and predicted ring density. However, the sampling strategy in this study was designed to capture the natural variability found in black spruce unmanaged stands by sampling across a range of tree ages and sizes. By modelling the parameters of the individual density profiles, we were able to reduce these errors when compared with initial model fits to the pooled dataset. The low fit indices for the b ₁ and b ₃ parameters were most likely a result of the large variability in the initial decline near the pith.

The high density of wood adjacent to the pith in black spruce has been attributed to high earlywood density and low earlywood proportion (Alteyrac et al. 2005; Koubaa et al. 2005). Both earlywood density and latewood proportion tend to decline rapidly in the first few rings, leading to a decline in overall ring density. The wood adjacent to the pith, termed ‘flexure wood’ by Telewski (1989), is thought to be formed in response to complex loading patterns to which trees are subjected at a young age. This high density wood also has steeply oriented microfibrils that confer low stiffness, and together, these characteristics ensure that trees are flexible enough to withstand high stresses without breaking. Differences in local conditions between individual trees may lead to different initial values of density and, therefore, different rates of decline.

As expected, ring density was negatively correlated to the average ring width of the mature wood. In our analysis, the best predictions were obtained between rings 25 and 60. However, these must not be taken as precise limits. Rather, they are indicative of an overall negative relationship between radial growth rate and wood density in the part of the profile described by the second term of our equation. Accordingly, dominant trees tended to have lower density rings in the period preceding the decline than the trees in other dominance classes. This is in accordance with observations in Sitka spruce that, for a large range of ring widths, latewood width remains relatively constant between rings, so that variation in ring width is closely related to earlywood width (Brazier 1970). Similar results have been found in black spruce (Zhang et al. 1996) and jack pine (Pinus banksiana Lamb.) (Savva et al. 2010).

A gradual decrease in ring density after approximately 60 years was observed in around 24 % of the samples, which, to our knowledge, has been reported in very few studies (e.g. Schär and Schweingruber 1987). The latter study reported a similar density decline in subalpine fir (Abies lasiocarpa (Hook.) Nutt.). This was associated with a reduction in latewood density after approximately 50 years, whilst earlywood density remained relatively constant. Elsewhere, evidence of a positive correlation between growth rate and wood properties has been reported, but such correlations appear to be mainly related to water availability. For example, in young loblolly pine (Pinus taeda L.), higher ring density associated with faster radial growth was related to the increased growth of latewood in wetter growing seasons (Cregg et al. 1988). Conversely, a simultaneous decrease in ring width and wood density in Norway spruce coincided with increased early season rainfall, which resulted in the increased production of earlywood (Wimmer and Downes 2003).

In the current study, the observed decline in wood density typically spanned several years, with a faster rate of decline at higher cambial ages. For these reasons, it is highly unlikely that it was related to water availability. The decline in wood density near the bark was associated with higher cambial ages and slower growth rates in the mature wood, i.e. in the period before the onset of the decline. Comparable studies of black spruce wood density have tended to sample either young or dominant trees (e.g. Zhang et al. 1996; Alteyrac et al. 2005). Measures of mean growth rate and the slope of mean growth rate (‘growth level’ and ‘growth trend’, respectively) have been used to predict tree mortality in Norway spruce (Bigler and Bugmann 2003, 2004), although in our study, growth trend was not significantly correlated with the decline in wood density. Trees with little or no latewood—a phenomenon referred to as ‘starvation wood’ by Brazier (1977)—have also been reported growing in poor stands in Scandinavia. Examination of the microscopy samples in our study tended to confirm the existence of a similar type of wood in samples with a density decline, which contained a low proportion of latewood. It is possible that, in situations where resources are severely limited, black spruce trees can only maintain radial growth by producing lower density wood (Swenson and Enquist 2007). The observed decline in wood density may therefore imply a reduction in aboveground resource allocation. Such a change in resource allocation may help increase the longevity of old, declining trees (Robichaud and Methven 1993).

Although height along the stem had a significant influence on all the parameters of Eq. 7, its overall effect varied with cambial age. Other studies on Norway spruce have found a weak or insignificant pattern of variation in wood density with position in the stem (Mäkinen et al. 2007; Jyske et al. 2008). In the current study, wood density tended to increase with height in the stem up to a cambial age of approximately 45 years, which is in agreement with the findings of Alteyrac et al. (2005). It has been suggested that this is due to the physiological ageing of the cambium, so that juvenile wood formed higher up the stem has more mature traits than wood produced at young cambial ages lower down the stem (Larson 1969; Alteyrac et al. 2005). However, the fact that the trend was reversed beyond ring 45 suggests that more work is necessary to fully understand the factors driving the variations in wood density with height along the stem.

5 Conclusions

This study investigated wood density variations in black spruce stems growing in unmanaged stands. A two-stage modelling approach was justified since it helped separate the interrelated effects of cambial age, growth rate and height along the stem on wood density. In mature wood, radial growth rate was found to be negatively correlated with wood density. However, declining trends in annual ring density near the bark were observed, which were found to be more common in old and slow-growing trees. We hypothesise that such trends reflect a gradual reduction in tree vigour over the life of the tree.

The observed decline in wood density has implications for forest management because silvicultural interventions designed to reduce the effects of the decline by promoting vigorous growth, such as commercial thinning, may need to be considered. An appropriate rotation age may also be used as a means of avoiding the decline in wood density. In addition to any considerations about its possible effect on end-product properties, the occurrence of the decline could also have implications for aboveground biomass estimations, which could be improved by including the conditions for when a decline might occur. Future work should be conducted to investigate the possible impacts of this phenomenon on carbon allocation to other parts of the tree.