Wood is a renewable resource that is used as a raw material in the construction, furniture-making, paper pulp manufacturing and other industries. Many of the uses rely on information being available about wood properties. Given the wide range of timber species (Zanne et al., 2009) and their inter- and intra-specific differences (Zobel & Van Buijtenen, 2012), detailed study of wood properties of different species is essential for optimising use of this raw material.

Characterisation of wood is complicated by the anisotropic and hygroscopic nature of the material, in which basic variables (e.g. weight and dimensions) vary depending on environmental conditions. Indeed, the variations in the dimensions derived from changes in air humidity limit the end uses of wood (Zobel & Van Buijtenen, 2012; Glass & Zelinka, 2021).

Studies of wood properties should include analysis of density, which is a good indicator of wood quality, because of its positive relationship with wood shrinkage and also its positive effect on hardness and mechanical strength (Zobel & Van Buijtenen, 2012; Senalik & Farber, 2021). Density is also important as a proxy variable affecting hydraulic conductance phenomena: as density increases, the wood capacitance tends to decrease (Savi et al., 2017) and embolism resistance tends to increase (Li et al., 2016). Wood density is also related to the ecological strategies of woody plants, by determining the amount of carbon invested per unit of stem volume (Martínez-Cabrera et al., 2009).

The dimensional stability of wood is of particular interest regarding the use of wood as a structural and construction material, and it also plays an important role in furniture-making. Shrinkage and swelling are limiting factors for many end purposes of wood, as distorted or shortened beams or swollen parquet may strongly reduce the quality of the installed product.

Volumetric shrinkage and, by extension, shrinkage in all three main dimensions of wood (longitudinal, radial and tangential) are known to be related to wood density (Schulgasser & Witztum, 2015). As wood shrinkage is evaluated using not replicable methods (due to sorption hysteresis in wood) (Hou et al., 2022), indirect methods that predict shrinkage variables on the basis of other easily measurable variables, such as density, are very useful. Predictive models for estimating shrinkage variables have been constructed for a range of different species of hardwoods (e.g. Riesco Muñoz & Rodríguez Campos, 2010). Also, predictive models for intra-tree, inter-tree and inter-plot variations have previously been described (e.g. Zobel & Van Buijtenen, 2012; Glass & Zelinka, 2021; De Lima et al., 2022).

Linear mixed effects models are useful for predictive purposes beyond determining the simple linear relationship between two or more variables. Linear mixed effects models are constructed by adding one or several random terms to the equation with the aim of systematising the random error inherent in linear regression. This avoids two common problems in biological data: lack of independence of the measurements, and hierarchical structures within the data (Calama & Montero, 2004; Winter, 2013). Thus, linear mixed effects models can be used to determine the relationship between shrinkage and density in terms of goodness of fit, as has been found for some softwoods (e.g. Gapare et al., 2008; Grekin & Verkasalo, 2010).

The laurel (Laurus nobilis L.) belongs to the Lauraceae family and grows in Spain and other Mediterranean countries. The species is a remnant of the laurel forests that covered the Mediterranean Region during periods of more humid climate. In the Pliocene Era, laurel forests were displaced by the more drought-tolerant sclerophyllous plants that grow today (Woodward, 2009). Laurel grows usually on the fringe of deciduous forests, bordering agricultural areas or in coastal areas (Fig. 1), and it is more abundant in forests with Mediterranean traits (Pemán García et al., 2006).

Laurel commonly occurs as tree forms in mature deciduous stands (Ruiz de la Torre, 1991). The Cortegada Island laurel forest (Galicia, NW Spain) is one of the largest laurel forests in Europe, exhibiting complete and continuous crown cover, with some individuals reaching heights of more than 26 m after a natural regeneration process lasting 80 years (Lamas & Rozas, 2007).

Despite the importance of laurel and its wide area of distribution, descriptions of laurel wood are scarce. As far as we know, the only studies related to the physical properties of this wood species are those by Pulgar Lorenzo & Riesco Muñoz (2018) and Riesco Muñoz & Pulgar Lorenzo (2022). According to these studies, laurel wood is homogeneous, moderately heavy, compact and rather hard, and checks sometimes appear in the wood during drying. Laurel wood is much cheaper than other domestic hardwoods and it is used to craft marquetry objects and in carving and turnery.

The hypothesis of this study is that the dimensional stability of laurel wood in response to changes in moisture content can be predicted from density. The main objectives of the study were to construct models that use density to predict shrinkage variables and to determine how variations in the density-shrinkage relationship are affected by factors such as individual tree and position in the tree stem.

Seventeen individual Laurus nobilis L. trees were selected and felled in a plot (latitude: 43.624552, longitude: -7.893340, WGS 84; altitude: 197 MASL) located in Galicia (NW Spain). The sample was formed by trees of average age, 32 years (range, 24-41 years), breast height diameter, 21 cm (range, 14-29 cm), and total height, 8.6 m (range, 5.1-13.1 m).

Cross sections of the trunk (discs or slices) were obtained from the base of each felled laurel, at heights of 1 and 2 m in the stem and sometimes higher (when the stem diameter was large enough: diameter over bark higher than 10 cm). The slices, of approximate thickness 15 cm, were sawn to produce two types of defect-free prismatic specimens with longitudinal (axial), radial and tangential target dimensions of respectively 40 mm × 20 mm × 20 mm and 20 mm × 40 mm × 40 mm (Figs. 1 and 2).

Analysis of small defect-free specimens enabled detection of relationships between wood variables. Such relationships would not be easily detected in wood exhibiting defects or unusual features, which would add undesirable sources of variation (Riesco Muñoz & Remacha Gete, 2011).

The mass and the three main dimensions of the specimens were measured at three successive wood moisture contents: saturation, at approximately 12% wood moisture (equilibrium moisture content in the laboratory environment), and at 0% moisture (oven-dry wood). The mass of each specimen was measured to an accuracy of 0.01 g, and the volume was calculated (assuming each specimen to be a perfect prism) by multiplying the value of each dimension measured between centers of opposite faces with a micrometer (for dimensions close to 20 mm) and a calliper (for dimensions close to 40 mm) to an accuracy of 0.01 mm. The density was determined according to the Spanish standards UNE 56528:1978 (IRANOR, 1978a) and UNE 56531:1977 (IRANOR, 1977a), and according to the harmonised standard UNE-EN 408:2011 (CEN, 2019b). The moisture content was determined according to the harmonised standard UNE-EN 13183-1:2002 (CEN, 2019a). The air-dry density (d ₁₂) was expressed relative to a moisture content of 12%, by applying the correction equation included in standard UNE 56531:1977 (IRANOR, 1977a). The oven-dry density (d ₀) was calculated by dividing the mass of the specimens with zero moisture content (after an oven-drying period) by the corresponding volume. The basic density (bd) was calculated as the ratio between the oven-dry mass and the saturated volume of the specimens.

The shrinkage variables were determined according to the Spanish standards UNE 56528:1978 (IRANOR, 1978a) and UNE 56533:1977 (IRANOR, 1977b). The linear shrinkage in the three main grain directions (longitudinal, radial and tangential) and the volumetric shrinkage were calculated as the difference in dimension at saturation moisture content and in the oven-dry state and expressed as percentages of the dimension in the oven-dry state. Longitudinal shrinkage was only determined in the 40 mm × 20 mm × 20 mm specimens (712 specimens), to detect the small variations inherent in the axial direction on the longest dimension. The low degree of longitudinal shrinkage, together with deformation of the specimens during drying, frequently yields errors in estimation of this parameter. Similarly, radial and tangential shrinkage were only measured in the 20 mm × 40 mm × 40 mm specimens (246 specimens), to detect variations in transverse direction on the longest dimensions. The 40 mm testing length is established in the standard UNE 56533:1977 (IRANOR, 1977b).

The width of annual growth rings (rw) was calculated for each saturated specimen as the dimension of the longest radial segment in the cross section divided by the number of complete growth rings in the segment and averaging the measurements made on both faces.

The average number of specimens sampled per tree was 56, with a minimum of 14 specimens in a particular tree.

The descriptive statistics of the physical properties were obtained after validating and refining the test results.

Correlation techniques and graphical summaries were used to describe the relationships between variables. Those variables indicated to be related by significant values of the Pearson’s correlation coefficient (r) were included in linear regression models based on the ordinary least squares method. The model selection criteria involved calculating the adjusted coefficient of determination (R² _adj ) and the residual standard error (RSE). Results were considered significant at p < 0.05 and highly significant at p < 0.01.

The error term in linear or linearizable models accounts for the deviations from predicted values due to external factors not related to the covariates. In order to reduce the variance and the error term, a random component in the form of a factor or qualitative variable was added to the systemic part of the linear mixed effects model (i.e. fixed effects, in this case the intercept and the parameter estimate for the variable oven-dry density). The random component or factor was added in order to consider the possible reasons for the high variability in wood grown under similar conditions. Thus, while the fixed effects model covered the predictable overall influence in the response, the random effects model splits the data in order to establish some structure in the randomness, thus reducing the variance of the error.

Some linear mixed effects models have been constructed using the regression estimates as a starting point, to overcome problems related to autocorrelation and to minimise the variance in the error term inherent in fixed effects regression models, with the aim of explaining some of the inter- and intra-tree variation in stem properties (Bolat et al., 2023; Saygili & Kahriman, 2023).

In the present study, three factors were initially considered random effects: the tree of origin of the wood (“tree”), height in the stem (portion or “slice”) and radial position (“radial”). Tree included 17 levels, one for each felled tree. Slice was derived from the three different heights in the stem (0.15, 1 and 2 m) where the slices were sampled. Finally, radial included three transversal positions in the stem (inner, middle, and peripheral). Specimens sampled above a height of 2 m in the stem were excluded from fitting in order to achieve greater homogeneity in the sample.

The nature of the remaining variables must be considered before proceeding with further calculations. As it is not possible for a single slice to originate from different trees, the factors were nested, or rather slice was nested in tree. Accordingly, the level of hierarchy was not the same for the variables. This affected the way the random effects were ordered in the linear mixed effects model. Thus, with these possible random effects, the behaviour of the linear mixed effects model was tested by considering several scenarios:

In all these possibilities, the coefficient of the variable oven-dry density (d ₀) remained as the only fixed effect.

The expressions of the mixed models tested using the lme4 R package notation (Bates et al., 2014) were as follows:

The first model (1), shown in Fig. 3, was the simple linear regression of the volumetric shrinkage (vs) against the oven-dry density (d ₀); the whole expression being vs = a + b*d ₀. The second (2) and third (3) models included a random effect term (tree and slice, respectively), with varying intercept and constant slope. The expression “(1|tree)” denotes the tree variable as a random effect, with variable intercept and constant slope (represented by the “1” in “1|tree”). A more complex relationship between the two random effects was incorporated in model (4). Thus, although the slices were taken at the same heights in the 17 trees, each triplet of slices was taken from a specific individual tree, so that the variable slice was considered nested within the variable tree. Therefore, the model accounted for the variability caused by the tree factor, and the variability produced by the slice factor, nested in tree. This was expressed by the addend “(1|tree/slice)”, equivalent to “(1|tree) + (1|tree:slice)”, where “tree:slice” accounted for the interaction between the two factors, which was also considered a factor. Finally, both random intercept and random slope (the random slope introduced by d ₀ in the expression “d ₀|tree/slice”) were added to the previous model, which had random intercept and a constant slope, to produce model (5).

In the model selection process, the values of the residual standard error (RSE), Akaike information criterion (AIC) and Bayesian information criterion (BIC) were considered, along with the goodness of fit in the graphical representation of the models.

The Likelihood Ratio Test was used to assess the significance of the fixed effect (d ₀) in the fitted models.

Statistical calculations and graphical representations were carried out in R (Core Team R, 2013) by using the RStudio environment (RStudio Team, 2015) and several base packages as well as the broom (Robinson, 2017), ggplot2 (Wickham, 2016) and lme4 (Bates et al., 2014) packages.

One anatomical trait (ring width), three density variables and four shrinkage variables (shrinkage in volume and in the three main directions of wood) were measured in the sample. The main statistics of the studied physical properties are shown in Table 1. The number of specimens complied with the Spanish standard UNE 56528:1978 (IRANOR, 1978a) in relation to sample size required for testing the physical properties of wood on small defect-free specimens.

The partial correlation matrix (Table 2) shows the variables with potentially valuable correlations (for modelling purposes), from among the desirable predictor variables (ring width and densities, which are easy to obtain) and the four shrinkage variables. The most notable finding in the correlation matrix was the linear relationship between the oven-dry density and the volumetric shrinkage (r = 0.70**; R ² = 0.49**), in addition to the obvious relationships between density-related variables.

Therefore, the correlations between the tested variables (listed in Table 1) indicated a possible relationship between volumetric shrinkage (vs) and oven-dry density (d ₀). Plotting those variables (Fig. 3) confirmed the positive relationship. Therefore, the following models were considered:

The explanatory power of the different proposed linear or linearizable models for the relationship oven-dry density vs. volumetric shrinkage is shown in Table 3. The increase in adjusted coefficient of determination and the reduction in residual standard error was not large enough in linearizable models, i.e. equations (6) to (12), from the simplest linear model, equation (1). For the sake of simplicity, and as all linear models tested (volumetric shrinkage vs. oven-dry density) were very similar in terms of their explanatory power, model (1) was selected to describe the statistical relationship between volumetric shrinkage and oven-dry density. According to model (1) (Fig. 3), the estimated volumetric shrinkage is high beyond an oven-dry density of 569 kg m^-3, and the values were above this level in 99% of the specimens tested.

The first step in selecting random terms for the linear mixed model was to apply analysis of variance (Table 4) to check for differences in the volumetric shrinkage variable caused by the factors tree, height in the stem and radial position. In the present study, use of inter- and intra-tree variation models was possible due to the nature of the sample. Further research is required in order to analyse the variation in different plot locations.

Grouping tree and slice produced significantly different groups, although not in the case of radial position, for which a high p-value (0.504) was obtained (Table 4). The radial position factor was discarded as having a possible random effect as it had a negligible influence on the response variable. This may be explained by the very uneven sample size in each group: very small numbers of specimens were assigned to pith and bark groups relative to the number assigned to the middle stem section position.

The unavoidable limitation derived from the sample size and sample diversity under analysis logically did not allow extrapolation of the specific models obtained in this work to all existing populations of laurel. However, the observed patterns of variation (significant variation in the models between trees, the absence of significant, able-to-be-modelled variation with depth in the stem and marked variation along the stem) were very prominent and repeated in all sample trees tested. Therefore, it is possible that similar trends will be repeated in other stands in the age range considered here.

The comparative statistics for the five linear mixed effects models are shown in Table 5, in which the best results appeared when the variable tree was included: models (2), (4) and (5). For comparison of models with random slope and constant slope, model (5) usually provided a better fit than model (4), and model (5) was therefore selected. Model (5) is a linear mixed effects model that explains the relationship between volumetric shrinkage and oven-dry density, taking into account the variation carried by the factors tree and height in the stem, with both random slopes and intercepts (RSE = 1.039; AIC = 2977.980; BIC = 3021.764). Finally, Table 6 presents the summary of the chosen model, i.e. the mixed effects model (5).

In Fig. 4 several plots related to the residuals of model (5) are presented. In the residuals vs. fitted plot, the residuals appear to be fairly randomly scattered around the zero line. It suggests that the model does not have major issues with non-linearity. The spread of the residuals seems relatively constant across the range of fitted values, although there might be a slight indication of increased spread at both the lower and higher ends of the fitted values. However, it is not very pronounced. There are a few residuals that are relatively far from the zero line, indicating potential outliers or influential points. It is important to check these points to detect if they affect the model significantly.

As for the qqplot, it indicates that while the residuals in the center of the distribution might be approximately normal, there are departures from normality in the tails (both lower and upper). This suggests the presence of outliers or heavy tails, meaning that the residuals have more extreme values than what would be expected under a normal distribution. The scale-location plot suggests that the residuals do not have constant variance across all levels of fitted values. There is some evidence of heteroscedasticity, especially at the lower end of the fitted values, where the spread is larger.

The graphical representation of model (5) showed a good fit (Fig. 5). The coefficients calculated for each level of the random effects terms (tree, and slice nested in tree), and the resulting fixed effects (intercept and the coefficient of d ₀) were all negative, ranging from -3.199 to -1.892 for the intercept, and all were positive, in the range 0.028, 0.034 for the coefficient of d ₀.

The Likelihood Ratio Test assessed the significance of the fixed effect d ₀ in the new model by comparing the following two alternatives:

The ANOVA (Table 7) indicated a highly significant difference (p < 0.001) between the models, so that the contribution of the variable d ₀ remained important in the proposed mixed effects model.

The mean width of the growth rings (Table 1) indicates that laurel is a slow growing species, according to the Redington classification (1999). The air-dry density (at 12% moisture content) falls in the range of the moderately heavy hardwoods (510-700 kg m^-3) following the UNE 56540:1978 interpretation criterion (IRANOR, 1978b). The values of the three density variables tested were very homogeneous, with coefficients of variation of around 6-7%, and the relative error in the estimation of density variables was, as a consequence, always much lower than 5%, the upper limit recommended in standard UNE 56528:1978 (IRANOR, 1978a) (Table 1). The low coefficients of variation of the density variables were lower than expected for hardwood species with porous growth rings (according to the database of species-specific coefficients of variation in wood density provided by Wagenführ & Wagenführ, 2022).

A negative and highly significant linear correlation between the average ring width and density of the specimen was observed, although the correlation coefficients were negligible (Table 2). The high significance was attributed to the large sample size (n = 920). The fact that the density was homogeneous and the ring width was heterogeneous in the sample explains the weak relationship between density and ring width at stand level (Table 2).

Physical and anatomical properties of wood partly depend on the cambial age (measured as the number of rings from the pith), and the relationship between wood density and ring width occurs as an indirect effect of cambial age. Cambium age should therefore be used as a predictor covariable to improve the accuracy of ring width vs. density models (Adamopoulos et al., 2010). However, this would be impractical, because cambium age is not easy to measure by counting rings on sawn pieces. On the other hand, height in the stem and ring width can be readily obtained as predictor variables because the length of the merchantable bole is managed in harvesting operations and mean ring width can be obtained if the stand age is known.

The wood tends to shrink severely and to produce drying cracks, according to the interpretation criteria in standard UNE 56540:1978 (IRANOR, 1978b). The longitudinal shrinkage was highly variable, which is characteristic of immature trees, in which juvenile wood is abundant (Zobel & Van Buijtenen, 2012; Senalik & Farber, 2021). The transverse stability of wood (radial and tangential shrinkage) was moderate, with relatively low variability. The relative errors (Table 1) were lower than the threshold proposed in standard UNE 56528:1978 (IRANOR, 1978a), except for the variable related to longitudinal stability (because of its high coefficient of variation).

The noticeable correlation between oven-dry density and volumetric shrinkage (Table 2) yielded higher r values for some trees, considering the subsample corresponding to each individual tree, although the r values for almost one third of the trees were lower than for the whole sample. This is why development of individual-tree models was not pursued. As the sample population originated from the same plot, no models were constructed to test between-stand variation, pending further research.

The models with the nested variable slice in tree (Table 5) yielded the lowest values of RSE, AIC and BIC, and models (4) and (5) were therefore the most suitable for describing the relationship between the variables involved. Of these models, the nested variables model, equation (4), with both random intercept and slope had a higher RSE, higher AIC and a slightly lower BIC.

The model finally selected, model (5), was a linear mixed effects model that explained the relationship between the two selected quantitative variables (volumetric shrinkage and oven-dry density), taking into account the variation caused by the qualitative variables tree and height in the stem, and with both random slopes and intercepts.

Once the linear mixed effects model was constructed, the next aim was to determine whether the fixed effect for the variable d ₀ was significant. In other words, the objective was to determine whether the oven-dry density was an explanatory variable for the volumetric shrinkage behaviour, while considering the variation due to the random effects. The significant difference (p < 0.001) between the models with and without the fixed effect (Table 7) confirmed the value of including d ₀ in the linear mixed model. Thus, calculation of a linear mixed model such as that proposed here is recommended, including individual tree and height in the stem as factors, and oven-dry density as an independent variable.

The moderately high density and volumetric instability of laurel wood may represent technological restrictions regarding the industrial use of this raw material, although the uniformity of these wood physical properties across all trees at stand level is a positive trait. A linear mixed effects model with individual tree and height in the stem as factors and oven-dry density as an independent explanatory variable can accurately predict the volumetric shrinkage of laurel wood. Estimation of the expected shrinkage of the stems is helpful for selecting the most stable logs for industrial purposes. The model considers the variation due to the factors tree and height in the stem, with both random slopes and intercepts. The pith-to-bark variation was negligible in the sample of laurel stems analysed. The accuracy of the model is limited to young trees, which were the only type of trees sampled because of the rarity of larger and older individuals in the study area. Further research is required to analyse the inter-plot variation and the performance of mixed effects models in samples of older trees.