Optimizing the debarking and cutting schedule of cork oak stands

Abstract

Key message

Optimal management of cork oak forest stands was analyzed for different site indices and cork growth rates. Optimal debarking intervals varied during the rotation and were sometimes shorter or longer than the officially recommended range of 9–14 years.

Context

Quercus suber L. is one of the most important multipurpose tree species in the Mediterranean area. Its main product is cork, appreciated for its elasticity, impermeability, and thermal insulation properties. Cork oaks are debarked at constant intervals, which vary from 9 to 14 years depending on the area. However, since the growth rate of cork is not constant during the rotation, it may be optimal to use variable debarking intervals.

Aims

This study optimized the debarking and cutting schedules of Quercus suber stands and analyzed the influence of economic and stand-related factors on optimal management.

Methods

The study employed a simulation system where the existing growth and yield models for Quercus suber were used with a non-linear derivative-free optimization algorithm. Discount rates and cork prices were tested as economic factors and cork growth rate and site productivity as stand-related factors.

Results

The optimal debarking interval varied during the rotation. Increasing cork growth rate increased the optimal number of debarkings and shortened their interval. Decreasing discount rate increased the optimal number of debarkings during rotation while decreasing cork price decreased the number of debarkings.

Conclusion

The profitability of the management of cork oak stands depends on site fertility and stand density; management is not profitable on poor sites or at high discount rates. This study is the first that simultaneously optimizes the cutting and debarking schedule of cork oak stands, allowing the debarking interval to vary.

1 Introduction

The importance of non-wood forest products (NWFP) has gained attention in the past years within the scientific community (Calama et al. 2010; Wolfslehner et al. 2013) and elsewhere because of the importance of NWFP for rural development. Management for sustainable production of NWFPs is especially relevant for Mediterranean forests, which are characterized by their multifunctionality and the diversity of products they provide (Bravo et al. 2011). Among these NWFPs, cork is one of the most representative for Mediterranean forests, providing annual revenues of 325 million of euros at the European level, Portugal and Spain having the largest areas of cork oak forest with 713,000 ha in Portugal and 475,000 ha in Spain (Pereira 2007). Cork is the outer bark of cork oak. It has good properties in terms of elasticity, impermeability, and thermal insulation (Pereira 2007). Cork extraction is a sustainable process because cautious cork extraction does not damage the tree, and following the extraction, new bark grows. The properties of cork make it appropriate for a variety of applications. Moreover, it is natural, renewable, and recyclable material. Most of the cork is used as bottle stoppers. However, cork is also used for flooring, paneling, heat and sound insulation, decorative purposes, transportation industry, and sorption applications, among other things (Duarte and Bordado 2015).

Two types of Quercus suber stands are distinguished in Spain by the ecological characteristics (climate and soils) and management of the stand. One is the “dehesa,” which is an open woodland silvipastoral system (60–80 trees ha⁻¹), where cork oaks provide acorns and shelter for grazing animals (Joffre et al. 1988). Dehesa corresponds to the Portuguese “montado.” Revenues are derived not only from cork production but also from cattle. The other types of stands are called “cork oak forests” (280–400 trees ha⁻¹) (Pasalodos-Tato et al. 2018). They are pure forest systems and common for instance in Catalonia (northeastern Spain) and the provinces of Cádiz and Málaga (southern Spain).

Current recommendations for the management of cork oak forest propose regeneration by the uniform shelter wood method when the trees reach 150–180 years’ age and up to five thinnings during the rotation (Montero et al. 2015). Formation pruning is recommended at 8, 14, 27, and 36 years, and two pre-commercial thinnings are also recommended (Ovando et al. 2010). Cork extraction starts when around 50% of the stems in the stand reach 70 cm circumference over bark. The cork obtained from the first and second debarking (called virgin and secondary cork, respectively) is of low quality, suitable only for disintegration industry. The following debarking operations take place every 9–14 years, depending on climatic and soil characteristics (Montero and Cañellas 2003; Pereira and Tomé 2004). These later debarkings yield high-quality cork suitable for wine stoppers. The minimum debarking interval accepted by official recommendations is usually 9 years, but shorter intervals may be authorized in some Spanish regions.

Cork growth presents some peculiarities, which may put the traditionally accepted regular debarking schedules into question. Increment in cork thickness is not closely correlated with the growth rate of trees (Sánchez-González et al. 2007a) and the amount of cork obtained in a debarking operation depends on the number of times the tree has been debarked earlier (Machado 1944; Tomé et al. 1999). These factors may lead to the conclusion that uneven debarking intervals may increase cork yield. For example, it may be optimal to use shorter debarking intervals in the middle of the rotation when the growth rate of cork is fast. Moreover, economic conditions such as discount rate and cork price may also influence the timing of debarking and silvicultural operations of cork oak stands. These factors have not been included in the optimization of debarking schedules at the stand level.

Optimizing the management of cork oak stand is a complicated problem since the profitability of management depends not only on debarkings, but also on the thinning treatments and rotation length of the stand. Cork oak wood may not be economically valuable, but thinning treatments and stand density affect the growth of trees, which in turn affects the growth of the circumference of tree stems and therefore also cork yield. Thinning treatments also determine how the cork production is distributed among different trees in the stand. Moreover, as the growth rate of cork thickness varies during the rotation, optimal debarking intervals are not necessarily equal.

Most of the earlier studies on optimal cork oak forest management are focused on the landscape level (Borges et al. 1997; Falcão and Borges 2005; Costa et al. 2010; Borges et al. 2014; Palma et al. 2015). Falcão and Borges (2005) and Costa et al. (2010) used numerical approaches to optimize the management of cork oak forests subject to even-flow constraints. Palma et al. (2015) integrated the effects of climate change on the optimal management of cork oak and holm oak systems. They simulated different treatment schedules for cost oak stands and found their optimal forest-level combination using mixed integer programming. However, since stand management was constrained by forest-wide considerations, it is likely that the management schedules in the solution may not coincide with optimal stand-level prescriptions. The situation is the same in several other studies that use combinatorial optimization to propose stand-level prescriptions in forest level management planning context (e.g., Borges et al. 2014). Since the potential number of combinations of debarking years, thinning years, thinning intensities, and rotation length is huge, stand-level optimization may provide information useful to management planning. It may provide information that helps to reduce the number of alternative prescriptions to be considered in forest level management planning. At the stand level, Paulo and Tomé (2017) used simulation to compare debarking intervals of 9–14 years but debarking intervals and cutting schedules were not optimized.

Methods available for stand management optimization include dynamic programming, non-linear programming methods (for instance Hooke and Jeeves), and population-based methods (Hyytiäinen 2003; Pukkala 2009; Pasalodos-Tato 2010). The choice of the method depends on the level of detail in the simulation of stand development so that non-linear optimization is often used with individual tree growth models. For example, Brodie et al. (1978) used dynamic programming with stand-level growth models and de Miguel et al. (2014) employed non-linear programming with distance-independent individual tree growth models. Since individual tree models allow the most detailed simulation of stand development, they were used in this study, together with non-linear programming. There are several alternatives for non-linear programming. This study employed the method of Hooke and Jeeves (1961), which has been used widely in forestry (see Hyytiäinen (2003) and Pasalodos-Tato (2010) for references). It has also been applied in studies dealing with the joint production of timber and different NWFPs (Palahí et al. 2009; Miina et al. 2010; de Miguel et al. 2014; Pasalodos-Tato et al. 2016). The method has performed well when alternative optimization algorithms have been compared (Pukkala 2009, Arias-Rodil et al. 2015).

No studies have been conducted that simultaneously optimize the debarking and cutting schedules of cork oak stands, allowing both debarking and cutting interval to vary. Therefore, the objective of the present study was to optimize the management of cork oak stands, considering the debarking schedule, thinning treatments, and rotation length. Existing growth and yield models for cork oak (ALCORNOQUE, Sánchez-González et al. 2007b) were used in a simulation software. To have a sufficient set of models for simulating the dynamics of cork oak stands, a self-thinning model was fitted based on data from National Forest Inventory plots. The models on the dynamics of cork production and tree growth and survival were linked with an optimization algorithm, which was used to find optimal management schedules for a set of cork oak stands under different discount rates and cork prices.

2 Material and methods

2.1 Study area and analyzed stands

The study area is the natural park Los Alcornocales in South Spain. It is one of the most important cork production areas in Spain and can be considered representative of Spanish cork oak forests. The area is located in the southern part of Andalusia, and it is characterized by a mild Mediterranean climate, the mean annual temperature being 16–18 °C and the annual precipitation around 1200 mm. Most rains occur between autumn and spring with a dry period in summer.

The analyzed stands differed in terms of cork growth rate (later referred to as cork index, which is defined as cork thickness 9 years after debarking) and tree growth rate (site index, which is dominant height of trees at 80 years). The analyzed stands cover the existing range of variation in the growth rate of cork and trees in the Los Alcornocales region.

Three stands (stands 1, 2, and 3 in Table 1) that differ in cork index but have the same site index were analyzed to test the influence of cork growth rate on the optimal management of Quercus suber stands. The cork indices were 20.47, 29.52, and 38.67 mm at 9 years since debarking, representing cork index classes of 20, 30, and 40 mm, respectively. Secondly, two stands with nearly the same cork index class (30 mm) but different site indices (stands 2 and 4 in Table 1, with site indices 14 and 8 m, respectively) were used to analyze the influence of site index on optimal management. Only the cork produced by stems was included in the analyses.

Table 1 Summary of the characteristics of the analyzed stands

2.2 Models for simulating tree growth

Stand development was simulated on individual tree basis. Site productivity was described by site index, which is equal to dominant height at 80 years. The following model was used to calculate the site index and to predict the height growth of dominant trees (Sánchez-González et al. 2005):

$$ H=\frac{20.7216}{1-\left(1-\frac{20.7216}{SI}\right)\bullet {\left(\frac{100}{T}\right)}^{1.4486}} $$

(1)

where H is dominant height (m) and T is stand age (years).

Annual diameter increment was predicted from the following (Sánchez-González et al. 2006):

$$ idu=0.18+0.79\frac{1}{N}-1.02\frac{1}{SI}+2.45\frac{1}{du} $$

(2)

where idu is the annual diameter increment under cork (cm), N is the number of trees per hectare, SI is the site index (m), and du is the diameter at breast height under cork (cm).

Tree height was calculated with the model developed by Sánchez-González et al. (2007c):

$$ h=1.3+\left(H-1.3\right){\left(\frac{du}{D}\right)}^{0.4898} $$

(3)

where h is the total tree height (m), H is the dominant height (m), du is the diameter at breast height under cork (cm), and D is the dominant diameter under cork (cm).

Based on data from the second National Forest Inventory (NFI2), a self-thinning model was developed using data from plots measured in cork oak forests. The plots that were assumed to be at the self-thinning limit were chosen, and a model that related the maximum stand density to the mean diameter was fitted:

$$ {N}_{max}=\exp \left(12.581-1.8987\bullet \ln \left({d}_q\right)\right) $$

(4)

where N_max is the maximum possible number of trees per hectare, and d_q is quadratic mean diameter of trees over cork (cm). The fitted self-thinning model (Eq. 4) was modified following the method proposed by Vanclay and Sands (2009) resulting in a model that predicts some mortality already below the self-thinning limit (Guzmán et al. 2012):

$$ {N}_{max}=\frac{N_0}{{\left(1+{\left(\left(\frac{\pi {N}_0}{40000}\right)\frac{d_q^2}{37}\right)}^{1.3}\right)}^{\frac{1}{1.3}}} $$

(5)

where N₀ is the initial stocking (number of planted trees per hectare). Taking the statistically fitted self-thinning line (Eq. 4) as a reference, the model of Vanclay and Sands (2009) was adjusted so that it joined with the self-thinning line (Fig. 1).

Fig. 1

Self-thinning frontier of cork oak forest and its modification according to Vanclay and Sands (2009) for three different initial stand densities (600, 1200, and 1800 trees ha⁻¹)

2.3 Models for simulating cork yield

Cork thickness at breast height in the first debarking was calculated with the following model (Sánchez-González et al. 2007d):

$$ {ct}_1=1.2677\left[1.8763\left(\frac{h-1.3}{h}\right)-2.7015{\left(\frac{h-1.3}{h}\right)}^2+2.2734{\left(\frac{h-1.3}{h}\right)}^3\right]\bullet d $$

(6)

where ct₁ is the cork thickness at breast height in the first debarking (mm), h is the total height of the tree (m), and d is the diameter at breast height over cork (cm).

The GADA formulation derived by Krumland and Eng (2005) from the Richards model was used to calculate the cork thickness after the first debarking (Sánchez-González et al. 2008):

$$ ct= CI{\left(\frac{1-{e}^{-0.04\bullet t}}{1-{e}^{-3.60}}\right)}^{\frac{2.43}{\frac{1}{2}\left(\ln (CI)+0.682+\sqrt{\ln {(CI)}^2+1.364\bullet \ln (CI)+9.365}\right)}} $$

(7)

where ct is the cork thickness at breast height (mm) t years after debarking, and CI is the cork index (mm). Figure 2 shows the increment of bark thickness as a function of CI and time since debarking according to Eq. 7.

Fig. 2

Increment in cock thickness after debarking for cork indices 20, 30, and 40 mm. Cork index is the thickness of the cork at breast height 9 years after debarking. The dashed horizontal line is the minimum cork thickness required for wine stoppers

Based on Natividade (1950), Montero and Vallejo (1992), and Sánchez-González et al. (2007d), cork thickness was assumed to vary along stem height as follows:

$$ {ct}_i= ct+\left(1.3-{h}_i\right)\bullet 3.33 $$

(8)

where ct_i is the cork thickness at distance h_i (m from ground). Debarking height was based on the instruction that debarking height should be equal to stem circumference at breast height multiplied by 1.5 (Montero et al. 2009):

$$ DH=0.015\bullet 2\pi d $$

(9)

where DH is the debarking height (m), and d is the diameter at breast height over cork (cm).

The thickness of the cork was also assumed to depend on the number of times the tree was debarked following an ascending–descending pattern so that the fastest growth rate was achieved after the 7th debarking (Machado 1944; Tomé et al. 1999). The bark thicknesses obtained from Eq. 7 were multiplied by the empirical coefficients presented by Tomé et al. (1999) (dots in Fig. 3).

Fig. 3

Variation of cork growth rate as a function of the number of debarkings (adapted from Tomé et al. 1999)

The weight of cork was calculated from:

$$ w=\pi \cdotp du\cdotp 0.5\left({ct}_{top}+{ct}_{bottom}\right)\cdotp DH\cdotp \rho $$

(10)

where w is the weight of the cork extracted from a tree (kg), du is the diameter at breast height under cork (m), ct_top is the cork thickness (m) at debarking height (DH), ct_bottom is the cork thickness at ground level, DH is the debarking height (m), and ρ is the density of the cork, assumed to be 251 kg m⁻³. Thicknesses ct_top and ct_bottom were calculated with Eq. 8.

In the case of the first and second debarkings, all the cork was of lower quality (quality 2). From the third debarking onwards, quality 1 cork was obtained if cork thickness was 27 mm or more (27 mm is the necessary thickness for stoppers). The volumes of quality 1 and 2 corks were calculated as illustrated in Fig. 4.

Fig. 4

Calculation of cork volumes in quality classes 1 and 2. Left: cork thickness at debarking height is greater than 27 mm. Right: cork thickness is greater than 27 mm only for a certain proportion of the debarking height

2.4 Optimization

The optimization algorithm employed in this study was the Hooke and Jeeves direct search method (1961). A set of decision variables was optimized. Decision variables are controllable variables that define the management schedule. In this study, the management schedule was defined by cork extractions, thinnings, and rotation length. The number of commercial thinnings tried was 0, 1, 2, and 3. In addition, it was assumed that a pre-commercial thinning was conducted 14 years after planting (Table 2). Preliminary optimizations showed that a higher number of thinning never increased profitability, which was due to the low density of the initial stands. Thinnings were uniform, and final fellings were shelter wood fellings in two stages with a time step of 10 years in between. The decision variables of this study were as follows:

• For debarkings:

  • Stand age at first debarking or number of years since previous debarking (for subsequent debarkings)

• For thinnings:

  • Stand age at first thinning or number of years since previous thinning (for subsequent thinnings)

  • Thinning intensity (% of number of trees)

• For final felling

  • Number of years since the last thinning

Table 2 Costs of silvicultural treatments in cork oak stands

The objective function maximized in optimization was the soil expectation value (SEV, in € ha⁻¹). Although the initial stands were 20 years old, all costs and revenues were discounted to the beginning of the rotation (planting year). It was assumed that exactly similar rotations are repeated to infinity. Management costs included planting, shrub clearing, pruning, and pre-commercial thinnings. The silvicultural operations that were assumed in calculations were based on Ovando et al. (2010) (Table 2). The used discount rates were 1, 3, and 5%.

Cork prices used in the study are based on COSE (2014). Cork price was 0.24 € kg⁻¹ for quality 2 cork and 1.2 € kg⁻¹ for quality 1 cork. The prices were varied to analyze the effect of cork price on optimal management. Quercus suber wood is rarely used as timber. Therefore, it was assumed that all the wood of cut trees was sold as firewood, with a stumpage price of 30 € per ton of dry wood. An entry cost of 160 € ha⁻¹ was assumed for every cutting.

Data availability

The data analyzed in the present study comes from the Second National Forest Inventory of Spain; the different files are available at: http://www.mapama.gob.es/es/biodiversidad/servicios/banco-datos-naturaleza/informacion-disponible/ifn2_descargas.aspx

3 Results

3.1 Optimal management of cork oak stand

The optimal management schedule included two commercial thinnings when cork index (CI) was 20 mm and no thinnings (in addition to pre-commercial thinning) when CI was 30 or 40 mm. Stands with higher CI had longer rotation lengths (Table 3; Fig. 5; Supplementary Table 1). Since cork generated much more income than wood, the SEV of the optimal schedule was closely related to CI (Table 3).

Table 3 Effect of cork index on the optimal management schedule when discount rate is 3% and site index is 14 m

Fig. 5

Optimal management schedules for stands with cork indices 20, 30, and 40 mm when discount rate is 3%, and the prices of quality 1 and 2 corks are 1.2 and 0.24 € kg⁻¹, respectively. The dots represent the debarking operations, and the dashed line represents the last regeneration cutting that was not optimized

The number of debarkings increased with CI. Debarking intervals tended to be shorter with higher CI, and the first and second debarkings took place earlier (Supplementary Table 2; Fig. 5). The interval between debarkings varied during the rotation, most clearly in CI 40 mm. The interval was the shortest in the middle of the rotation when the growth of cork thickness was fast. The shortest interval was 5 years and the longest as long as 25 years. The average debarking interval was 16.6 years for CI 20 mm, 10.5 years for CI 30 mm, and 8.2 years for CI 40 mm. The proportion of quality 1 cork increased with CI; the proportion of quality 1 cork was 24.4% with CI 20 mm and 62.2% with CI 40 mm. The mean annual cork yield over the whole rotation also increased with cork index, ranging from 0.10 t ha⁻¹ year⁻¹ (CI 20 mm) to 0.40 t ha⁻¹ year⁻¹ (CI 40 mm).

3.2 Effect of discount rate on optimal stand management

Decreasing discount rate led to longer rotation lengths while increasing discount rate increased the number of thinnings, one-thinning schedule being optimal with a 5% discount rate (Table 4; Fig. 6; Supplementary Table 3). A discount rate of 5% resulted in a SEV of − 1890 € ha⁻¹ while the SEV was 320 € ha⁻¹ with a 3% discount rate, suggesting that the internal rate of return of cork oak plantations was between 3 and 5%.

Table 4 Effect of discount rate on the optimal management schedule when cork index is 30 mm and site index is 14 m

Fig. 6

Optimal management schedules for stand 3 (cork index 40 mm) for different discount rates (1, 3, and 5%) when the prices of quality 1 and 2 corks are 1.2 and 0.24 € kg⁻¹, respectively. The dots represent the debarking operations, and the dashed line represents the last regeneration cutting that was not optimized

With a high discount rate (5%), the mean annual cork yield decreased in both quality classes as compared to optimal management with lower discount rates (Table 4). Increasing discount rate also led to fewer debarkings and longer debarking intervals (Supplementary Table 4; Fig. 6).

3.3 Effect of cork price on optimal stand management

When cork price was decreased by 30%, final felling occurred earlier while a 30% increase in cork price led to longer optimal rotation length (Table 5; Fig. 7; Supplementary Table 5). The mean annual cork yield increased slightly with increasing cork price. The number of debarkings increased with increasing cork prices and decreased with lower cork prices (Supplementary Table 6; Fig. 7).

Table 5 Effect of cork price on the optimal management schedule when cork index is 30 mm, site index is 14 m, and discount rate is 3%

Fig. 7

Optimal management schedules for stand 3 (cork index 40 mm) for cork prices of quality 1 and 2 corks when discount rate is 3%. The dots represent the debarking operations, and the dashed line represents the last regeneration cutting that was not optimized

3.4 Effect of site index on the optimal stand management

Site index (SI) had a clear effect on the optimal management of cork oak stands. Although SI does not affect the growth in cork thickness, it does influence the cork weight production via the growth of stem circumference. The management of stand 4 (SI 8 m) was not profitable under the assumed economic parameters and silvicultural operations (SEV was − 2684 € ha⁻¹). The optimal management of this stands involves longer rotation length and one thinning treatment (Table 6; Fig. 8; Supplementary Table 7). Besides, the growth of the trees was so slow that the first debarking operation took place as late as in year 101 (Supplementary Table 8; Fig. 8), allowing only few debarking operations during the rotation. The mean annual cork yield was 0.14 kg ha⁻¹ year⁻¹, of which only 27.1% was quality 1 cork (Table 6).

Table 6 Effect of site index on the optimal management schedule when cork index is 30 mm, site index is 14 m, and discount rate is 3%

Fig. 8

Optimal management schedules for stand 4 (SI 8 m) and stand 2 (SI 14 m) with cork index 30 mm when discount rate is 3%, and the prices of quality 1 and 2 corks are 1.2 and 0.24 € kg⁻¹, respectively. The dots represent the debarking operations, and the dashed line represents the last regeneration cutting that was not optimized

3.5 Comparison with official recommendations

The shortest debarking interval obtained in optimization was 9 years for CI 20 mm, 7 years for CI 30 mm, and 5 years for CI 40 mm. Since some of the optimal intervals for cork indices 30 and 40 mm were shorter than the shortest interval of the official recommendations (9 years), we analyzed how much profitability would decrease due to the adoption of the minimum debarking interval of 9 years. This comparison was done by using a 3% discount rate and the baseline cork prices.

The results showed that shorter minimum debarking interval improved profitability; the contribution of cork to the total SEV increased from 265 to 320 € ha⁻¹ in CI 30 mm when the constraint for 9-year interval was removed and 986 € ha⁻¹ to 1084 € ha⁻¹ in CI 40 mm (Table 7). The optimal rotation length and mean annual cork yield were not much affected (Table 7; Supplementary Tables 9 and 10; Fig. 9).

Table 7 Effect of the constraint for the shortest allowed debarking interval (4, 7, or 9 years) on the optimal management schedule when cork index is 30 mm, site index is 14 m, and discount rate is 3%

Fig. 9

Comparison between the optimal management obtained when the official recommendation for minimum debarking interval (9 years) is followed and when the shortest allowed debarking interval is 7 years with CI 30 mm and 4 years with CI 40 mm. Site index is 14 m and discount rate is 3%. The dashed line represents the last regeneration cutting that was not optimized

4 Discussion

The optimal production of NWFPs from the economic point of view has been a focus of attention in last years. The positive contribution of NWFPs (mushrooms, berries, pine nuts, etc.) to the profitability of forestry has been proved by several studies that have used stand level optimization techniques similar to the one applied in the present study (Palahí et al. 2009; Pasalodos-Tato et al. 2009; Miina et al. 2010; de Miguel et al. 2014; Pasalodos-Tato et al. 2016). The present study focused on the management optimization of cork oak stands. Our results show, among other things, that it is not optimal to debark at regular interval. Optimal debarking intervals are site-specific, and they may be shorter or longer than the currently recommended range of 9–14 years. Paulo and Tomé (2017) used constant debarking intervals (from 9 to 14 years) to examine the profitability of different management schedules of cork oak stands without testing the possibility of altering the debarking interval during a rotation. Their results showed that profitability was lower on poor sites independently of debarking interval, and the use of high discount rate decreased the most profitable debarking interval.

The results of the present study are in line with those found for other NWFPs: high price of the NWFP often leads to delayed thinnings and longer rotations, the effect resembling that of decreasing discount rate (Pasalodos-Tato et al. 2009, 2016). The optimal schedule depends on cork prices, discount rate and cork index. The number of debarkings increases with increasing cork index (Table 3; Fig. 5), and the interval between debarkings is shorter (Supplementary Table 2; Fig. 5). Decreasing discount rate increases the number of debarkings (Table 4; Supplementary Table 4; Fig. 6). Increasing cork price leads to longer rotations and higher number of debarkings (Table 5; Supplementary Tables 5 and 6; Fig. 8). However, the mean annual cork production is not much affected. Stands growing on poor sites have longer optimal rotation lengths. Cork oak forestry is rarely profitable on poor sites (Table 6; Supplementary Tables 7 and 8; Fig. 8). These stands should be managed in a more extensive way and for other types of NWFP, for instance acorns to feed wild game or cattle.

To assess the profitability of the whole cork oak system, other benefits traditionally related to cork oak stands may also be considered, namely grazing, hunting, mushroom picking, or honey production. However, lack of quantitative information on these products prevented their inclusion in the present study. In Spain, Ovando et al. (2010) used cost-benefit analysis to explore the profitability of afforesting marginal lands with cork oak. They included other sources of income apart from firewood and cork, namely hunting, grazing, land renting for dry crops and afforestation, and silvicultural subsidies. The management proposed in that study used a more intensive management schedule (higher number of thinning and pruning treatments) and fixed treatment prescriptions. Ovando et al. (2010) found that the profitability of afforestation options was highly dependent on subsidies.

The model for cork thickness (Sánchez-González et al. 2008) was fitted for debarking intervals between 9 and 12 years. Because the optimal debarking intervals were often shorter than 9 years or longer than 12 years, there is a need to refit the model using data with more variation in debarking interval. It would also be useful to analyze whether short debarking intervals would increase damages to the trees. Although this topic has already been analyzed, no definitive conclusions can be drawn based on earlier results (Oliveira and Costa 2012; González-García et al. 2013; Davila et al. 2014).

5 Conclusions

The management of cork oak forests is a challenging task since many factors influence the optimal management. The present study optimized the joint production of wood and cork to improve the profitability of the management of cork oak stands. The results show that both economic factors (discount rate and cork prices) and stand properties (cork index, site index) affect the optimal management of cork oak stands. The optimizations showed that adopting debarking intervals shorter than official recommendations may improve the profitability of the management of some cork oak stands. Debarking intervals should vary during the rotation so that shorter intervals are used when cork growth is fast. However, the effect of these changes on the growth and survival of cork oak trees should be studied before implementing substantial changes in large areas. The results also suggest that a shift from a static management system towards optimized and more flexible management may improve profitability.