Modeling size-density trajectories for even-aged beech (Fagus silvatica L.) stands in France

Abstract

Key message

We studied the size-density trajectories of pure even-aged unthinned beech stands in the ranges of 625–40,000 trees per hectare initial densities and of 12–33 years of age. A new piecewise polynomial function family was fitted to the trajectories, giving way to various applications. Initial number of stems per hectare (N ₀) and 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 line. The two former parameters were tied by a linear relationship, which allowed the prediction of trajectories not considered in this study. Furthermore, the generic trajectory equation fitted the trajectories of thinned stands not used in the estimation of the parameters.

Context

This paper models the size-density trajectories of pure even-aged beech stands, including the early development stage, which is not as well documented as are the later stages.

Aims

The work reported in this paper concerns the development of a novel approach to size-density trajectories, considered as a mortality model to provide references to managers of beech forests.

Material and methods

A 33-year-old beech spacing trial beginning at 12 years of age provided the opportunity to study the size-density trajectories of unthinned stands of this species. The beech data helped us to develop a new piecewise function to model these trajectories. The model we chose was a polynomial segment smoothly joining two linear functions.

Results

The fits of this model allowed us to estimate the parameters of the size-density trajectories of all stands, which were the quadratic mean girths at mortality onset and at maximum density. A linear relationship between these characteristics allowed us to reduce the number of parameters needed to fit the trajectories and made it possible to predict a stand trajectory from any initial density not observed in the experimental stands.

Conclusion

A single-parameter function family could be used to fit the size-density trajectories of beech stands. The predicted trajectories have various applications in beech silviculture and growth simulators.

1 Introduction

Tree competition can be defined as the sharing of the limited resources of a stand between the trees. For a given stand, competition has two correlated effects. These effects are mortality and a reduction of the individual growth of the surviving trees. Size-density relationships quantify this correlation, with mortality expressed as the number of surviving trees per unit area and the impact of competition on growth expressed as the mean individual size of the trees. The result is that size-density relationships can be expressed by equations that predict the number of stems per hectare as a function of stand mean sizes. They express how mortality due to overcrowding is intermingled with the effects of competition on growth. These relationships provide a convenient way to study and describe mortality in undisturbed stands. Size-density trajectories can be used on their own to establish management rules (Drew and Flewelling 1979), as well as to simulate mortality when used jointly with a growth model.

The maximum size-density line of a tree species almost always appears as a linear logarithmic relationship between the number of stems per unit area and a mean stand characteristic such as the mean diameter at breast height or the mean volume. In most ecological studies, the works of Yoda et al. (1963) are often the earliest ones credited for the use of these relationships. Nevertheless, Reineke (1933) earlier observed a close linear logarithmic relationship between mean stand diameter at breast height (DBH) and the maximum number of live standing trees in undisturbed even-aged stands and recommended its use as a means to define a stand density index (SDI). Later, Curtis (1970) stressed the connection between most density measurements, including SDI. Thus, from the mid-1980s to the present, size-density relationships have been extensively investigated, and their use has been regularly reported in the forestry literature to model self-thinning, to measure stand density, and to provide insight into forest stand dynamics (Turnblom and Burk 2000). Smith and Hann (1984, 1986) were probably the first authors to consider the nonlinear part of size-density relationships, which is of particular importance to model mortality of the early development stage of stands.

After a period of irregular mortality, i.e., mortality due to various undetermined causes, following tree planting, a forest stand is subject to regular mortality (Lee 1971), also referred to as competition-induced mortality. At this stage of development, it takes more or less time—depending on the initial density—for this competition-induced mortality to appear in the stand. Stand density is not necessarily at its maximum value when mortality begins. Competition for site resources often begins before competition-induced mortality. After tree mortality has appeared, stand density still gradually increases until stands with the same mean girth at breast height—or other mean characteristics—and more stems per hectare rarefy to an extreme degree. This is considered as the maximum density stage.

Let Ln be the Neperian logarithm, Cg the quadratic mean girth at breast height, and N the number of live stems per hectare.¹ When the stand is located in the (Ln(Cg), Ln(N)) coordinate system, it follows a trajectory that reflects the process described above. The size-density trajectory is a continuous curve comprised of three parts. Starting with a zero slope for a range of Cg values (it is a horizontal line segment on this range), the slope of the trajectory gradually decreases, i.e., the trajectory becomes a concave curve, and its slope eventually reaches a minimum value that does not change with increasing values of Cg. This last part of the trajectory is a line with a constant negative slope, representing a stage of development commonly known as the maximum size-density relationship.

The research reported here originated from the need to model the size-density trajectories of unthinned beech stands of various initial densities in the two-dimensional (Ln(Cg), Ln(N)) coordinate system. Instead of the SDI (Reineke 1933; VanderSchaaf and Burkhart 2012), we used the relative density index (RDI) to measure stand density (Curtis 1970, 1982; Dhôte et al. 2000). For a stand with quadratic mean girth at breast height Cg and a mean area per tree s, we may consider the mean area per tree s ₀ in a stand with the same quadratic mean girth at breast height and with the maximum possible number of stems per hectare. The RDI of the given stand is the ratio s ₀/s. For a pure even-aged stand with mean girth at breast height Cg and N the number of stems per hectare, RDI = N / ((e ^a)(Cg ^b)), where a and b are the parameters of the maximum size-density relationship for the considered species (Curtis 1970, 1982). The SDI and the RDI differ by a constant factor. For unthinned stands at the last stage of stand development, RDI is maximum, equal to 1, and the relationship between Ln(Cg) and Ln(N) is linear.

There is no single analytical expression that makes it possible to exactly fit the complete size-density trajectory of a stand from an early development stage up to the stage of maximum density when the size-density relationship is a straight line. Analytical expressions only make it possible to approximate the whole size-density trajectory with the asymptotic behavior of a curve (Smith and Hann 1984), most often turning a size-density relationship upside down (Puettmann et al. 1993; Tang et al. 1994). At contrast, segmented functions have the advantage to model explicitly the nonlinear and linear parts of the relationship. Different examples can be found in the literature: Cao and Dean (2008) used a three-segment polynomial model and a reduced form of it as a two-segment model to fit the size-density trajectories of permanent plots. However, this model did not include the maximum RDI stage (which, incidentally, resulted in trajectories that crossed each other). VanderSchaaf and Burkhart (2008, 2010) used a four-segment polynomial model including the linear part of the trajectory at maximum RDI.

2 Materials and methods

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

2.1.1 Study site and experimental design

All of the data were taken from an experimental plantation established in the spring of 1981 in the Lyons-la-Forêt National Forest in northwestern France (Colin et al. 2012; Ningre and Colin 2007). The study site was located on flat terrain, at an altitude of 190 m, with a homogeneous leached brown soil with variable levels of mottling in the B horizon, and developed in clay with flints at depths ranging from 45 to 65 cm. The climate was oceanic with slight continental tendencies; the average annual temperature was 10.5 °C, while the average annual precipitation was 883 mm, as quantified by the nearest meteorological station for the 1992–2000 period. The dominant height at 33 years of age for the site was 18.5 m, corresponding to stands of good fertility (Pilard-Landeau and Simon 2008).

Beech (Fagus silvatica L.) seedlings were 3 years old, all from the same local origin (state forest of Arques in Normandy, northwest France). The seedlings were planted in an open area at six different densities, thus covering a wide range of stand densities (625, 2500, 4444, 5000, 10,000, and 40,000 seedlings per hectare).

The experiment covered a total area of 3 ha and was set up as a randomized complete block design, with partially replicated treatments within blocks. It was comprised of two blocks (Fig. 1), each consisting of nine rectangular plots corresponding to the six different abovementioned stand densities, the densities of 625, 2500, and 4444 being replicated twice in each block. The dominant heights of the two blocks in 2010 (mean of dominant heights for all densities except the lowest one)—18.67 and 18.35 m at 33 years of age—did not significantly differ (ANOVA analysis), which confirmed the homogeneity of the site conditions of the two blocks.

Fig. 1

Observed size-density data—number of stems per hectare (N) as function of quadratic mean girth at breast height (Cg), in logarithm transformations—from the varying planting densities—625 to 40,000 trees per hectare—of the Lyons-la-Forêt beech (Fagus silvatica L.) experimental stands and fitted trajectories based on Eq. (5). The inflection points of the trajectories are located in such a way as to make it possible to postulate a linear relationship between initial density (N ₀) and mean girth at breast height (Cg ₀) when mortality first occurs in the stands

Each plot was surrounded by a buffer zone of 3 to 8 m in width, depending on the plantation density. The plot areas—buffer zone excluded—ranged from 1440 m² for the lowest density to 600 m² for the highest density. The plots were mechanically brushed until the age of 14 years after planting to maintain a strictly intraspecific competition in the early phases of stand development.

2.1.2 Density “treatments”

Irregular mortality (Fridman and Ståhl 2001; Lee 1971) due to transplant shock was generally quite low (except for one plot), with mortality less than an average of 14 %. Subsequently, 17 of the 18 experimental plots, with mortality ratios ranging from 2 to 23 %, were retained for the study. The actual initial densities observed in these plots, i.e., the number of live seedlings per hectare after transplant shock, ranged from 472 to 38,000 trees per hectare (Table 1). Six density treatments were defined—“625a,” “625b,” “2500,” “4444,” “10,000,” and “40,000”—corresponding to planting densities, and plots with comparable initial densities were grouped together. The plots with planting densities of 5000 and 4444 trees per hectare, which presented relatively close initial densities, were grouped together in treatment 4444, and conversely, the plots with a planting density of 625 trees per hectare were split into two groups, 625a and 625b, since they presented contrasting initial densities (Table 1).

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

A single thinning (thinning from below) was performed in 1995 for three plots in block #2 at initial densities of 4444, 10,000, and 40,000 trees per hectare (Table 1) in order to simulate a cleaning operation of varying intensity. In these plots, stand density was reduced to the mean density reached in 1995 by the plots in treatment 2500, i.e., approximately 2000 trees per hectare. Except for these three thinned plots, the evolution of the initial densities after the initial period of irregular mortality was dependent on self-thinning alone (or competition-induced mortality), also referred to as regular mortality (Lee 1971; Oliver and Larson 1996, p. 78).

2.1.3 Size-density data

At the time of the first measurement in 1984, a permanent sample of 50 regularly spaced trees, i.e., with a uniform spatial distribution, was selected in each plot. From 1984 to 1995, total height and girth at breast height were regularly measured between four and six times, i.e., at 1–3-year intervals, on all live trees from the original sample, depending on the plot. The quadratic mean girth (Cg) and the density—number of stems per hectare—of each plot were estimated from the tree sample data. The plot density for a given year was calculated as the product of the actual initial density and the percentage of all live trees in the initial sample for the considered year. The mean girth (Cg) was estimated as the quadratic mean girth of the live trees in the sample, with 0 for the girth of trees smaller than 1.30 m. Cg was estimated for the first time in 1989 when nearly all of the trees had reached the height of 1.30 m. At that time, regular mortality was very low.

After 1995, all live trees in each plot were measured for girth at breast height from four to ten times depending on the plot. The number of trees measured was 2693 at the first survey after 1995, decreasing to 1986 in 2013. For that period, plot density and mean girth were obtained from the complete inventory of live trees per plot.

2.2 Size-density curve modeling and statistical analysis

To work with a smaller number of parameters, we used a three-segment polynomial model including the initial plateau of a trajectory preceding the onset of mortality and the last linear part after reaching maximum density. The initial number of stems per hectare and the inflection point at the shift from a horizontal line to a concave (polynomial) curve were parameters of this model. Data analysis then made it possible to postulate an alignment of the inflection points of the beech stand trajectories which was used to reduce the number of parameters needed to fit the trajectories.

Let Cg ₀ be the quadratic mean girth at breast height of a stand at the onset of regular mortality and N ₀ its number of stems per hectare prior to this onset. N ₀ will be referred to as the initial stand density. Let Cg ₁ be the stand quadratic mean girth at breast height when the maximum RDI is reached (since maximum RDI is reached after the onset of mortality, we always have Cg ₀ ≤ Cg ₁). Let a and b be the parameters of the maximum size-density relationship for the species considered. In the (Ln(Cg), Ln(N)) coordinate system, the size-density trajectory of an unthinned stand is composed of three parts. For 0 ≤ Cg ≤ Cg ₀, the number of stems per hectare is constant (and equal to N ₀). It then decreases with an inflection of the trajectory curving downwards at the point (Ln(Cg ₀), Ln(N ₀)). Finally, for Cg ₁ ≤ Cg, the trajectory is again linear and follows the maximum size-density relationship Ln(N) = a + b Ln(Cg) (Fig. 1).

We postulated that the part of the trajectory located between the occurrence of tree mortality and the beginning of maximum density could be approached by a second-degree polynomial with coefficients p, q, and r, in the equation, y = p + qx + rx ², where x = Ln(Cg) and y = Ln(N). If such a piecewise trajectory has a continuous and smooth equation (i.e., continuous, with a continuous derivative), Cg ₁ is fully determined by Cg ₀ and N ₀ as follows:

$$ Ln\left(C{g}_1\right)=\frac{2\left(Ln\left({N}_0\right)-a\right)-bLn\left(C{g}_0\right)}{b} $$

(1)

For this trajectory, the evolution of the number of stems per hectare, from the early sapling stage to any stage after maximum RDI has been reached, is modeled by the following family of continuous smoothed function f of Cg, depending only on the four parameters Cg ₀, N ₀, a, and b, as illustrated below:

$$ Ln(N)=f\left(C{g}_0,{N}_0,a,b,Cg\right)=\left\{\begin{array}{llll}\hfill & Ln\left({N}_0\right)\hfill & \hfill & \mathrm{if}\ 0<Cg\le C{g}_0\hfill \\ {}p+\hfill & qLn(Cg)\hfill & +r{\left(Ln(Cg)\right)}^2\hfill & \mathrm{if}\kern0.5em C{g}_0<Cg\le C{g}_1\hfill \\ {}a+\hfill & bLn(Cg)\hfill & \hfill & \mathrm{if}\ Cg>C{g}_1\hfill \end{array}\right. $$

(2)

Let RDI₀ be the stand RDI at the point (Ln(Cg ₀), Ln(N ₀)) of inflection of the trajectory. On the basis of the definition of the RDI, it follows that

$$ Ln\left(RD{I}_0\right)=Ln\left(\frac{N_0}{e^aC{g_0}^b}\right)=Ln\left({N}_0\right)-a-bLn\left(C{g}_0\right) $$

(3)

Equation (2) only depends on the parameters Cg ₀, N ₀, a, and b because the continuity constraints on this equation and its derivative imply the following relationships:

$$ \left\{\begin{array}{ccc}\hfill p=\hfill & \hfill Ln\left({N}_0\right)\hfill & \hfill +\frac{b^2Ln{\left(C{g}_0\right)}^2}{4Ln\left(RD{I}_0\right)}\hfill \\ {}\hfill q=\hfill & \hfill -\frac{b^2Ln\left(C{g}_0\right)}{2Ln\left(RD{I}_0\right)}\hfill & \hfill \hfill \\ {}\hfill r=\hfill & \hfill \frac{b^2}{4Ln\left(RD{I}_0\right)}\hfill & \hfill \hfill \end{array}\right. $$

(4)

Equation (2) was fitted to the beech data using a mixed-effects model, where b is a fixed parameter common to all treatments. Graphical analysis suggested that parameter a had a mean value common to all treatments with an added normal random effect due to variability from stand to stand. The parameters that depended on the treatment were Cg ₀ and N ₀, the first as a fixed parameter and the second as a parameter with a normal random component due to stand effects. The resulting statistical equation is the following:

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

(5)

where the indices i, j, and k stand 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. The within-group residual errors ε _ijk are independent normally distributed random variables with mean 0, while the random effects a _ij and N _0ij are such that the vectors (a _ij, N _0ij)^T are independently distributed normal random variables with mean 0 and variance-covariance matrix ψ (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′, cov(a _ij, N _0ij) = ψ, according to Pinheiro and Bates 2000). These hypotheses were subjected to validation.

Furthermore, as a result of this fit, graphical and statistical analysis showed that the data were reasonably consistent with the assumption of aligned inflection points (Ln(Cg ₀), Ln(N ₀)). It could therefore be formally postulated that two parameters, a ₁ and b ₁, exist, satisfying the following equations:

$$ Ln\left({N}_{0_i}\right)={a}_1+{b}_1Ln\left(C{g}_{0_i}\right),i=1,\ldots,n $$

(6)

Finally, Eq. (5) was fitted to the beech data with the constraint of the Eq. (6) system (i.e., by using, in Eq. (5), in place of Cg _0i, the function of N _0i obtained by solving the system of Eq. (6) for Cg _0i). To investigate the possible application of this fitted equation to thinned stands, the data of the thinned stands of initial densities of 4444, 10,000, and 40,000 stems per hectare were plotted together with the trajectory predicted by this equation for stands of 2000 stems per hectare.

Graphical and statistical analyses were performed using R 3.1.1 (R Development Core Team 2012). The nonlinear mixed model was fitted using the nlme package.

3 Results

Starting at various density levels, the observed data for all stands were distributed around trajectories that revealed similar shapes. Equation (5) was satisfactorily fitted to these data (Fig. 1) using a mixed-effects model. Parameter estimates and fit statistics are given in Table 2. The usual diagnostic plots, per stand, of the standardized residuals were performed with the following results. The plots of the standardized residuals vs. fitted values (Fig. 2) showed that the fit was unbiased, with no heteroscedasticity in the within-group errors (this assessment was based on the comparison of these plots to values independently drawn from the standard normal distribution plotted against the fitted values of each stand). The adequate linearity of the QQ plots of the standardized residuals (Fig. 3) supported the normality hypotheses of these error terms. The 1-year lag plots of the standardized residuals did not show evidence of year-to-year correlation patterns of the within-group errors. Moreover, the QQ plots of the random effects (a _ij and N _0ij) also supported the normality hypotheses for these effects.

Table 2 Parameter estimates and statistics of Eq. (5) fit to the beech experimental stand data

Fig. 2

Scatter plot of the standardized within-group residuals vs. the fitted values from the size-density fit of Eq. (5) to the Lyons-la-Forêt beech experimental stand data, giving evidence for the unbiasedness of the fit and the homoscedasticity of the residual errors (each panel is identified by the initial stand density and the stand number separated by a point)

Fig. 3

Probability plots of the standardized residuals of the fit of Eq. (5) to the data of the Lyons-la-Forêt beech experimental stands. The linearity of the relationships provides evidence for the normality of the within-group errors (the ε _ijk in Eq. (5))

Since a linear relationship seemed to exist between transformed initial stand density (Ln(N ₀)) and mean stand girth at breast height (Ln(Cg ₀)) at the onset of mortality (Fig. 1), Eq. (5) was fitted to the stand data with the constraint of the Eq. (6) system. The results appear in Fig. 4 and Table 3, which gives the parameter estimates and statistics of this fit. Since the previous residual analysis did not contradict the normality hypothesis, an analysis of variance could be adequately performed. This ANOVA (Table 4) showed that the hypothesis of alignment of the points of inflection of the trajectories could be accepted at the 1 % significance level.

Fig. 4

Observed size-density data from the varying planting densities—625 to 40,000 trees per hectare—of the Lyons-la-Forêt beech experimental stands and fitted trajectories based on Eq. (5), with the constraints of the Eq. (6) system between the parameters N _0i and Cg _0i for all treatments i

Table 3 Parameter estimates and statistics of Eq. (5) fit to the beech experimental stand data with the constraints of the Eq. (6) system

Table 4 Analysis of variance table for Eq. (5) fits: the equation alone vs. the same equation with the added constraints of the Eq. (6) system

Table 5 gives the RDI and estimated mean girths at breast height (Cg ₀) when mortality first appears in the beech experimental stands and estimated mean girths at breast height (Cg ₁) when the stands reach maximum RDI. Figure 5 shows the thinned stand size-density data together with the trajectory predicted by Eq. (5) with parameters satisfying the constraints of Eq. (6) for stands of 2000 stems per hectare.

Table 5 Stand characteristics estimated by Eq. (5) with the added constraints of the Eq. (6) system when fitted to the beech data: RDI and mean girth at breast height (Cg ₀) when mortality appears, mean girth at breast height (Cg ₁), with corresponding density, when the stands reach maximum RDI

Fig. 5

Thinned size-density data of the Lyons-la-Forêt beech experimental stands and the trajectory predicted by Eq. (5) with parameters satisfying the constraints of Eq. (6) system for stands of 2000 trees per hectare

4 Discussion

The Lyons-la-Forêt density trial with its wide range of densities allowed us to analyze and fit the stand-density trajectories of beech from very early stages of development without any competition-induced mortality until the stage where stands follow the maximum size-density line, with particular attention to the onset of mortality corresponding to the inflection of the size-density trajectories.

4.1 Alignment of the inflection points of the plot density trajectories

The hypothesis of alignment of the inflection points of the stand trajectories was incorporated into Eq. (5) through the constraint of the Eq. (6) system. This hypothesis was not rejected at the 1 % probability significance level, and fitted in this way, Eq. (5) allowed the prediction of a trajectory from initial stand density (the inflection point being determined by the initial density). According to Fisher (1934), it is usual in statistical hypothesis testing to use the 5 % probability level as the standard threshold to accept (or reject) a hypothesis. As compared to the 5 % probability level, the 1 % level increases the probability of accepting a false hypothesis (type II error) but reduces the probability of rejecting a true hypothesis (type I error). This is exactly what we needed because the rather low p value (0.014 in Table 4) is due to treatment 625a, since the inflection points of the other trajectories are relatively straightly aligned. However, there is no reason why the inflection point of the trajectory of treatment 625a should not behave identically to that of the other treatments, apart from the fact that this trajectory is less well-defined than the others because the point where maximum RDI is reached is the farthest from the data range. With the constraint of the alignment of the inflection points, a closer estimation of the initial stand densities (N ₀) could be obtained, particularly for the two highest density treatments (Tables 2 and 3).

4.2 Mortality onset and initial density

The slope of the line locus of the inflection points of the stand density trajectories characterizing mortality onset was greater than the slope of the maximum size-density line. In other words, mortality began to occur at a lower relative density (RDI) for lower initial stand densities (0.31 for stands at 625 trees per hectare compared to 0.58 for stands at 40,000 trees per hectare; see Table 4). This is in line with the results obtained for red pine (Pinus resinosa Ait.) by Smith and Hann (1984) and for loblolly pine (Pinus tæda L.) by VanderSchaaf and Burkhart (2012). However, it is contrary to the results obtained for red alder (Alnus rubra Bong.) by Puettmann et al. (1993), Smith and Hann (1984), and Hibbs and Carlton (1989) and for Douglas fir (Pseudotsuga menziesii [Mirb.] Franco) by Puettmann et al. (1993). For these species, mortality began at a constant RDI value (varying between 0.31 and 0.5 for red alder and close to 0.6 for Douglas fir). It is possible that beech trees with larger crowns to sustain are prone to die when their growing space is reduced since tree photosynthesis would then be insufficient to satisfy their respiration demand (Oliver and Larson 1996). However, the slope of the line locus of the inflection points found for beech may slightly differ when additional data is obtained at Lyons-la-Forêt, especially for the stands planted at the lowest density that have not yet reached the maximum size-density line. The same inconsistency was encountered by Puettmann et al. (1993) and Smith and Hann (1984).

4.3 Maximum size-density line

The constant and slope parameters of the maximum size-density line obtained in this study for the Lyons-la-Forêt experiment (12.997 and −1.388, respectively) slightly differ from the ones previously obtained in a study aimed at establishing the maximum size-density line for beech at a larger scale and based exclusively on data from stands at maximum density (13.68 and −1.566, respectively, for the constant and slope parameters; Le Goff et al. 2011). This may be due not only to the possible variability of the maximum size-density line parameters among stands with differing site conditions but also to the fact that stands at low initial density (625a and 625b) in the Lyons-la-Forêt experiment had not yet reached maximum density, which could give rise to some uncertainty in the parameter estimation. This gives a strong incentive to continue the experiment, almost up to the point when the low-density plots will have reached the maximum size-density line in terms of density.

4.4 Size-density trajectories

4.4.1 Mortality simulation

The two-parameter curve model that proved useful in this study should be adequate—in the absence of an individual mortality model (Monserud et al. 2005)—to model mortality prior to any management treatment in a beech stand simulator, as well as possibly after a pre-commercial thinning, given the result that thinned stands followed the trajectory of unthinned stands of the same number of stems per hectare. This later result is reminiscent of the hypothesis of Pienaar and Turnbull (1973), which was supported by observed data, that “the growth rate of a thinned stand is identical to that of an unthinned stand of the same age and the same basal area as the thinned stand.”

4.4.2 Stand management

The parameters of the size-density trajectories in Table 3 make it possible to obtain the RDI curve at the onset of mortality (RDI_self):

$$ {RDI}_{\mathrm{self}}=\frac{N_{\mathrm{self}}}{N_{\max }}=\frac{\raisebox{1ex}{${e}^{a_1}$}\!\left/ \!\raisebox{-1ex}{$C{g}^{b_1}$}\right.}{\raisebox{1ex}{${e}^a$}\!\left/ \!\raisebox{-1ex}{$C{g}^b$}\right.}=\frac{1}{\left({e}^{0.199}C{g}^{0.244}\right)} $$

(7)

RDI_self is the lower limit of self-thinning—or competition-induced mortality. In stand density management diagrams (SDMDs) commonly used in North America (Drew and Flewelling 1979; Hibbs and DeBell 1994; Newton and Weetman 1994; Penner et al. 2006) and, more recently, in Europe (Pretzsch 2009; Sales Luis and Fonseca 2004), the RDI_self line (in log-log scales) is one of the three management boundary lines positioned, whereas the other two characterize the onset of competition and full site occupancy (VanderShaaf and Burkhart 2012). RDI_self allows the estimation of the loss of growing stock when the stand exceeds this upper limit of density and then provides useful information as when to carry out thinnings during stand rotation. Since, for beech, the self-thinning line does not appear to be parallel to the maximum size-density line and has a steeper slope, therefore, the initial stand density must be considered in order to predict the RDI level at which competition-induced mortality will start (the higher the initial stand density is, the higher the RDI at which mortality begins).

SDMD diagrams could be developed for beech with additional developments to characterize the particular RDI lines previously mentioned and not already established and to visualize alternative silvicultural regimes for beech in terms of RDI trajectories with quantification of the growth and yield of stands using available growth models for beech (Le Moguedec and Dhôte 2012).

4.4.3 Application to other species

The size-density trajectory model constructed in this study could be transposed and applied to other species. This could be the case in France for the following well-represented species for which long-term data are available, including sessile and pedunculate oak (Quercus petraea and Quercus robur), maritime pine (Pinus pinaster), Douglas fir (P. menziesii), and Corsican pine (Pinus nigra subsp. Laricio). A national network, the “Cooperative for data on forest tree and stand growth,” which experiments with various thinning regimes (including “no thinning”), exists for these species (Bédéneau et al. 2001).

5 Conclusion

This research contributes to the analysis of size-density trajectories, using a unique experiment to test the effects of a wide range of planting densities on beech stand development over a long period of time. Our results show that the inflection points of the trajectories are aligned and that mortality begins at decreasing RDI values with decreasing initial densities. The maximum predicted size-density line is close to the one established earlier at a larger scale (Le Goff et al. 2011), giving a wider applicability to the size-density trajectory equations fitted in this experiment. As a result, two useful lines for the management of even-aged beech stands could be drawn in Ln(N)-Ln(Cg) diagrams (SDMDs): one indicating conditions in which mortality begins and one indicating conditions of maximum stand density. Incorporating size-density trajectory equations in stand growth simulators would also improve the prediction of tree mortality. Finally, the segmented model fitted here to the beech density trajectories could be easily adapted to other species if appropriate data are available, giving rise to the comparison of trajectories between species and, more specifically, to the onset of mortality, which reflects the ability of each species to survive strong inter-tree competition.