Using the SUBER model for assessing the impact of cork debarking rotation on equivalent annual annuity in Portuguese stands

Aim of study : Use the SUBER model to evaluate the influence of the cork debarking rotation period (CDR) on equivalent annual annuity (EAA) value. Area of study : Nine simulated stands, varying in site index (14.4, 15.6, 17.1) and cork quality characteristics (high, medium, low). Material and methods : EAA values were computed considering CDR periods varying from 9 to 14 years, two contrasting structures of cork prices (high and low cork price scenarios), and three discount rate values (0.5%, 2% and 5%). Main results : For discount rates of 0.5% and 2% the impact of different CDR on the EAA is similar. In stands characterized by high to average site index values or high to medium cork quality characteristics, CDR of 9 and 11 years are associated with similar values of EAA. The variation of the CDR in stands characterized by low site index values and/or low cork quality characteristics did not have a relevant effect on the variation of EAA. For the simulations carried out with a discount rate of 5% the EAA decreases with the increase of CDR, indicating that the minimum legal value of 9 years for CDR should be applied. Research highlights : In stands characterized by high to average site index values or high to medium cork quality characteristics, a delay in the debarking may result in a significant increase of cork thickness and, as a result, of cork price. Detailed knowledge of cork and stand characteristics and updated information on cork prices structure and values are essential for the best usage of management tools such as the SUBER model, which can contribute to the decision-making process concerning the debarking operation.


Introduction
Because of the importance of cork oak trees in montado and dehesas ecosystems where cork oak (Quercus suber L.) is present, several management options are regulated by national legislation in some countries along its distribution area. One example is the interval between two consecutive cork extractions occurring in the same tree, here designated as cork debarking rotation (CDR), which, according to Portuguese legislation, should be of at least 9 years. Other countries have established different values for this management option. For example in Spain, as in Portugal, the debarking can take place after a minimum period of 9 years, but in Italy and Tunisia this minimum period is extended to 12 years.
The cork industry transforms this raw material into different products, but the economic feasibility of the whole sector is dependent on the production of stoppers of natural cork for use in the bottling of wines (Pereira, 2007;Corticeira Amorim, 2013). Therefore, the suitability of the raw material for this product establishes its commercial value. Two features are important when defining the value of the cork for the industry of cork stoppers: the thickness of the cork plank and the cork quality. Cork thickness can be measured either before or after boiling, but the latter is the reference for the industry and will be considered throughout this research (see details of this industrial process in Pereira, 2007). Cork quality is visually evaluated by an operator who, considering the cork porosity and the presence of cork defects, classifies the cork into 7 classes: from quality class 1 (best quality) to quality class 7 or refuse (worst quality). Due to the great variation of cork thickness (e.g. Paulo et al., 2016) and cork porosity (e.g.  between trees in the same stand, and since different combinations of these two variables result in different prices, average cork price is not a suitable variable for characterizing the cork produced in a particular stand. Instead, we suggest the usage of a two -dimensional structure of cork prices that includes the proportions of both cork thickness and cork quality classes that characterise the stand for determining revenues associated with cork extraction . In practice three aspects mainly influence the decision of the forest manager concerning the debarking: evolution of the structure of cork prices, labour costs and weather conditions during the cork growth period. For instance, an increase in the prices of thin corks (< 22.6 mm) may lead to the anticipation of cork extraction, while an increase in the price of thicker corks (> 40.6 mm) will make landowners postpone cork debarking one or two years, in the expectation of increasing the percentage of thick cork. In the years between 2000 and 2010, cork price structure experienced considerable fluctuations, generally showing a significant decrease in the value for most of the cork assortments. This was more evident in cork classified as poor and medium quality classes (APCOR, 2014). Even if in recent years the prices for some cork quality classes have increased and have been more stable, there is an increasing number of forest managers extending the CDR in their stands in an attempt to increase the value of the extracted product and therefore increase the resulting income.
In the matter of cork oak stands, there are three available models that can be used for this task: the SUBER model (Tomé, 2004;Paulo, 2011;Tomé & Faias, 2014) and the CORKFITS model (Ribeiro & Surový, 2011) for Portuguese cork oak stands, and the ALCORNOQUE model (Sánchez-González et al., 2007;Tomé & Faias, 2014) for Spanish cork oak stands. The SUBER model has been selected because: • the model is included in the sIMfLOR platform (Faias et al., 2012) accessible for download at www.isa.ulisboa.pt/cef/fctools • the model simulates annual cork growth based on tree level cork growth index values • the model considers tree cork quality at tree level • the model includes an economic module that uses cork prices structure The objectives of the present work are: (i) to analyse the impact of site index, cork quality and thickness distribution and cork price structure on the equivalent annual annuity (EAA); (ii) to analyze the effectiveness of the extension of the CDR on the increase of EAA values for stands with different characteristics.
Although the simulation unit of the SUBER model is the tree, for a time step of one year, its outputs are produced at the stand level. Since forest management is made at this scale, this approach is considered adequate for the purpose of the present work. A multistand approach addressing cork oak stands in Portugal has been carried out by Borges et al. (1997) using previous simpler versions of the SUBER model that did not allow the simulation of different CDR. Palma et al. (2015), using the present version of the SUBER model, focused on an approach for the establishment of adaptive management plans to climate change on a regional approach, but the solutions found for the optimal CDR calendar are, in practice, difficult to implement in the current land management scenario. These studies are relevant at a strategical level, e.g. definition of forest policies or regional prediction of cork availability, but are not appropriate to support the landowner's decisions.
Revenues and costs resulting from pastoral or agricultural activities are not accounted for either, since the tree and cork growth modules included in the SUBER model do not take into account the influence of these practices on tree and cork growth. These have been considered by Campos et al. (2007aCampos et al. ( , 2007bCampos et al. ( , 2008, Costa et al. (2010) and Fragoso et al. (2011) using alternative simulation tools that do not account for the increase of cork thickness with increasing cork age, the existence of cork prices structures or the differences in cork growth and quality between trees in the same stand.

Tree data
Tree data was collected on three existing permanent plots established in three young cork oak plantations. The stands present similar tree density and management operations up to the moment of the measurement (Table  1). A characterization of the stands at the year of the measurement is presented in Table 1, including site index value for a base age of 80 years estimated with the model developed by Sánchez-González et al. (2005).
The three selected plots are part of a larger database that includes measurements made in 42 permanent plots established in young plantations in Portugal, which present site index values ranging from 11.5 to 18.4 m (reference age of 80 years), with an average value of 16.6 m.

Cork data
The SUBER model needs the input of the empirical distribution values for the cork thickness and for the cork quality variables. These two distributions may be applied to the model from two alternative sources: i) at the tree level: information collected during a forest inventory and cork sampling procedures; ii) at the stand level: information derived from forest inventory or other source but not only available at the stand level (see details in Paulo, 2011). In the latter option the model uses the distribution values to assign cork growth index and cork quality class values to each tree in the stand, using Monte Carlo simulation (e.g. Gilks et al., 1996). In the present research the information provided refers to frequency distributions.
As no relationship has been found between cork quality and cork growth index (Pereira, 2007) the distributions considered as model inputs were selected from a large data set including 187 cork inventories made in different stands for estimation of the cork value before the extraction. Each inventory included more than 30 cork samples. The average cork price for each of the 187 cork inventories was computed using the real cork prices structures from 2009 (Associação dos Produtores Florestais de Coruche, personal communication). The ones resulting in a higher, medium and lowest cork price value were used for the simulations and assumed   (Table 2). It is important to notice that passing from high to low cork prices generally represents a decrease in the fraction of good quality corks (1, 2 and 3 cork quality classes) and/or thick corks (>41 mm cork growth index class), and simultaneously an increase in the fraction of low quality corks (6 and refuse cork quality classes) and/ or thinner corks (< 23 mm cork growth index class). In detail the variation of cork price may be more complex, since it results from the simultaneous combination of the different classes of thickness and quality classes.
The SUBER model includes the  system of equations for predicting the evolution of individual tree mature cork caliper over time, and the Paulo & Tomé (2010) system of equations for estimating cork biomass production depending on cork age. Cork quality class of each tree is assumed to be constant along the years of the simulations and regardless of the cork age at the debarking by the model, as no relation between cork quality and age has been demonstrated in the literature. This assumption is reasonable if cork defects, for instance due to insect presence (e.g. Coroebus undatus, Crematogaster scutellaris), are discarded. These assumptions allow the model to determine the cork value, assuming different cork debarking rotation periods on each simulation.

Forest management options and silvicultural operations
The main aspects regarding the management of cork oak stands in Portugal are regulated by national legislation (Decree-Law n.º 169/2001 and Decree-Law n.º 155/2004). These determine: • The minimum size of the tree perimeter at breast height over cork (p) after which the first debarking is allowed: 70 cm.
• The minimum interval between two consecutive cork extractions on the same tree, designated by the cork debarking rotation (CDR): 9 years.
• The maximum values for debarking coefficient (ratio between debarking height and tree perimeter over cork): 2 for the first debarking, 2.5 for the second debarking and 3.0 for the following debarking operations.
Regarding the forest management operations defined for the simulations, the assumptions take into account both national legislation and common practices for high natural and cultural value systems management in Portugal. Details on the forest management operations are presented in Table 3.
The maximum value of 14 years, considered for the interval between two cork extractions, is based on the fact that this is the longest period recorded in  the database of the Forchange research group, that was used by  and Paulo & Tomé (2010) models included in the SUBER. Intervals longer than 12 years between debarking operations are only encountered in very particular situations and in a reduced number of stands.

Economic data
Two real cork prices structures were considered for characterizing two distinct moments in the cork market in recent decades: the highest and lowest average cork prices, which occurred in the years 2000 and 2009, respectively. These distinct cork prices structures are the basis for the definition of two distinct scenarios: the first characterizing a favourable management scenario and the second an unfavourable management scenario. An average price scenario was not considered since this is difficult to define when considering a cork price structure instead of one single average cork price. For instance, the definition of a cork price structure characterized by an average cork price for all the cork thickness and quality classes would be completely unrealistic and not useful for fulfilling the objectives of the present study. These market fluctuations have originated from a combination of drivers, such as the evolution of the world wine market, the development of new cork applications and the improvement of methods of conservation of the raw material in the industry in order to cope with the fluctuations of material supply.
Although the classification of the cork is made by the industry considering five thickness classes and seven quality classes, cork prices are not differentiated for the resulting 35 combinations. For example, cork quality classes 1, 2 and 3 are typically attributed the same price, as the use given to this cork material is generally the same (e.g. . This is reflected in Tables 4 and 5. The values presented in these tables already include the costs from the debarking operation. Virgin cork price is equal to the price attributed to refuse cork. When comparing the cork price values from the two tables, although referring to two different years, the differences regarding the two cork price structures and values are clear. Poor and medium cork quality classes were strongly affected by the decrease in price in 2009, but good cork even revealed some increase in the prices from 2000 to 2009. Another important aspect derived from these tables is the decrease in the number of cork price categories, going from 9 to 6 between 2000 (Table  4) and 2009 (Table 5).
In order to compare the simulations, the equivalent annual annuity (see for example Davis et al., 2001) is computed by the SUBER model as EAA = (r NPV) / [ 1 -(1+r) -n ] Prices already include the costs from the debarking operation (see Table 6). Prices already include the costs from the debarking operation (see Table 6).
where EAA is the equivalent annual annuity (€ ha -1 year -1 ), NPV is the net present value (€ ha -1 ), r is the discount rate assumed as 0.5%, 2% or 5% depending on the simulation, and n is the number of years of the simulation (80 years for all the simulations). NPV (see for example Davis et al., 2001) was calculated as where R t is the revenue for year t, C t are the costs for year t, r and n are as previously defined.
NPV value includes the revenues resulting only from cork production. No revenues resulting from wood selling, cattle raising, or any additional activity (e.g. subsidies) were considered.
Costs from all silvicultural operations previously defined in section 2.3 are included in the C t value. Costs are considered constant for the period of the simulation and were obtained in the CAOF (Comissão de Acompanhamento das Operações Florestais) reference matrices for 2013/2014, available in http://www.anefa. pt/. Since the objective of the paper is the evaluation of the impact of site productivity, cork growth index and quality distributions and cork prices structure on cork debarking rotation for long-term management of new cork oak plantations, the stands were assumed to be similar regarding site conditions, namely slope. Table 6 presents the value for each considered operation.
Costs and revenues associated to sanitary thinning operations are both considered equal to zero. Although some revenue may come from the selling of woody material with larger diameter, this is often consumed in domestic use. Other biomass residuals (small branches and leaves) are frequently discarded in the field, sometimes incorporated in the soil.

Simulations
Through the combination of tree data collected in three young existing stands (undebarked and for this reason without information on cork characteristics), with real cork quality and cork growth index frequencies distribution, nine simulated stands were defined for simulations. Growth and cork production of each stand was simulated until the stand reaches 80 years of age, considering two contrasting cork price structures and cork debarking rotations (CDR) ranging from 9 to 14 years. Simulations were carried out using the SUBER 5.0 model (Paulo, 2011;Tomé & Faias, 2014). A total of 324 simulations were carried out (3 stands x 3 cork distributions x 2 cork price structure x 6 CDR alternatives x 3 alternative discount rate).

Results
Figures 2, 3 and 4 present the results from the simulations carried out considering discount rate values of 0.5%, 2% and 5%, respectively. They present the values of the equivalent annual annuity (EAA) obtained for each simulated stand and each cork price structure considered, as a function of the cork debarking rotation (CDR). The analysis of the effects of different site indices, cork quality and thickness distributions and cork price structures on the CDR that maximizes equivalent annual annuity (EAA) are presented in the discussion section.

The impact of different site index of the stand on the EAA
The impact of site index on EAA can be observed by comparing the three curves presented on each of the graphics in Figures 2, 3 and 4 separately. As expected the EAA values are smaller in sites with lower site indices. This reduction results from lower tree growth rates in sites with lower site index value, according to Tomé et al. (2006) diameter growth model that is used by the SUBER model. Since cork production (cork weight) is related to tree perimeter at breast height (Paulo & Tomé 2010), stands with lower site index are characterized, for the same age, by trees with lower perimeter and consequently lower values of cork produced for a same CDR. These results point out the importance of the Table 6. Cost of the silvicultural operations considered along the simulations (see Table 3). Includes all valid legal taxes. Source: Comissão de Acompanhamento para as Operações Florestais (CAOF) reference matrices for 2013/2014, available in http://www.anefa.pt/

Operation Unit Cost
Mechanical control of the understory vegetation with mounted knifes or chains € ha -1 195.57 Formation pruning € tree -1 0.45 Cork debarking * € @ -1 3.50 * Value already deducted in the cork prices presented in Tables 4 and 5 location of cork oak stands, since the relationship of local edafoclimatic variables and tree productivity has been demonstrated on several occasions (Hidalgo et al., 2008;Lacambra et al., 2010;Paulo et al., 2015). When looking at the relationship between the curves, the ones denoting high and medium site indices are nearly parallel in all the graphics, indicating that for stands located in good sites the site index does not have a strong impact on the definition of the CDR that maximizes EAA. For scenarios using discount rates of 0.5% and 2% (Figures 2 and 3) CDRs of 9 or 11 years provide similar results. For scenarios using the discount rate of 5% (Fig. 4) a decreasing tendency of the EAA with CDR is observed, showing that CDR of 9 years should be implemented if EAA maximization is the main objective of forest management.
The simulations of the low site index stands show that, despite the discount rate considered, EAA values were low and not very sensitive to the variation of the CDR. For these stands the importance of the multifunctional management of the montado is clear, since the income resulting from other products such as game, agriculture production, wild mushrooms picking etc, not accounted for in the SUBER model, will have an important input in the final economic return for landowners (Campos et al., 2007a). Subsidies derived from the common agricultural policy are also expected to impact economic return (Fragoso et al., 2011). Without this economic support some montado areas may be driven to land use conversion, depletion of oak trees due to tree aging and lack of tree renewal, or even land abandonment (Campos et al., 2007b;Campos et al., 2008).

The impact of cork growth index and cork quality distributions
The effect of the cork quality and cork growth index distributions can be assessed by comparing the same type of curve along each of the three graphics named with prefix 'H', 'M' or 'L', separately in the upper and lower lines of graphics presented in figures 2, 3 and 4. As expected a decrease in the fractions of thicker and good quality cork, observed when moving from graphic H to M and then to L, results in a significant decrease of EAA values. This result is related to the fact that these two variables are the ones determining cork price, and consequently the income resulting from the cork extraction.
The impact of the cork growth index and cork quality distributions on the CDR that maximizes the EAA was more evident in the simulations carried out with discount rates of 0.5% and 2% (Figures 2 and 3). For high and medium average cork an increase of the CDR to 11 years may result in very similar or even higher EAA values. This suggests that postponing the extraction of the cork for one year is not compensated by an increase of the cork thickness, since this does not create a shift in the price to a higher cork price class. This is only obtained when the debarking is delayed by two years, increasing the CDR to 11 years.
Stands characterized by low and medium cork are characterized by a reduction of the EAA value for increasing CDR. Although the increase of the CDR gives rise to an increase in cork thickness, low quality cork characteristics are related to lower cork price values, in particular cork quality classes 5, 6 and refuse (see tables 4 and 5), the ones more represented in the simulated stands presented in graphics 'L'. For a discount rate of 5% (Fig. 4) the impact of extending CDR to 11 or more years is reduced since a decrease or stabilization of the EAA value is observed in all curves, suggesting that carrying out debarking operations are best practised in periods of 9 years.
Results demonstrate the importance of the cork growth index and quality distributions on the definition of CDR. In order to provide accurate data to the SUBER model and improve the accuracy of the simulated outputs, prior forest inventory and cork sampling should be carried out (cork sampling). This is already a standard operation in some of the major cork production regions of Portugal, carried on with specific sampling designs that optimise the information about cork value .

Impact of different cork price structures
The effect of the cork price structures on EAA can be observed by comparing the same type of curve (full, broken and dot lines) of the upper and lower graphics (suffix '2000' and '2009') in each of the Figures 2, 3 and 4.
As expected, and regardless of the discount rate taken for the simulation, the values of the EAA are strongly dependent on the cork price structure, with a clear reduction in the resulting values when moving from cork prices structure for the year 2000 to the one for the year 2009 (high and low cork price scenarios respectively). It is also observed that the same curves presented in the upper (2000) and lower (2009) graphics are nearly parallel. This shows that the CDR values that maximize the EAA are not related to the cork price structure.

Decision on cork debarking rotation
Increasing income resulting from cork extraction may be achieved by: i) reduction of labour costs; ii) increasing cork price; iii) increasing the amount of cork extracted by increasing debarking coefficients and consequently debarking intensity; iv) increasing cork thickness and/or the percentage of good quality cork. The evolution of the first two is market driven, and forest managers are limited in their influence. Still, due to the impact of cork price structures on the final income, annual assessments considering the historical records of cork price structures should be carried out when cork age is 7 or more years, in order to provide additional information relevant to the decision regarding CDR. This can be achieved using the SUBER model, which allows the realization of simulations considering cork thickness and cork quality data at tree level, different CDR and different cork price structures.
In the matter of the amount of cork extracted, it is common practice not to debark trees up to the maximum value allowed by law. Although the relationship between debarking intensity and tree growth and resilience is still an area in need of research (Oliveira & Costa, 2012), landowners are usually conservative concerning debarking intensity, and even adapt it to the tree growing and sanitary conditions. Decision on CDR should consider not only the variables considered in this research, but also growth conditions, taking into account the importance of cork harvesting and climate for stem radial growth of Quercus suber L. (Oliveira et al., 2002;Leal et al., 2008). It is frequent, when the price of thinner cork is low, for landowners to decrease the debarking height since cork thickness is known to decrease along the stem height (Montero, 1992) and simultaneously cork extraction costs are known to increase when debarking at upper tree heights and along tree main branches.
Since the relationship between management practices and cork quality is still not well established, it is difficult to influence cork quality. Although cork quality is known to be largely related to genetic characteristics (e.g. Teixeira et al., 2014), further research on the development of reproduction material is needed. However, cork thickness is a result of cork age (determined by the applied CDR) and influenced by climatic conditions during the growing period, such as precipitation and temperature (e.g. Caritat et al., 2000;Paulo et al., 2016).
The large range of variables influencing cork production implies that any decision regarding CDR is a complex task. Decision support tools such as the SUBER model can provide additional insight during this process. It is suggested that this should be made in a two-stage procedure: 1. Initial assessment of the evolution of the EAA for extending CDR periods from 9 to 14 years at tree level: this should be made using the SUBER model, following a similar approach to the one presented here. This long-term simulation will give a first definition of the best CDR according to the stand characteristics. This analysis can also be incorporated in forest management plans. 2. Around 7 to 9 years of cork age, using updated data on cork quality and thickness distributions and recent cork prices structure information, simulation of cork thickness evolution and expected EAA value, for CDR values ranging from 9 to 14 years. This can also be made using the SUBER model. It is expected that this simulation will contribute to the decision on whether to debark or not. If the decision is made to postpone the debarking, new simulations can be done the following year if up-to-date data on cork prices structure is available.

Conclusions
The three considered drivers (site index, cork quality and thickness distribution, and cork price structure) influence the EAA value.
The two most influencial ones, which allow maximizing the EAA, are site index and the cork quality and cork growth index distributions. The variation of EAA due to the CDR extension was evidenced in scenarios considering discount rates of 0.5% or 2% of. For these scenarios similar values of EAA were obtained by postponing the debarking for CDR values to 11 years for stands characterized by high or medium site index and high or medium cork characteristics (quality and growth index). Stands characterized by low site index or low cork characteristics presented lower, nearly constant, and in some cases negative EAA values, irrespective of the CDR considered in the simulations.
For 5% discount rates the extension of the CDR period systematically originates a reduction of the EAA value, suggesting an increase in the number of cork extractions instead of an increase in cork thickness as a result of the CDR postponement.
Growth and yield models such as the SUBER model can be helpful tools for decision makers in the decision to undertake or postpone a debarking operation. The need for accurate and updated data on cork stand characteristics and cork prices is crucial for the accomplishment of similar studies.