Self-thinning dynamics in cork oak woodlands: providing a baseline for managing density

Aim of study : The study aims to evaluate the maximum potential stocking level in cork oak ( Quercus suber L.) woodlands, using the ecologically-based size-density relationship of the self-thinning law. Area of study : The study area refers to cork oak forests in mainland Portugal, distributed along its 18 districts from north to south. Material and Methods : A dataset with a total of 2181 observations regarding pure cork oak stands was collected from the Portuguese Forest Inventory (NFI) databases and from research plots. The dataset was subjected to two filtering procedures, one more restrictive than the other, to select the stands presenting the higher stocking values. The two resulting subsets, with 116 and 36 observations, from 16 and 10 districts of mainland Portugal, respectively, were then used to assess and describe the allometric relationship between tree number and their mean diameter. Main results : The allometric relationship was analysed and modelled using the log transformed variables. A slightly curvilinear trend was identified. Thus, a straight line and a curve were both fitted for comparison purposes. Goodness-of-fit statistics point out for a good performance when the data is set to the uppermost observed stocking values. A self-thinning line for cork oak was projected from the estimated relationship. Research highlights : The self-thinning model can be used as an ecological approach to develop density guidelines for oak woodlands in a scenario of increasing cork demands. The results indicate that the recommendations being applied in Portugal are far below the maximal potential stocking values for the species. It is therefore of the utmost importance to review the traditional silvicultural guidelines and endorse new ones.


Introduction
Cork oak (Quercus suber L.) is a species of great ecological, economic and social value with a natural distribution limited to the western Mediterranean basin. It occurs in north Algeria, Morocco and Tunisia in Africa, and in France, Italy, Portugal and Spain in Europe. Quercus suber is the major softwood species in Portugal, occupying around 737 000 ha (ICNF, 2013), which represents 23% of the forested area of the mainland. According to the Portuguese Cork Association, Portugal is the leading producer of cork, with 49% (around 100,000 tons) of the world production. Cork oak forests are responsible for 14 M tons of CO 2 fixed per year; the environmental services of these forests are credited to be around 100 € per hectare per year (APCOR, 2016). The whole system sustainability is currently severely threatened by plagues and diseases, such as the insect Coroebus undatus or the fungi-like Phytophthora spp. Weather conditions, with the drought being a limiting factor, and inadequate practices may contribute to increase its vulnerability to biotic agents. Ongoing practices that involve the intensification of soil disturbances (soil tilling), to control shrub invasion and/or to promote pastureland, have also a negative impact on the success of tree natural regeneration (Simões et al., 2016). Additionally, other warning issues concerning the sustainability of these ecosystems have been noticed. Climate change scenarios forecast rising temperatures and decreasing precipitation in the dry period, with consequences in the environmental stress that promote the propagation of plagues and diseases. These factors are expected to be associated with a modification of the potential area of distribution of the species in the inland (Pereira et al., 2006). In parallel, a decrease in the crown cover and in the stocking values of the stands was reported (Paulo et al., 2015), based on data from the Portuguese National Forest Inventories. Meanwhile, the cork oak demands from industry are increasing. The reduction of stand density and the increase in demand for cork might lead to an overexploitation of the existing trees. As stated in Marques et al. (2000), referring to cork oak woodlands in the North part of Portugal, forests with lower values of density were associated with an overexploitation of the trees, namely at the crown level. According to the authors, the solution encompasses increasing the density of the woodlands artificially or through the promotion of natural regeneration. Ribeiro (2015) also state that the sustainability and the resilience of the cork oak forest are firmly founded in the structure and density of trees for maintaining a continuous cover. Solutions for a sustainable management of the species cannot therefore be alienated from a comprehensive knowledge of the optimum values of stand density. This led to the research question of whether the traditional stand density guidelines correspond to "an optimum" or instead, are placing at risk the sustainability of the cork oak system. Natividade (1950) and Montero & Cañellas (2003) point out that the density of the cork oak stands is a sensitive issue, being quite difficult to reach conclusions regarding "optimal" or acceptable densities. Contrarily to other species, the particular features of cork oak and of its management, often as silvopastoral systems, make it difficult to follow spacing guidelines. Furthermore, the stands originated by spontaneous regeneration are quite heterogeneous with a non-uniform spatial distribution. Typically, the recommended values of density are empirically based. Table S1 [supplementary] presents the number of trees per hectare for given values of circumference over bark at 1.30 m height level (C) and for the corresponding diameter values (d = C/π), according to the bibliography. Spacing among trees such as 4m × 3m, 8m × 2m or 4m × 4m (625-825 trees. ha -1 ) have been proposed for use in the plantation of new cork oak stands (Costa & Pereira, 2007). The stands are managed with prescription of heavy thinning in young stages to reduce the density in one half. The thinning is performed on stands with less than 20 years, i.e. before the trees reach a minimum size for debarking. The density after thinning reaches the 300-400 trees. ha -1 , approaching the values pointed out by Natividade (1950) for a girth over bark interval of 60-70 cm ( Table  S1 [supplementary]).
The optimal densities for a given forest system can be assessed from the allometric relationship between density (Y) and size (X) of the trees, expressed as a power law (Y = aX r ), with a and r being allometric coefficients from the ecologically-based concept of self-thinning law. The concept considers that, in the evolution of monospecific even-aged populations of plants experiencing complete crown closure, as trees develop, average space available for growth decreases, and the number of living trees decreases. According to the concept, mortality in pure stands is expected to occur due to the increasing intraspecific competition effect. Although density and size can both be represented by different metrics, the preferable variables are the ones originally considered by Reineke (1933): number of trees (N, trees.ha -1 ) for density, and mean tree diameter (dg, cm) for size. In this context, a biological interpretation of the allometric coefficient r arises. In pure stands, r indicates the intensity of intraspecific competition (self-tolerance), defined as the proportion of trees, ∂N/N, eliminated by a certain increase in average diameter, ∂dg/dg. Because the intensity decreases with species tolerance, r varies among tree species, usually between 1.4 and 1.9 (Zeide, 1985). Tolerant species tend to show a higher intercept (Harper, 1977) and a lower slope (Zeide, 1987) due to the variation in self-tolerance. The effect of self-tolerance has been mentioned to influence the shape of the trajectory, with some authors suggesting it to be a curve or a set of curves instead of assuming a straight line (e.g. Tang et al., 1994;Río et al., 2001 andCharru et al., 2012). Theoretically, it is expected that a deviance from the straight line should occur at both ends, with a concave shape: in young stands, before the onset of competition begins, and later on, in older stands, due to the decrease in self-tolerance (Zeide, 1987). According to Zeide findings, the decrease in self-tolerance is responsible for the fact that the canopy gaps cannot be filled as promptly as they occur (Zeide, 2005). Moreover, in older stands, insects, diseases, and other disturbances may reduce stocking more rapidly than the residual stand is able to re-occupy growing space (Shaw, 2006).
The maximum size-density relationship, also known as self-thinning line, supports the definition of the density metric, Stand Density Index (SDI) by Reineke. The SDI has demonstrated to be an interesting methodology to evaluate density and to provide an ecological baseline for forest management. Specific points in stand development, other than the maximum size-density relationship, such as crown closure or imminent competition-mortality, can be defined in terms of tree size and stand density, through the use of SDI.
In Portugal, Luis & Fonseca (2004) developed a Stand Density Management Diagram (SDMD), with isolines for SDI and volume variables, to help in the management of pure stand of Pinus pinaster. The structure of the SDMD allows the users to assess the levels of occupancy of the stands, measured by SDI, in the stocking diagram, and to decide for intervention (thinning) based on the levels of SDI and on the management goals. The adequacy of the methodology has been validated for other forest systems in other regions. Cruz & Bascuñan (2014) developed a silvicultural management scheme for Nothofagus dombey to be used as a density guide, considering a combination of the SDI with the potential available area and the variation in diameter growth. The SDI was used to identify the minimum density and the onset density for tree crown competition, defining a range of adequate density values to optimize the diameter growth of the trees with minimum interventions. Aigbe & Omokhua (2015) adopted a similar methodology based on the self-thinning law to develop a stocking chart and a stand density management model for tropical rainforests.
The present study aims, specifically, to provide results about the limiting number of trees that can stand alive in the cork oak woodlands, per unit of area, based on the concept of self-thinning law, as introduced by Reineke (1933). The authors hypothesize that the occupancy level achieved with the traditional density guidelines remain below the suitable reference values for the optimum growth-density interval of the species. In case this hypothesis holds true, silvicultural guidelines for the species should be revised.

Material
The data used in this study includes information collected during two assessments of the Portuguese Forest Inventory (NFI) -NFI4 and NFI5, carried out in 1995-1996 and in 2005-2006, respectively. Additional data from research plots (RP) was provided by N. Ribeiro (University of Évora, Portugal). These RP are located in Machuqueira do Grou, in Ribatejo region. A synthesis of the RP data can be found in Aronson et al. (2009), pp. 57-58, as an example of a site profile corresponding to the privately owned formations common in southern Portugal. The three data sources refer to pure stands of cork oak. The plots spread over the country (Fig. 1) throughout the 18 districts of mainland Portugal. Table  1 summarizes the characteristics of the sampling data.

The allometric model and the Stand Density Index
The allometric relationship, which is ecologically based, assumes that the variation in the number of trees is proportional to their growth in diameter Integrating both sides of Equation (1) produces the general expression for the allometric model, given by where "k" and "r" are allometric constants. As the tree diameter presents a high correlation with tree crown size, the metric dg can be assumed to be a surrogate of crown cover. Applying a log transformation (ln) to Equation (2), a linear limiting relationship known as "self-thinning" law by Reineke (Reineke, 1933) is defined as Specific points in stand development other than the maximum size-density relationship, such as crown closure or imminent competition-mortality, can be defined in terms of tree size and stand density, through the use of the Stand Density Index (Equation (4)).
To account for curvature in the self-thinning trajectory, the bend can be estimated as a space plan by generalizing Equation (3) to a polynomial of second degree (Equation (5)).
where ¯ ln dg represents the mean value of the transformed variable ln dg.
The self-thinning line can be considered, therefore, as an upper bound corresponding to the maximum density achievable by the species.
For any stand, relative density can be described directly by the self-thinning line, where the maximum size-density line represents 100% of density and other size-density combinations are expressed as a percentage of the maximum, or can be based on the calculation of the SDI metric.

Estimation of the self-thinning line
The exact determination of the self-thinning trajectory for any population presents difficulties related with the composition of the stand -pure vs mixed (Río et al., 2016) -and, for a given species, with the source of data and the estimation method. For a given dataset obtained from a forest inventory, only part of the sampled stands are expected to be in a true self-thinning mode (Shaw, 2006). In such cases, it is necessary to apply methods that are insensitive to observations under stocked conditions or to filter the data for upper boundary plot selection (e.g. Luis & Fonseca, 2004).
To filter the given dataset of 2181 observations (Table 1), two procedures were implemented involving analytical algorithms. In the first procedure, for plots with the same value of N, the algorithm selected the one with the largest value of dg, and, for plots with an equal value of dg, it selected the one with the highest N value. These filtering criteria eliminated observations in the three datasets, along the range of cork oak geographical distribution, retaining 116 observations (plots), hereafter designated as Sample 1. It contains 50 observations from NFI4, 63 from NFI5 and 3 from the RP dataset. The distribution of Sample 1 plots includes 16 out of the 18 Portuguese districts.
A detailed analysis of Sample 1 detected the existence of plots with similar values of N that presented a large variability in the dg values and plots with close values of dg presenting great differences on N. Thus, it was decided to apply additional filtering criteria: for observations with a difference in N equal to 1 tree.ha -1 , the one with the highest value of dg was retained; secondly, for observations with a difference in the value of dg up to 1 cm, the plot with the highest value of N was selected. This procedure resulted in a second subset -Sample 2with 36 observations, being 17 from the NFI4 dataset and 19 from NFI5. The distribution of the plots retained in Sample 2 includes 10 out of the 18 Portuguese districts, as shown in Fig. 1.
The two subsets of data, Sample 1 and Sample 2, are characterized in Table 2.
The straight linear model, Equation (3), and the quadratic model, Equation (5), were fitted for each data sample, using least square regression procedures, and residual analysis was conducted to evaluate the adequacy of the model. The variance inflation factor (Myers, 1990) was used to measure linear dependencies among regressors in Equation (5). Goodness-of-fit was based on the statistical criteria coefficient of determination adjusted (R 2 adj) and root mean square error (RMSE).
Both models, Equations (3) and (5), describe an average limiting relationship between plant density and average plant size. As stated by Shaw (2006), it is expected that the sample data includes a part of stands developing at the maximum growing stock and other being less than optimally stocked, namely at young ages. Hence, after fitting, the estimated lines should be moved up vertically to be set above the observed density-size pairs, in order to represent a potential maximum limiting line. The methodology selected to fix the upper limit was the one proposed by Luis & Fonseca (2004). It consists on raising the intercept term to the upper limit of the confidence interval at the 95% one tail level [α = 0.05; t (1 -α, n -p df), where p represents the number of the estimated parameters]. The statistical analyses were conducted with the JMP® (v. 10.0) software of SAS® Institute Inc.. Fig. 2 shows the distribution pattern of density-size for both sets of data used in the study, and the average fitted trajectory obtained with the considered models -Equations (3) and (5).

Results
The parameters (and the standard errors) of the estimated models for defining the self-thinning trajectory are shown in Table 3. Analysis of the residuals for the fitted models did not present evidence of violation of regression assumptions, although a high dispersion of residuals was found for Sample 1 (Fig.  3). The curvilinear models did not present problems of multicollinearity, as measured by the variance inflation NFI4 and NFI5 refer to data from the National Forest Inventories, RP refers to research plots located in an oakland forest in the central part of Portugal, and Obs. indicates the number of plots. The variable N represent the number of trees per hectare (trees.ha -1 ) and dg represents the quadratic mean diameter over bark (cm) measured at the height level of 1.30 m. The symbol sd refers to the statistic standard deviation.  All the estimated coefficients are significant at the p = 0.05 level. The goodness-of-fit statistics R 2 adj and RMSE for the proposed models are also shown in Table 3.
Considering the results obtained for Samples 1 and 2, the estimates for the self-thinning trajectories were based on data from Sample 2.
The intercept for the functional forms corresponding to Equations (3) and (5), were raised to the upper limit of the confidence interval at the 95% one tail level [α = 0.05; t = 1.17]. The estimates are, respectively: lnN = 12.302 -1,806 ln dg (6) ln N = 13.096 -1.964 ln dg -0.590 (ln dg -3.661) 2 (7) A graphical comparison of the density values resulting from the proposed self-thinning trajectories with the silvicultural guideline values presented in Table S1 [supplementary], according to the values by Natividade (1950), Lamey (1893), Montoya (1985Montoya ( , 1988, and Montero & Cañellas (1999, 2003, is provided in Fig. 4. Fig. 4 also presents the density lines corresponding to three values of Stand Density Index (expressed in percentage): SDI = 35%, SDI = 45% and SDI = 60%. SDI = 35% corresponds to the lower limit of full occupancy, SDI = 45% is in the range of a full occupancy, and SDI = 60% corresponds to the density level linked to a probable start of mortality by competition. These three lines were drawn from the self-thinning model expressed in Equation (6).

Adequacy of the methodology and represen tativeness of the self-thinning trajectory
The methodology used, based on the self-thinning concept and on the stand density index by Reineke, is recognized as adequate and important to help evaluate Table 3. Coefficients and fit statistics of the estimated models to support the definition of the self-thinning line for cork oak.

Model (Equation)
Subset  stand density in different species forest systems. The research developed in this field and the documented case studies that are found in the literature include pure (even and uneven aged) and mixed stands. To the best of the authors' knowledge, there are no reasons to doubt the suitability of the methodology for the particular case study, involving cork oak forests. The algorithms used in the sampling procedure retained the 36 highest values of density-dimension (N, dg) pairs, which constitute a representative group to define the self-thinning line for the species (Sample 2). The estimated models for the self-thinning trajectories (Equations (6) and (7)) explain the high proportion of variability found in the relative variation of the number of trees per unit of area (R 2 adj > 0.8). The fact that the quadratic term was found to be statistically significant (Equation (7)) is in accordance with results provided in the literature for other species and is moreover theoretically expectable. The curvature from the straight line is likely to occur, as explained by Zeide (1987), both in young stands, before the onset of competition, and in older stands, when the canopy gaps can no longer be filled sufficiently fast, due to a decline in self-tolerance. The maximum density values obtained directly by the straight line model (Equation (6)) and the curvilinear model (Equation (7)), as shown in Fig. 4, indicate the possibility of promoting the development of the cork oak woodlands with densities greater than the usual prescribed values. The selfthinning trajectory defined by the curvilinear model (Equation (7)) describes a density-size relationship that is beyond the ones presented by previous studies. Although there is no statistical evidence (at p < 0.001) to refute the validity of the quadratic self-thinning relationships, the authors consider that the straight linear model provides a good approximation to the self-thinning trajectory for cork oak forests. The gain in terms of percentage of explained variation provided by Equation (7), comparatively to the variation explained by Equation (6), is 6.4%, as measured by the difference in the values of the R 2 adj statistic. In a study performed on a set of different species, Charru et al. (2012) concluded that the curvilinear model was significantly more accurate for 7 temperate and Mediterranean tree species. Nevertheless, for the three oak species included in the study (pedunculate, pubescent, and sessile oaks), the quadratic term was not statistically significant, which led the authors to retain the linear relationship pattern for those oak species.
The self-thinning model expressed by Equation (7) might be considered, and it will be so in the present study, as an upper bound corresponding to the maximum possible density for the species. The estimated allometric coefficient r found by Charru et al. (2012) for oak species was in the range of -1.911 to -1.758. For cork oak, in this study, the allometric coefficient was estimated as -1.806, which is remarkably similar to the value presented in Charru et al. (2012) for the pubescent oak (-1.809) and different from the fixed value of r = -1.605, suggested by Reineke (1933).  Natividade (1950), Montoya (1985Montoya ( , 1988, Montero & Cañellas (1999, 2003, and Lamey (1893) (Table S1 [supplementary]), to the straight linear model (Equation (6)) and to the quadratic model (Equation (7)). Hypothesized density values of SDI, calculated from Equation (6) are also shown. Density is represented by number of trees (N, trees.ha-1) and size refers to circumference over bark at 1.30 m height (C, cm).

Optimal density values
Assuming the self-thinning line as defined by the allometric relation stated in Equation (6), it is possible to define density lines corresponding to percent values of maximal occupancy. Based on Long (1985) metrics for relative density values of SDI, an SDI of 25% has been associated with the transition from open grown to populations at the beginning of competition, an SDI close to 35% defines the lower limit of full occupancy, and SDI values of 55-60% have been associated with the lower limit of self-thinning. The selection of the density value will depend on the objectives of the landowners and on the type of system to promote: agrosilvopastoral or mainly forest. It is noteworthy that the density values referred by Montoya (1985Montoya ( , 1988 and the indicators corresponding to the lower bound defined by Montero & Cañellas (1999, 2003 are in the region between the lower limit of full occupancy of the stands and the onset of self-thinning, approaching the relative density of 45%. Hence, the figure of SDI equal to 45%, obtained from Equation (6), might be considered a reasonable option for silvopastoral systems. Higher values of density for natural cork oak ecosystems are expected to occur. Furthermore, regarding forest management mainly for cork, a value of SDI close to 60% is an interesting option to consider, without any evident reason for concern. The trajectory defined by Equation (6) fits very closely to the values reported by Lamey (1893) for the Esterel region in France. Those values correspond to the higher known densities recommended for the species. Natividade (1950) criticized the values presented by Lamey (1893), referring that for those densities the crowns could cross or overlap and suggesting that the values should be reduced by 10 to 20%. This recommendation was empirically based. In denser stands the trees present narrower angles of insertion of the branches. In fact, the shape and the size of the crown can be affected not only by the density but also by silvicultural practices such as pruning. Differences in the architecture of the trees from the South to the North regions of the country were also reported in Fonseca & Parresol (2001), which were due to the type (and density) of the woodlands but also to the silvicultural practices associated with managing them as silvopastoral systems or as forests. Furthermore, the effect of the topography, namely of the slope, might be of importance. For a given density, measured by the number of trees per unit of area, a stand in a hilly terrain will possess a greater number of trees than a stand in a plain area. This is due to the standardization of the definition of area, as it stands for the horizontal projection.

Aspects to be addressed
The bark of cork oak is a highly valued non-wood forest product (NWFP) with major importance in the Portuguese forestry sector. Regarding potential effects of tree density on the thickness of cork, Ribeiro (2015) refers that the cork oak stands can be managed with stand basal areas greater than 12 m 2 ha -1 , without compromising the productivity or the thickness of the cork. A basal area of 12 m 2 ha -1 corresponds to the traditional size-density guideline proposed for the species in Portugal (Natividade, 1950;DGRF, 2006), which, expressed as relative density (percent value of SDI), is associated with the lower limit of full site occupancy. The optimal density for cork yield without greatly affecting the thickness of the cork can be extended, based in the results, to relative densities around 55% of SDI. Although it is beyond the scope of this study, the authors hypothesize that increasing stand density might affect the growth of phellogen, which would subsequently require stripping cycles longer than the typical 9-12 years. This might not, however, interfere significantly in the total productivity, as the stripping surface would increase at stand level, with a greater number of trees involved in the production.
In the context of silvicultural guidelines, it is important to attend to the ecological dynamics of the species at earlier ages. The seedlings growing under the canopy of the mature trees are tolerant to the shadow but present a fair development (Montero & Cañellas, 2003). Hence, a conflict between the promotion of very high densities and the promotion of the natural regeneration might happen if the management favours to retain the forest with densities near the self-thinning line. The systems that promote grazing are not favourable to the success of the natural regeneration. Measures to help natural regeneration are required in those systems.
The research developed concerns the study of selfthinning dynamics in cork oak woodlands, to assess the maximum density attainable by the species and to provide a scientific baseline for managing the density of the species in Portugal. The implementation of the results in new management guidelines requires awareness of landowners, forest technicians, stakeholders and of the Minister of Agriculture, Forestry and Rural Development in order to promote discussion and review the policies that have been followed for the species. It is noteworthy that, according to Portuguese legislation (DR, 2015), the density of reference to support forestry activities such as plantation, seedling and promotion of natural regeneration, is just 450 trees per hectare. Evaluation of actual stocking levels of the forests, using the ecologically-based self-thinning concept, will noticeably confirm that new guidelines are needed.

Conclusions
The noticed reduction in the forest crown cover during the last decades and the expected decrease in forest area, driven by climate change, will severely impact the economic, social and ecological benefits of cork oak lands. The present work focused on the study of optimal density with an ecologically based approach. The estimation of the self-thinning trajectory found for the species in this study is supported by theoretically accepted ecological assumptions. A curvilinear trend was noticed, which agrees with the self-tolerance dynamics of the trees.
Managing higher densities is desirable, in response to the increasing cork industry demands, and also attainable, as indicated by the results of this study.
The comparison of relative densities, based on the stand density index by Reineke, with traditional silvicultural proposals, has shown that the guidelines must be rethought. Classic orientations have led to woodlands kept with densities considerably reduced, far below the competition levels associated with mortality induced by competition. This can be explained by excessive caution when stating density regulation guidelines for the species in Portugal, and also by a generalization of the recommendations provided by Natividade (1950) for all the cork oak systems, independently of promoting or not multiple uses. The proposed optimal densities are closer to the values mentioned for the French forest in the Esterel region, by Lamey (1893), in the late 19th century. This is a strong argument in favour of the credibility of the study's results. It was shown that it is possible to manage the species with densities higher than those traditionally recommended and followed for decades, without compromising the natural equilibrium of these ecosystems. Promoting stands with high densities might also be a valuable adaptive action in case the potential area for the species is reduced.