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 ₂ 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 /p), 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 and Charru 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

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)

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 factors ( VIF ): the VIF was 1.016 for Sample 1 and 1.116 for Sample 2.

All the estimated coefficients are significant at the p 0.05 level. The goodness-of-fit statistics R ² 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 [a = 0.05; t = 1.17]. The estimates are, respectively:

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.

T. Fonseca wishes to express her gratitude to Prof. Carlos Pacheco Marques from the University of Trás-os-Montes e Alto Douro, for his commitment in promoting the research of cork oak stands in the Northern region of Portugal, namely through the Project PAMAF 4007.

The authors express their appreciation to ICNF for providing the data of the National Forest Inventories and to Prof. Nuno Ribeiro from the University of Évora for sharing the data collected in the research plots.