Size-density trajectories for even-aged sessile oak (Quercus petraea (Matt.) Liebl.) and common beech (Fagus sylvatica L.) stands revealing similarities and differences in the mortality process

Abstract

Key message

We studied the size-density trajectories of pure even-aged unthinned experimental sessile oak ( Quercus petraea (Matt.) Liebl.) stands in the ranges of 994–135,555 trees per hectare initial densities, observed from the ages of 5 to 38. We compared them to unthinned beech ( Fagus sylvatica L.) stands from the same experimental area. An original piecewise polynomial function was fitted to the trajectories, giving way to various applications. For each species, the initial number of trees per hectare ( N ₀ ) and the mean girth at breast height at the onset of mortality (Cg ₀ ) were parameters of the trajectory model, in addition to the parameters of the maximum size-density lines. The two former parameters (Cg _0, N ₀ ) were tied by a linear relationship, which allowed the prediction of trajectories for initial densities not included in the study data. For oak and beech, mortality onset occurred at a constant relative density (RDI), for all initial stand densities, respectively, 0.35 and 0.29. The comparison of the size-density trajectories of oak and beech allowed to establish that oak needs more space than beech for comparable mean girth, and then is less efficient than beech in its space requirements.

Context

This paper models the size-density trajectories of pure even-aged sessile oak stands, including the early development stage. It compares the oak results with those on common beech on the same site from a previous study.

Aims

A novel approach to size-density trajectories, with an original polynomial piecewise function previously used for beech stands on the same site, was satisfactorily used again as a mortality model to provide references to managers of oak forests.

Material and methods

A 38-year-old oak spacing trial, re-measured from year 5 to year 38, provided the opportunity to study the size-density trajectories of unthinned stands of this species.

Results

The fit of the piecewise polynomial function allowed us to estimate the parameters of the size-density trajectories of all stands, which were the initial number of trees per hectare (N₀) and the mean girth at breast height at the onset of mortality (Cg₀), in addition to the intercept (a) and slope (b) of the maximum size-density line. A linear relationship between Ln(N₀) and Ln(Cg₀) (where Ln is the Neperian logarithm) allowed us to reduce the number of parameters needed to fit the trajectories and made it possible to predict a size-density trajectory from any initial density not observed in the experimental stands. Moreover, this later line appeared to be parallel to the maximum size-density line, and new data allowed to establish that this was also the case for the beech stands on the same site. This parallelism feature translates to the onset of mortality occurring at the same relative density for stands of every initial density that is 0.35 for oak and 0.29 for beech.

Conclusion

Given the parameters of the maximum size-density line, a single-parameter function family could be used to predict the size-density trajectories of oak stands. The predicted trajectories have various applications in oak silviculture and growth simulators. The oak data and new data for beech stands on the same site allowed to compare the two species and draw conclusions on similitudes and differences concerning mortality and space requirements of both species.

1 Introduction

Size-density relationships provide an insight into the dynamics of undisturbed forest stands from the point of view of the correlated variations of a mean stand characteristic and the number of trees per hectare (number of trees is sometimes used as an abbreviation). Usually, the mean characteristic of the stand is the mean quadratic diameter at breast height (1.30 m) or the mean quadratic girth at breast height (mean girth is used as an abbreviation). We will consider pure undisturbed even-aged stands, only subject to intra-specific competition-induced mortality (coined “regular mortality” by Lee 1971). It has been observed that, as the mean quadratic girth at breast height (Cg) of these stands increases, the number of trees per hectare (N) does not decline until a critical value (Cg₀) of Cg, which marks the onset of competition induced tree mortality, is reached (first stage of stand development). Then, for greater values of Cg (second stage of stand development), N progressively decreases (Smith and Hann 1984, 1986), and eventually, its relative variation becomes proportional to Cg relative increment (third stage of stand development). This later relationship¹ translates to Ln(N) depending linearly on Ln(Cg), where Ln is the Neperian natural logarithm. This linear relationship was described and studied—probably for the first time—by Reineke (1933) and became a main subject of investigations on various forest species, afterwards, as attested to by many observations of forest stand data from Reineke (1933) to—for instance—Pretzsch and Biber (2005) and more recently Trouvé et al. (2017). Many other references were given in Le Goff et al. (2011) and Ningre et al. (2016a, b). Thus, the third stage of stand development is reached when stand density is at its maximum value, and the linear relationship between the maximum number of individuals per unit area and their mean size, on a log-log scale, is called the maximum size-density relationship. We do not consider a later old growth stage, due to tree decline, when overstory trees die in an irregular fashion (Oliver 1996). Maximum size-density relationships exist also for plants in general (Yoda et al. 1963; Deng et al. 2012). For forest trees, it has been used to define relative stand density measures (Reineke 1933; Curtis 1970, 1982) and contributed to establish management diagrams (Drew and Flewelling 1979; Jack and Long 1996, Lopez-Sanchez and Rodriguez-Soalleiro 2009; VanderSchaaf and Burkhart 2012; Long and Vacchiano 2014). Useful discussions of self-thinning of forest stands and maximum size-density can be found in chapter 6 of Johnson et al. (2009) and in Lhotka and Loewenstein (2008), with references on past works on the subject.

For a given forest stand, the size-density trajectory is the succession of the observed values of N and Cg, and, by extension, of the logarithms of these values, over a given time period. In this paper for oak, we used a new kind of continuous smooth piecewise polynomial function—introduced in a previous work on beech (Ningre et al. 2016a)—to model and fit the mean size-density trajectories of oak stands, expressing the number of trees (Ln(N)) as a function of the mean girth (Ln(Cg)) over the three development stages considered above. The number of trees per hectare (N₀) and the mean girth (Cg₀) at the onset of mortality were the parameters of this spline function, which allowed the explicit estimation of these two characteristics. Other uses of segmented polynomial functions to model size-density trajectories appear in VanderSchaaf and Burkhart (2008), VanderSchaaf (2010), and Cao and Dean (2008). With the function we used, all the stand characteristics at the critical development stages, namely (N₀), (Cg₀), and (Cg₁), the mean girth at which the stands reach the maximum density, could be explicitly estimated.

In this paper, we report a study of the size-density trajectories of unthinned oak stands. These oak stands are a part of the same experimental site (Lyons-la-Forêt) as the beech stands whose data were used in a previous study of size-density trajectories (Ningre et al. 2016a). The two species shared a common experimental design. The same three segments piecewise polynomial function defined in the beech study was used for the oak stands; thus, comparisons could be made between the two species. The oak stand measurements covered a 33-year period, beginning at an early development stage, and with a time interval between measurements of one to 3 years, exceptionally 4 years. These data featured a wide range of initial densities. As a result of a functional dependence between the initial oak stand density (N₀) and the mean girth (Cg₀) at the onset of mortality, all size-density trajectory equations reduced to a one-parameter function family indexed by the initial stand densities.

There are few studies on self-thinning and maximum size-density of oak stands, due to the scarcity of data needed (Johnson et al. 2009). Most are based on mixed species forest data in the USA (Williams 2003, for instance). In Europe, Barrio Anta and Álvarez González (2005) considered stand density management measures related to the stand density index of Reineke (1933), for a limited density range of pedunculate oak (Quercus robur L.) stands. For the same purpose, Fonseca et al. (2017) obtained maximum size-density relationships for cork oak stands (Quercus suber L.) with up to 3350 trees ha⁻¹. The review of Pretzsch and Biber (2005) on maximum size-density lines for European pure even-aged forest stands included oak stands with a rather limited range of densities. In France, maximum size-density lines were established for sessile oak on forest inventory data by Charru et al. (2012) and on data from unthinned experimental stands by Le Goff et al. (2011). Thus, the oak data set of the Lyons-la-Forêt long-term experiment addressed this knowledge gap.

Throughout this work, we used a relative density measure (Curtis 1970), the relative density index (RDI), to measure stand density (Drew and Flewelling 1979). For a pure even-aged stand with quadratic mean girth at breast height Cg, N trees per hectare, and mean area—or average growing space—per tree s = 1/N, the stand RDI is the ratio s_m / s, where s_m is the mean area per tree of a stand with the same quadratic mean girth at breast height and having the maximum possible number of trees per hectare (N_Max) permitted by this mean girth. Thus, we have: RDI = N / N_Max = N / (e^aCg^b), where a and b are the parameters of the maximum size-density—linear—relationship for the considered species (Ln(N_Max) = a + bLn(Cg)), respectively, the constant at the origin and the slope.

Further references and details, particularly concerning the origin and the development of forest tree stand trajectories and their applications, are to be found in the previous beech study (Ningre et al. 2016a). For this study, we used nonlinear mixed models, leading more easily to detailed statistical analysis than “frontier methods” considered by Zhang et al. (2005) and Vanclay and Sands (2009).

2 Material and methods

2.1 The Lyons-la-Forêt experiment and the oak plot survival and mortality data

2.1.1 Study site and experimental design

The data of this study came from an experimental plantation established in the spring of 1981 in the Lyons-la-Forêt National Forest in northwestern France (49° 28′ 19″ N, 1° 33′ 52″ E). This experiment was set up at the same time and in the same site as the experimental beech (Fagus sylvatica L.) plantation used to model size-density trajectories of beech stands (Ningre et al. 2016a). The soil and climate of the site were very favorable for beech and for oak. The dominant height² of oak at 36 years of age was 18.7 m, which corresponds to an excellent site index for oak (Sardin 2008).

The acorns used to regenerate the experimental site were all from the same origin (state forest of Roumare in Normandy, northwest France) and were either directly sown on the site or planted as seedlings. The oak seedlings were 1 year old at the date of plantation, in the spring of 1981, at which stage they were 25-cm tall (Colin et al. 2012). The seedlings were planted in an open area at three different densities (1333, 2667, and 5333 seedlings per hectare), and the acorns were sown at the density of 50 trees per m² (500,000 trees per hectare) (Table 1). Irregular mortality (Ningre et al. 2016a) due to transplant shock was rather low and planted seedlings were replaced after mortality, more or less successfully, in the year following plantation (Colin et al. 2012).

Table 1 Description of the Lyons-la-Forêt experiment, established in 1981

The planted plots were arranged in a randomized complete block design, covering a total area of 1.3 ha, whereas a single sown plot of 500 m² was installed close to the planted plots. The design comprised four blocks, each consisting of three rectangular plots corresponding to the actual stand densities (1333, 2667, and 5333 seedlings per hectare), the plot size depending on initial density (Table 1). Each plot was surrounded by a buffer zone of 4 to 6 m in width, depending on the plantation density. To maintain a strictly intra-specific competition in the early phases of stand development, several vegetation and pest control measures were undertaken during the experiment (see Colin et al. 2012 for details).

The dominant heights of the 4 blocks in 2003 (mean of dominant heights for all densities)—12.19, 12.48, 12.75, and 12.42 at 24 years of age for blocks 1, 2, 3, and 4 respectively—were not significantly different (ANOVA analysis), which confirmed the homogeneity of the site conditions of the four blocks. The dominant height of the sown plot in the same year was slightly lower (11.3 m), probably due to the high level of inter-tree competition.

2.1.2 Density “treatments”

Irregular mortality, after initial replacement of dead planted seedlings, was relatively low, except for the plot at planting density of 1333 trees per hectare in block 4. Then, five density treatments were defined corresponding to the planting densities and the density of sowing, and the plots with comparable initial densities were grouped together. The plots with a planting density of 1333 trees per hectare were split into two groups, 1333a and 1333b, as they presented contrasted initial densities (Table 1).

Some of the planted plots were thinned one time. For those plots, the thinning was carried out in 2003 for the plots of block 2 and block 3 at initial planting densities of 1333, 2667, and 5333 trees per hectare. Stand density was reduced to a same level—600 to 656 trees per hectare—independently of the initial density, after selecting crop trees (between 30 and 45 trees selected, depending on plot area). Only the part of the size-density trajectory before thinning was used in the modelization process, while the part of the observed density trajectories after thinning was considered for comparison with the predicted trajectories obtained from the model.

2.1.3 Size-density data

The oak stand data, quadratic mean girth (Cg) and density (N)—number of trees per hectare—cover a 33-year period from 1984 to 2017 (Ningre 2019).

Planted plots

At the time of the first measurement in 1984, a permanent sample of 43 to 50 trees with a uniform spatial distribution was selected in each plot (i.e., a total of 580 trees). From 1984 to 1997 (i.e., in years 1984, 1990, 1994, and 1997), total height and girth at breast height were regularly measured on all live trees from the original sample. The quadratic mean girth (Cg) and the density (N)—number of trees per hectare—of each plot were estimated from the tree sample data (as it was done for beech, Ningre et al. 2016a). Cg was estimated for the first time in 1990 when nearly all of the trees had reached the height of 1.30 m.

From 2000 to 2017, all live trees in each plot were measured regularly for girth at breast height (nearly each year since 2004). The total number of trees measured, before and after some plots were thinned in 2003, appears in Table 2. For the period 2000 to 2017, plot density and mean girth were obtained from the complete inventory of live trees per plot. The reduction of the number of trees measured, apart from the reduction due to thinnings between 2003 and 2004 for “thinned plots,” was due to self-thinning alone (regular mortality: see Ningre et al. 2016a).

Table 2 Number of live trees measured in the planted plots of the Lyons-la-Forêt experiment since the first complete survey in year 2000

Sown plot

At the time of the first measurement in 1987, the tree survey was made on nine 1-m² sub-plots systematically distributed in the sown area of the experiment. All trees were higher than 1.3 m. For the following surveys and up to 2006, one plot of 50 m², centered in the 500-m² sown plot, was defined (Tables 1 and 3). A last survey in 2017 was done on a larger area of nearly 170 m² including the 50-m² sub-plot.

Table 3 Number of trees measured and calculated stand density in the sown plot of the Lyons-la-Forêt experiment, since 1987

2.2 Size-density curve modeling and statistical analysis

The three segments piecewise polynomial function defined below has been successfully applied to beech and Douglas-fir stands (Ningre et al. 2016a, b). It is a continuous and smooth—i.e., continuously differentiable—function. We will refer below to the definitions given in the introduction. The size-density trajectory equation for a given stand is defined as follows with the function f:

$$ {\displaystyle \begin{array}{l}\mathrm{Ln}(N)=f\left({\mathrm{Cg}}_0,{N}_0,a,b,\mathrm{Cg}\right)\\ {}\kern3em =\Big\{\begin{array}{l}\mathrm{Ln}\left({N}_0\right)\kern9em \mathrm{if}\ 0<\mathrm{Cg}\le {\mathrm{Cg}}_0\\ {}p+q\mathrm{Ln}\left(\mathrm{Cg}\right)+r\mathrm{Ln}{\left(\mathrm{Cg}\right)}^2\kern2.25em \mathrm{if}\ {\mathrm{Cg}}_0<\mathrm{Cg}\le {\mathrm{Cg}}_1\\ {}a+b\mathrm{Ln}\left(\mathrm{Cg}\right)\kern7em \mathrm{if}\ \mathrm{Cg}>{\mathrm{Cg}}_1\end{array}\end{array}} $$

(1)

The parameters p, q, and r appear in the continuity and derivability constraints which solve in Cg₀, N₀, a, and b to give the following values:

$$ \Big\{{\displaystyle \begin{array}{l}p=\mathrm{Ln}\left({N}_0\right)+\frac{b^2\mathrm{Ln}{\left({\mathrm{Cg}}_0\right)}^2}{4\mathrm{Ln}\left({\mathrm{RDI}}_0\right)}\\ {}q=-\frac{b^2\mathrm{Ln}\left({\mathrm{Cg}}_0\right)}{2\mathrm{Ln}\left({\mathrm{RDI}}_0\right)}\\ {}r=\frac{b^2}{4\mathrm{Ln}\left({\mathrm{RDI}}_0\right)}\end{array}} $$

(2)

RDI₀ being the stand relative density just before the occurrence of tree mortality, when the number of trees is still N₀, and the mean girth is Cg₀ (Ln (RDI₀) = Ln(N₀) − (a + bLn (Cg₀))). From the continuity and derivability constraints again, the mean girth at breast height (Cg₁) of the stand when it reaches the maximum density depends only on the previous parameters a, b, Cg₀, and N₀ and is given by the following equation:

$$ \mathrm{Ln}\left({\mathrm{Cg}}_1\right)=\frac{2\left(\mathrm{Ln}\left({N}_0\right)-a\right)-b\mathrm{Ln}\left({\mathrm{Cg}}_0\right)}{b} $$

(3)

Adequately, Eq. (1) is qualitatively suitable to depict the development stages of a pure even-aged forest stand reflected by its size-density trajectory, which is continuously evolving from an initial plateau to a linear maximum size-density relationship on a double logarithmic scale (Fig. 1).

Fig. 1

Graph of the size-density trajectory model, in log-log scales, defined by a piecewise polynomial function with Eqs. (1) and (2). Mortality onset occurs at (Cg₀, N₀) and the reach of maximum density at Cg₁. This function models the three stages of forest stand development in terms of size and density. Equation (4) is the statistical model used to fit this function

Equation (1) was fit to the oak data using a mixed effect model with b as a fixed parameter common to all treatments. As for the preceding beech study (Ningre et al. 2016a), a graphical analysis suggested that the parameter a had a mean value common to all treatments with an added normal random effect due to stand to stand variability and additional hypotheses. The parameters depending on the treatments were Cg₀ and N₀, the first as a fixed parameter, and the second having a normal random component due to stand effects. The resulting statistical equation was the following:

$$ {\displaystyle \begin{array}{l}\mathrm{Ln}\left({N}_{ij k}\right)=f\left(C{g}_{0i},{N}_{0i}+{N}_{0 ij},a+{a}_{ij},b,{\mathrm{Cg}}_{ij k}\right)+{\varepsilon}_{ij k}\\ {}i=1,\dots, n\\ {}j=1,\dots, {m}_i\\ {}k=1,\dots, {l}_{ij}\end{array}} $$

(4)

where the indices i, j, and k are for treatment, stand, and year, respectively; n is the number of treatments; m_i is the number of stands for treatment i; and l_ij is the number of years when stand j was measured in treatment i. Statistical treatments of the oak data were based on the following hypotheses of the mixed-effects models (Pinheiro and Bates 2000). The within-group residual errors ε_ijk are normally distributed random variables with mean 0; the random effects a_ij and N_0ij are such that the vectors (a_ij, N_0ij)^T are normal variables with mean 0 and variance-covariance matrix Ψ, independently distributed from one treatment or one stand to the other (i.e., cov (a_ij, a_i’j’) = cov(N_0ij, N_0i’j’) = cov (a_ij, N_0i’j’) = 0, if i ≠ i’ or j ≠ j’). Graphical analysis showed that for the residual errors (ε_ijk), a standard error (σ_ij) specific to each stand was needed. This later hypothesis could be handled with the preceding ones by the R Environment for Statistical Computing, with the mixed model package nlme that was extensively used during the study (Pinheiro and Bates 2000; Pinheiro et al. 2014; R Core Team 2014). These hypotheses were subjected to validation.

Furthermore, based on this fit, graphical analysis showed that the data were reasonably consistent with the assumption of aligned inflection points (Ln (Cg₀), Ln(N₀)), or, in other words, that there were two parameters a₁ and b₁ satisfying the following equations:

$$ \mathrm{Ln}\left({N}_{0i}\right)={a}_1+{b}_1\mathrm{Ln}\left({\mathrm{Cg}}_{0i}\right),i=1,\dots, n $$

(5)

Consequently, Eq. (4) was fit to the data with the constraints of Eq. (5). For this purpose, the mean girth (Cg) was replaced in Eq. (4) by its value as a function of the initial number of trees per hectare (N₀) obtained from Eq. (5), and Eqs. (1) and (2) were passed to the nlme fitting package as a R function.

3 Results

The fit of Eq. (4)—the statistical version of Eq. (1)—to the oak data clearly showed a linear relationship (Fig. 2) between the estimated values of Ln(N₀) and Ln(Cg₀) (Ln(N₀)=11.81–1.31Ln(Cg₀), R² = 0.99). Consequently, Eq. (4) was fit to the data with the constraints of Eq. (5) (that is, mean girth (Cg) was replaced in Eq. (4) by its value as a function of the initial number of trees per hectare (N₀) obtained from Eq. (5)). We assumed that there was a random effect due to treatment for the parameter a₁, without interaction with the stand effects. This random effect was needed to obtain estimates of a₁ and b₁ close to the values of the parameters of the line fitted to the inflection points estimated by the unconstrained fit of Eq. (4). Resulting from this fit, Fig. 3 shows the fitted trajectories for each treatment with the aligned inflection points and their limiting maximum size-density line, together with the observed stand data. The standardized—corrected to allow for plot discrepancies—residuals of the fit relative to the fitted values did not show time dependence, which is consistent with the graphics presentation of their empirical autocorrelation function (R ACF function, not shown) nor did they show heteroscedasticity (Fig. 4). Normality of these residuals was assessed with a QQ plot (Fig. 5). The normality of the random effects was also assessed by a QQ plot (except for the random effect associated with a₁, because there were too few values). The estimated values of the parameters a₁ and b₁ (a₁ = 11.87 and b₁ = − 1.331) were very close to the values obtained with the—least square—linear regression fit to the estimated values of the inflection points resulting from the fit of the unconstrained equation (i.e., Eq. (4) alone). Table 4 gives the parameters and statistics of the fit of Eq. (4) with the constraints of Eq. (5) (referred hereafter as trajectory model 1).

Fig. 2

Alignment of the inflection points of the oak stand size-density trajectories obtained by fitting Eq. (4), without constraint, to the oak stand data of the Lyons-la-Forêt experiment

Fig. 3

The oak stand size-density trajectories observed data of the Lyons-la-Forêt experiment and the model obtained by fitting to them Eq. (4), with the constraints of Eq. (5) (alignment of the inflection points). The maximum size-density line and the line where lie the inflection points are represented

Fig. 4

Standardized residuals versus the fitted values for the fit of Eq. (4) trajectory model, with the alignment of the inflection points constraints (Eq. (5)) to the oak stand data of the Lyons-la-Forêt experiment. Neither bias, heteroscedasticity, nor correlation appear. Each panel is identified with the format treatment/stand

Fig. 5

Normal (QQ) plot of the standardized residuals of the fit of Eq. (4) to the oak data of the Lyons-la-Forêt experiment with the linear constraints (5) on the inflection points. The linearity of the relationships provides evidence for the normality of the within-group errors (ε_ijk in Eq. (4))

Table 4 Parameter estimates and statistics of Eq. (4) fit to the oak experimental stand data with the constraints of Eq. (5)

Then, the attempt was made to fit Eq. (4) with the constraints of Eq. (5) and the hypothesis of the alignment of the inflection points parallel to the maximum size-density line (Fig. 6; Table 5), i.e., with the additional constraint (b₁ = b). The likelihood ratio test (Table 6) failed to reject this hypothesis at the 1% probability level (Wilks χ² test). The diagnostic plots of this fit (trajectory model 2) showed as good results as those previously performed without the hypothesis of parallelism, and the adequacy of the fitted model was fully assessed by displaying the fitted and observed values on the same plot (Fig. 7). Relying on this latest fit, we derived stand characteristics for the various treatments at the critical stand development stages corresponding to the onsets of mortality and maximum RDI (Table 7). The ages appearing in this table were obtained thanks to a mean girth-age nonlinear relationship fitted to the experimental data, with the plot initial densities (trees ha⁻¹) as a covariate (1.93 cm residual standard deviation for the mean girth estimation).

Fig. 6

The oak stand size-density trajectories obtained by fitting Eq. (4) to the data of the Lyons-la-Forêt experiment, with alignment of the inflection points parallel to the maximum size-density line (Eq. (5), with the additional constraint b₁ = b)

Table 5 Parameter estimates and statistics of Eq. (4) fit to the oak experimental stand data with the constraints of Eq. (5), and where, in addition, b₁ = b

Fig. 7

Observed and predicted size-density trajectories of the oak stands of the Lyons-la-Forêt experiment obtained with the fit of the model defined by Eq. (4) and the constraint of alignment of the inflection points parallel to the maximum size-density line. The plain line is the prediction at the treatment level which does not include the random effects due to stand-to-stand variations accounted for by the (complete) fit represented by the dashed line. The figure does not show any lack of fit

Table 6 ANOVA for the fit of Eq. (4) to the oak data with the constraints of Eq. (5), (1) without or (2) with the additional constraint b₁ = b

Table 7 Stand characteristics estimated by Eq. (4) with the added constraints of Eq. (5) and b₁ = b, when fitted to the oak data

Finally, the size density trajectory data of the stands that were eventually thinned to 621 trees per hectare were plotted together with trajectory model 2 for stands of initial density equal to the stand densities to which the thinned stands were lowered. Figure 8 allows comparing the observed trajectories with the trajectories predicted by this model.

Fig. 8

Oak stand size density data of the Lyons-la-Forêt experiment before and after thinning to a mean density of 621 trees per ha, compared to the trajectory predicted by Eq. (4) and the constraint of alignment of the inflection points parallel to the maximum size-density line, with initial density set to the residual density of the stands after thinning (621 trees per ha)

4 Discussion

4.1 Statistical treatment of the experimental data

One might think that a statistical test could be needed to assess the alignment of the inflection points. But, on one hand, in this context of a mixed model, as both random and fixed effects are involved, there is no statistical test currently available (Pinheiro and Bates 2000, pp. 82–91 § 2.4 Hypothesis tests and confidence intervals) to make the decision. The difference between the AIC values of the two models were large, but it did not seem advisable to resort to AIC in this case, as the aligned inflection point model was not embedded in the unconstraint model (because the structures of random effects in the two models were different), and the use of AIC is controversial in this case (Burnham and Anderson (1998), Ripley (2004)). On the other hand, we cannot think that there is no monotonic relationship between the initial densities of stands, and the mean girths at the onset of mortality. And, as it appears in Fig. 2, it does not seem that an alternative to a log-log linear relationship is possible.

Concerning the parallelism of the inflection point alignment to the maximum-size density line, the significance level of 1% probability was retained. This level was more prone to the risk of accepting a false hypothesis than the 5% level. But, Pinheiro and Bates (2000) warned that the statistical test used in this case is anticonservative (and they add “sometimes quite badly so”), that is, the resulting p value is underestimated. Moreover, no lack of fit of the model is appearing on Fig. 7, which does not cast doubt on the adequacy of the model.

4.2 Modelization of the size-density trajectories

The Lyons-la-Forêt density trial allowed us to analyze and fit the size-density trajectories of oak for a wide range of initial densities, as we did previously for beech on the same site (Ningre et al. 2016a). Similarly, the considered period covers the early development stage before competition-induced mortality appears, up to the stages where stands reach and follow the maximum size-density line (except for the lowest initial densities) with a particular attention to the onset of mortality (inflection point of the size-density trajectories). Moreover, we revised the size-density trajectories of beech because new data were available (see Appendix 1). Thus, the Lyons-la-Forêt density trial gave us the opportunity to compare the characteristics of the size-density trajectories of oak and beech on the same site conditions.

The size-density trajectory model, previously established for beech (Ningre et al. 2016a) and used for Douglas-fir (Ningre et al. 2016b), worked also for oak in the Lyons-la-Forêt trial. The piecewise polynomial function defined (Eq. (4)) proved adequate to represent the size-density trajectories of unthinned oak stands of various initial densities, from the early stages of development preceding the onset of mortality up to the maximum relative density. The components (Cg₀, N₀) of the inflection point of a given trajectory are parameters of this function, together with the specific parameters (a, b) of the maximum size-density line. The logarithms of these two components were linearly correlated for the oak stands of this study. Thus, as previously seen for beech and Douglas-fir, their size-density trajectories were uniquely determined by the initial stand density N₀.

4.3 Maximum size-density line

The constant and slope parameters of the maximum size-density line (a = 13.733 and b = − 1.652 respectively) obtained for oak in this study, with the constrained fit (trajectory model 2), were close to the ones previously obtained in a study dedicated to the comparison of maximum size-density lines of several broadleaved species on a larger scale in France (13.531 and − 1.566 respectively, Le Goff et al. 2011). These later values were not contained in the confidence intervals established for oak parameters in this study [13.56, 13.90] and [− 1.69, − 1.61] for the constant and slope parameters respectively, but close to the lower limit for the intercept and to the upper limit for the slope. The values obtained for oak by Charru et al. (2012) differed much more (14.87 and − 1.911 respectively), being largely outside the confidence intervals of the parameters obtained in this study: this may be due to the uncertainty of the maximum size-density line established by Charru et al. (2012) from data contained in a restricted diameter interval (15 to 30 cm) which, moreover, did not correspond to the interval of diameters of this study (1 to 25 cm). The estimation of the constant and slope parameters for oak in the Lyons-la-Forêt experiment will possibly be improved, as the plots at low initial planting densities (1333 trees per ha) will reach the maximum size-density line, that is when Cg and N will reach the values Cg₁ and N₁ respectively (see Table 7).

The constant and slope parameters of the maximum size-density line obtained for beech in the same experiment (Table 8; Appendix 1) compared relatively well with those of oak: the oak intercept (13.733) was contained in the beech confidence interval of the intercept coefficient [13.3, 13.9], while the slope parameter of oak (− 1.652) appeared slightly higher than that of beech (− 1.512) and outside the confidence interval of beech slope parameter (− 1.58, − 1.45) but not very far. Then, the maximum relative density of oak and beech in the Lyons-la-Forêt experiment occurred at about the same point for small mean girths, but for oak appeared to be lower for higher mean girths, the maximum size-density line of oak presenting then a slightly steeper slope than that of beech (Fig. 9a).

Table 8 Likelihood ratio test for the fit of Eq. (4) to the beech size-density data with the constraints of Eq. (5) and b1 = b

Fig. 9

Comparison of size-density trajectory features of oak and beech in the Lyons-la-Forêt experiment. a Characteristic size-density lines. b Average growing space of trees at maximum density and at mortality onset. c Size-density trajectories—log scale—for three initial densities (N₀ = 1000, 5000, 40,000 trees ha⁻¹) together with characteristic size-density lines. d Stand characteristics—age A₀ and dominant height H₀—at mortality onset, in relation to initial stand density (N₀)

The maximum size-density line (Ln(N) = a + bLn(Cg)) can be expressed in terms of the average growing space per tree (s_m), as s_m is proportional to 1/N with a factor depending on the surface unit. The variation of growing space with mean tree diameter (Fig. 9b) is such that oak occupies much space than beech, and its space need increases with age much more than it does for beech. This is in line with the results obtained by Pretzsch and Biber (2005) in the range of diameters of this study. However, for higher diameters, Pretzsch and Biber found that beech space requirements tended to meet and even exceed those of oak, in relation with the steeper slope of the maximum size-density line of beech as compared to that of oak observed in his case: this differed from our results (Fig. 9b). The maximum size-density lines established for oak and beech for a large range of diameters by different authors confirm the higher growing space requirement of oak (Rivoire and Le Moguedec 2012; Le Goff et al. 2011). Beech appears then more efficient than oak in its space requirements.

4.4 Mortality onset and initial density

The locus of the inflection points of the stand density trajectories corresponding to oak mortality onset was a line whose slope was not considered significantly different from the slope of the maximum size-density line. In other words, mortality began to occur at a same relative density (RDI₀ = 0.35), independently of initial stand density (Fig. 9a; Table 7).

As new inventory data were available, a re-evaluation of the locus of the inflection points of the size-density trajectories of beech in the Lyons-la-Forêt experiment was done (Table 9; Appendix 1). As for oak, it was found to be a line parallel to the maximum size-density line, mortality onset occurring at a constant RDI value equal to 0.29 (Table 10; Appendix 1). As a consequence, mortality of oak appeared to occur at a slightly higher relative density than that of beech (0.35 compared to 0.29). The occurrence of mortality onset at a constant RDI was observed for other species in foreign countries (see references in Ningre et al. 2016a): that was especially the case of red alder (Alnus rubra Bong.) (Puettmann et al. 1993).

Table 9 Parameter estimates and statistics of the fit of Eq. (4) to the data of the beech stands of the Lyons-la-Forêt experiment, with the constraints of Eq. (5) where b₁ = b (see text for details)

Table 10 Stand characteristics for the beech stands of the Lyons-la-Forêt experiment estimated by Eq. (4) with the constraints of Eq. (5) where b₁ = b. RDI, mean girth (Cg₀), and age at mortality onset

The trajectory Eq. (1) may be parameterized with RDI₀ (equation system (2) without further development of RDI₀). Fitted to the observed oak and beech data, this parameterized form of Eq. (1), with the constraints of Eq. (5), allowed estimation of the approximate 95% confidence intervals of RDI₀ for the two species, namely [0.23, 0.47] for the value 0.35 obtained for oak, and [0.24, 0.34] for the value 0.29 obtained for beech (the calculated confidence intervals are comparable with that of red alder [0.31, 0.5] obtained by Puettmann et al. (1993)).

The size-density line at mortality onset (Ln(N)=a₁+bLn(Cg)) can be expressed, as the maximum size-density line, in terms of average growing space per tree (s₀). At mortality onset, s₀ followed what was observed for the reach of maximum stand density (Fig. 9b). That is, the mean growing space (s₀) was higher for oak than for beech, which meant that less trees could stay alive in the case of oak and again that oak was less efficient than beech in its space requirements.

4.5 Size-density trajectories

To allow the comparison of oak and beech, the fitted parameters of the size-density trajectories of these species in the Lyons-la-Forêt experiment were used to predict size-density trajectories for a larger interval of initial densities.

4.5.1 Mortality charts

While mortality onset occurred at a relative initial density higher for oak, mean stand dimensions (Cg₀) at this stage appeared to be quite comparable for oak and beech. For a given initial density, however, after mortality onset, oak presented a higher mortality rate than beech, reaching its maximum “more rapidly,” that is within a narrower interval of mean girths, as the maximum size-density line was lower for oak than for beech, for the major part of stand development (Fig. 9c).

The age (A₀) at the onset of mortality was obtained with a relationship fitted to the age-Cg data, with stand initial density as a covariate (see “Results”). For oak, A₀ increased from ages 8 to 17 when initial density decreased from 40,000 to 1000 trees per ha (Fig. 9d). At the same time, dominant height (H₀) increased from about 3 up to 8 m. For beech, A₀ increased from about 11 to 20 years of age and H₀ from 4 to 9.5 m (Fig. 9d). So, while starting at about the same Cg₀, mortality onset happened earlier and at a lower dominant height for oak than for beech, for a given initial stand density.

These stand characteristics at mortality onset, depending on initial stand density, may help forest managers to carry out the first thinning before mortality begins. In the conditions of the Lyons-la-Forêt experiment, these thinnings should occur earlier for oak than for beech and at a lower dominant height despite a slightly greater height growth of oak (Fig. 9d).

4.5.2 Mortality simulation

As for beech (Ningre et al. 2016a), in the absence of an available individual mortality model, the size-density model established here for oak in the Lyons-la-Forêt experiment should be adequate to model mortality prior to any management treatment in stand growth simulators as “Fagacées” (Le Moguedec and Dhôte 2012).

However, unlike what was observed for thinned beech stands, for oak, the observed onset of mortality after thinning occurred several years after the predicted onset (about 8 years after, relying on the predicted trajectory of stands of initial density equal to the residual density after thinning; Fig. 8). This may be due to the type of thinning performed in the thinned oak plots at Lyons-la-Forêt where density was drastically reduced by removing all suppressed trees and a high proportion of codominant and dominant trees, leaving only about 600 dominant crop trees per ha when the initial stand density was comprised between 1200 and 5333 trees per ha. By removing all suppressed trees—those trees more prone to die—and by reducing strongly stand density, it took several years before the first trees among the remaining ones had their crowns sufficiently reduced by competition to die. Thus, mortality after thinning occurred when the stand density reached an RDI value between 0.6 and 0.75, which is relatively far from the RDI at the mortality onset observed in unthinned plots (0.35) and outside the confidence interval of this value.

4.5.3 Stand management

To help with the management of oak stands, the RDI corresponding to 10% mortality in number of trees was calculated (the 10% mortality level appears in forest management as a relatively low but nevertheless significant level of mortality). It appeared that this level is reached at a constant RDI value equal to 0.62, independently of the initial stand density N₀ (see Appendix 2).

Thus, in order to maintain the mortality of oak under 10%, the forest manager should keep stand density in the area comprised between the mortality onset line (RDI = 0.35) and the 10% mortality line (RDI = 0.62) in the (Ln(Cg)–Ln(N)) coordinate system. This is in accordance with the objective of maximizing net stand growth (gross stand growth minus mortality), which was observed at a RDI around 0.5 for sessile oak (Trouvé et al. 2019). Maintaining mortality under 10% can only be achieved with relatively heavy thinnings reducing stand density below the 0.35 RDI line, so as the density of the thinned stands remains between the 0.35 and 0.62 RDI lines after crown canopy closes again. This guideline might work fine for the first thinning. For subsequent thinnings, mortality might be delayed, as observed in our study, making the upper bound of the 10% mortality line following successive thinnings a point of uncertainty given data currently available. However, in our case, thinnings were very specific and differ from those currently applied in forest management.

Other mortality lines could be considered by choosing different mortality ratios (other than 10%) after calculating the corresponding RDI values (Table 11; Appendix 2). It can be noticed that when stands reach the maximum size-density line (RDI = 1), mortality already reaches a relatively high level: in the case of oak, 65% of trees have already died since mortality onset. This shows how important it is to include size-density trajectory models in stand growth simulators and to be careful when selecting temporary plots in order to establish the maximum size-density line of a given species (Le Goff et al. 2011).

Table 11 Values of relative stand density RDI corresponding to increasing cumulative mortality rates, between the onset of mortality and the advent of stand maximum relative density

Fig. 10

New fitted size-density trajectories (in logarithmic scales) of the beech plots of the Lyons-la-Forêt experiment with inflection points aligned parallel to the maximum size-density line (observed size-density data coming from the two blocks of the Lyons-la-Forêt beech design are also shown: see Ningre et al. 2016a)

Such characteristic density lines are at the basis of the construction of density management diagrams—called SDMDs—which are commonly used in the USA since Drew and Flewelling (1979) and more recently in Europe (Sales Luis and Fonseca 2004; Penner et al. 2006; Lopez-Sanchez and Rodriguez-Soalleiro 2009; Patricio and Nunes 2017; Minoche et al. 2017). It can be mentioned here that the characteristic size-density lines represented in these SDMD diagrams are often drawn somewhat arbitrarily parallel to the maximum size-density line, whereas the parallelism of these characteristic size-density lines for oak and beech could be really established in this study (Fig. 10; Appendix 2).

Fig. 11

Stand density management diagram (SDMD) for oak in the Lyons-la-Forêt experiment. Three parallel size-density lines are represented: the maximum size-density line (RDI = 1), the line of mortality onset (RDI = 0.35), and the line of 10% cumulative mortality rate (RDI = 0.62). Also represented on the figure is the zone where cumulative mortality is under 10% (shaded area)

5 Conclusion

The piecewise polynomial function used to describe the size-density trajectories of beech (Ningre et al. 2016a) and Douglas-fir stands (Ningre et al. 2016b) proved again useful in describing the size-density trajectory of oak in the Lyons-la-Forêt experiment and suited to a wide range of initial stand densities (from about 1000 to 140,000 trees per ha).

Our results show that the inflection points of the oak trajectories are aligned in parallel to the maximum size-density line. Then, mortality of oak starts at a constant RDI value, independently of initial stand density. It appeared to be also the case for beech, which was not established at first (Ningre et al. 2016a) but could gain evidence thanks to new additional data.

The maximum size-density line obtained for oak was close to the one established earlier at a larger scale (Le Goff et al. 2011). The comparison with that of beech in the same site conditions allowed to establish that oak needs more space than beech for comparable mean size. Then, beech appears more efficient than oak in its space requirements.

As a result of this study, a SDMD diagram composed of three useful lines could be established for the management of even-aged oak stands: one indicating the conditions of maximum stand density and the two others being relative to mortality. The first one corresponds to the beginning of competition-induced mortality and the other one to the reach of a 10% mortality level (in number of trees per ha). This level is considered as a limit not to be exceeded to reduce the production losses. Incorporating size-density trajectory equations in stand growth simulators would also improve the prediction of tree mortality. However, the conditions of renewed mortality after thinning operations still need to be better established.

It would be advisable, looking forward, to obtain additional data, more particularly for the low density oak plots of the Lyons-la-Forêt experiment, in order to better estimate the parameters of the maximum size-density line. More data should also be necessary to establish the onset of mortality after thinning for oak, as it appears from the Lyons-la-Forêt data that thinned oak stands do not behave as unthinned oak stands of the same initial density, as it was observed for beech (Ningre et al. 2016a).

Finally, the segmented model fitted here to the oak size-density trajectories, and already fitted to beech and Douglas-fir size-density trajectories, seems to be easily adapted to other species if appropriate data are available.

Appendices

Appendix 1. Update of the size-density trajectories of beech in the Lyons-la-Forêt experiment

As new inventory data were available for beech in the Lyons-la-Forêt experiment (Ningre 2019), the size-density trajectories previously established (Ningre et al. 2016a) were fitted again. This new fit, with a larger data set, allowed a better comparison with the size-density trajectories of oak. The density treatments considered in the updated data set were the same as those described in Ningre et al. (2016a), and the previous size density data were only supplemented with the new measurements made in 2016. We used the same piecewise polynomial function already used for beech in Ningre et al. (2016a), oak in this paper (Eqs. (1) and (2)), and for Douglas-fir (Ningre et al. 2016b).

The statistical model used to fit the new size-density relationships was the same as for Eq. (4), with two differences. There was no heteroscedasticity of the residual errors to be corrected, and the graphics of the empirical autocorrelation function of these errors showed an autoregressive correlation of order 1 that needed to be accounted for. These statistical hypotheses were validated with the diagnostic plots of the residuals.

As already observed for beech (Ningre et al. 2016a), and for oak in this study, with the new beech data, the inflection points—(Cg_0, N₀) on the log-log scale—of the size-density trajectories obtained by the fit of Eq. (4) still lined up. Contrarily to what happened with oak for the fit of Eq. (4) with the constraints of Eq. (5), there was no need to include a random effect on the parameter a₁ to obtain estimates of a₁ and b₁ close to the values of the—estimated—parameters of the line fitted to the inflection points estimated by the unconstrained fit of Eq. (4). As there was no random effect associated with the parameter a₁, Eq. (4) with the constraints of Eq. (5) and the additional constraint (b₁ = b) was a model embedded in the model of Eq. (4) alone, and passing from one model to the other involved only fixed effects. These conditions allowed performing a likelihood ratio test to compare these two models (Pinheiro and Bates 2000). No significant difference—at the 5% probability level—was found between the unconstrained and—doubly—constrained models (Table 8). The results of the fit of Eq. (4)—including the statistical hypotheses modified as stated for the beech case—with the constraints of Eq. (5) and the additional hypothesis (b₁ = b) are given in Table 9. The plots of the observed data and fitted curves of the constrained model did not show any lack of fit.

With this new fitting of the size-density trajectories of beech in the Lyons-la-Forêt experiment (Fig. 11), mortality onset took place at a constant relative density (RDI₀ = 0.29; Table 10), to be compared to the results found previously, that is, RDI values ranging from 0.31 to 0.58 for stand densities ranging from 625 to 40,000 trees per ha. Then, mortality appeared to start at a lower relative density than expected before, particularly for stands at high initial densities. However, the change in Cg₀ is very limited (3.4 cm compared to 4 cm for treatment 40,000, the highest initial density, and 61.6 cm, compared to 57 cm for treatment 625a, the lowest density). The change in stand characteristics at the onset of maximum relative density appears greater, particularly for stands of higher initial densities: mean girth Cg₁ and stem number N₁ are half the previous values for treatment 40,000 (Table 10 in this paper and Table 6 in Ningre et al. 2016a).

Appendix 2. Property of the constant RDI mortality lines of oak

To better control mortality in unthinned oak stands, it may be of interest to investigate a connection between their cumulated mortalities and RDI.

It appeared that from the onset of mortality (RDI₀ = 0.35) up to the advent of maximum relative density (RDI = 1), the relative number of surviving trees (N_t/N₀) of a given stand, where N_t is the number of trees per hectare observed at some time t, depended only on the relative density at this time, regardless of the initial density N₀ of the stand. The proof is sketched as follows.

Let Cg_t be the mean girth of a stand at time t, with relative number of surviving trees N_t /N₀, and let us consider the following logarithms:

$$ \mathrm{Ln}\left({N}_0\right)-\mathrm{Ln}\left({N}_t\right)=\Delta \mathrm{Ln}\left({N}_0\right) $$

(6)

Let Cg₀—as already specified—be the mean girth of the stand at the mortality onset. Using Eq. (5) with b₁ = b, it results that, for the nonlinear part of the size-density trajectory, described by Eqs. (1) and (2), the following equation holds:

$$ \mathrm{Ln}\left({\mathrm{Cg}}_0\right)-\mathrm{Ln}\left({\mathrm{Cg}}_t\right)=\frac{2}{b}{\left[\left(a-{a}_1\right)\Delta \mathrm{Ln}\left({N}_0\right)\right]}^{\frac{1}{2}} $$

(7)

Using again Eq. (5) with b₁ = b, and Eq. (7), and substituting Ln (Cg₀) for its value obtained by solving Eq. (6), we obtain the following:

$$ \mathrm{Ln}\left({N}_t\right)=c+b\mathrm{Ln}\left({\mathrm{Cg}}_t\right) $$

(8)

where $c={a}_1-\Delta \mathrm{Ln}\left({N}_0\right)+2\ {\left[\left(a-{a}_1\right)\Delta \mathrm{Ln}\left({N}_0\right)\right]}^{\frac{1}{2.}}$

The relative density RDI_t of a stand at the time t when natural mortality has reduced the stand density to N_t trees per hectare is RDI_t = N_t / N_Max,t where N_Max,t is the maximum number of trees per hectare of a stand of mean girth Cg_t. The logarithmic transformation of RDI_t gives the following:

$$ \mathrm{Ln}\left({\mathrm{RDI}}_t\right)=\mathrm{Ln}\left({N}_t\right)-\mathrm{Ln}\left({N}_{\operatorname{Max},t}\right) $$

(9)

As:

$$ \mathrm{Ln}\left({N}_{\operatorname{Max},t}\right)=a+b\mathrm{Ln}\left({\mathrm{Cg}}_t\right) $$

(10)

Using (8), (9) and (10), we obtain the following:

$$ \mathrm{Ln}\left({\mathrm{RDI}}_t\right)=c-a $$

That is:

$$ {\mathrm{RDI}}_t={e}^{c-a} $$

(11)

Then, the stands with a given relative number of surviving trees (N_t /N₀) have the same RDI value (Eq. (11)). The locus of the points corresponding to a given mortality rate for the trajectories of various initial stand densities is a line parallel to the maximum size-density line. In this study, we considered a cumulative mortality rate amounting to 10%, corresponding to a ratio of surviving trees of 0.9. For this value, the RDI is equal to 0.62 (Fig. 11). Other rates of mortality could be retained. The corresponding RDI values were calculated for increasing mortality rates (Table 11).