Modelling non-wood forest products in Europe : a review

Non-wood forest products (NWFP) like cork, edible mushrooms, pine nuts, acorns, resins, medicinal plants, and floral greens, among others, provide important recreational and commercial activities in the rural forested areas of the world, their production being in certain regions more profitable than traditional timber production. Despite the importance of non-wood forest products and services, forest management and planning methods and models in Europe have been traditionally wood production oriented, leading to a lack of tools that could help foresters to manage for optimizing these products. In the present work we’ll discuss the main factors and challenges limiting the development of classical empirical models for NWFP, and we will review the existing models for the main NWFP in Europe: cork, pine nuts, berries, mushrooms and resins.

beauty, biodiversity, berries and mushrooms were already regarded as the most important values of forests.Table 1 shows a general overview of the most important NWFP in Europe, together with their estimated utilization (harvested production) and commercial value.Despite the importance of non-wood forest products and services, forest management and planning methods and models in Europe have been traditionally wood oriented.This is due to the fact that most forest management models and silviculture schemes were developed to ensure a sustained and maximum yield of timber.Nowadays, it is well known that maximising wood production does not necessarily maximise the provision of relevant non-timber forest products and services (Calama et al., 2007;Bonet et al., 2008;Palahí et al., 2009;Miina et al., 2010); in some cases it is possible to identify some incompatibilities among silvicultural practices for promoting non-timber products and timber yield and quality (see Table 1).At the same time, forest management planning has developed approaches, models and tools that can explicitly and efficiently consider non-wood forest products as management objectives (Pukkala, 2002a(Pukkala, , 2002b)).Nevertheless, in order to apply such approaches and tools, there is a major bottleneck: models for assessing, in a quantitative way, the production of non-timber products in different forest situations and for different management schedules are required.One of the ways to develop such new models is based on statistical methods, using empirical data from permanent plots designed specially for inventorying and modelling non-timber products.In any case, there are actually very few models for non timber products in Europe, due to the lack of systematically collected data (Ihalainen and Pukkala, 2001), together with some challenges which make it difficult to develop predictive models for non-timber products.

The challenge of modelling NWFP
Together with the difficulty in getting adequate data for modelling NWFP, some other complexities emerge when developing predictive models for non-timber products.These difficulties concern both data requirements and modelling methodologies, and can be summarized as follows: -NWFP data are seldom normally distributed, showing asymmetrical distributions with a heavy, right 70 R. Calama et al. / Forest Systems (2010)  tail and abundance of zero and close-to-zero observations.Given the nature of the inventories, spatial and temporal correlations among observations are also common.
-High inter-annual and spatial variability, which require the systematic collection of long-term annual data series as well as a good characterization of weather/ climate conditions and soil attributes.
-Site index constitutes the main driver of empirical growth and yield models, devoted to timber production.However, traditional site index (defined as the dominant height of the stand at a given index age) does not usually indicate the potential productivity of NWFP.Other ecological attributes, such as topography, climate or soil conditions, are required to estimate the site potential productivity.
-Compared to timber, NWFP have a short harvesting period, and, in several occasions, the products perish soon after this period.In other occasions, these products are threatened by unauthorized collectors.Therefore NWFP data collection requires intensive monitoring systems, with a great effort of field work.
-In some NWFP, as mushrooms or berries, there exist a great variety of species that are harvested and marketed jointly, since they provide a similar end product, in spite of showing large differences in ecological and management requirements that should be taken into account in the modelling task.
-The production of NWFP like fruits, barks or resins involves several complex physiological processes such as allocation of carbohydrates to growth, flowering and fruiting, or other traits.These processes are poorly known, making it difficult to develop process based models for this purpose.
-When referring to NWFP, quality requirements are usually as important as the total amount of production.Therefore, models for NWFP should also consider quality attributes as response variables.
As a synthesis of the abovementioned challenges, Table 2 includes how these limitations can influence the modelling of the main NWFPs in Europe, together with the proposed requirements for data collection and modelling methodology.
In spite of the detected challenges and limitations for modelling NWFP (Tables 1 and 2), recent efforts have been carried in Europe in order to develop models for the main NWFP (Table 3).Modelling strategies have varied from pure empirical to process-based models, attaining spatial and temporal scales ranging from the tree to the region, and from a year to the rotation length.Modelling techniques go from ordinary regression analyses to multilevel mixed models, geostatistics and time-series analyses.In addition, in certain cases, the lack of empirical measurements has been overcome with expert modelling.

Existing models for the main European NWFP's
The following paragraphs are devoted to briefly describe the existing models for some of the most important NWFP in Europe: cork (Quercus suber), Pinus pinea nuts, edible mushrooms, berries and resin.

Cork
Cork is a thick and continuous layer of suberised cells, produced by the meristematic cork cambium (or phellogen), which makes up the external envelope of the stem and branches of the cork oak trees.It is a natural product that can be extracted without compromising tree vigour, since the tree is able to rebuild a new cork layer in reaction to cork extraction.Cork is usually commercialized on the basis of weight, the weight of cork extracted in one harvest (from then on designated by cork weight) being the ultimate objective of cork modelling.Other variables that are important for cork modelling are cork growth and cork thickness.
Cork oak is debarked during the growing season, and the growth of new cork starting shortly after the extraction.This is the reason why the first and last growth rings are not complete.Figure 1 shows an example of a cork with 9 years of age (8 complete growth rings).The first cork growth ring, located adjacent to the cork back, is usually smaller than the subsequent rings as it corresponds to the growth of cork in the growing season during which cork was extracted, after debarking.Therefore it is not a complete growth ring and is usually called the 1st half ring.Similarly, the last cork growth ring corresponds also to a partial growing period, before cork extraction, and is called the last half ring.The production of new cork is more intense in the years immediately after debarking and decreases from then on.The caliber (thickness) of a cork plank with tc years (cork age is designated by tc in order to distinguish from stand/tree age that is usually designated by t) is equal to the sum of the thickness of 72 R. Calama et al. / Forest Systems (2010)    the tc-1 complete rings plus the thickness of the two half rings (initial and final) and of the cork back thickness.For the sake of clarity, the term cork caliber (cc) is used here for the total cork thickness while the thickness of complete rings is designated by thickness of cork complete rings or simply by cork thickness (ct) as some authors did (Sánchez-González et al., 2008).The first cork weight models were simple or multiple linear models developed in Portugal in the 50's (Natividade, 1950;Guerreiro, 1951;Alves, 1958;Alves and Macedo, 1961;Ferreira et al., 1986) using fresh cork weight as dependent variable and tree diameter at breast height, measured without cork, stem height, crown width, cork thickness, number of debarked branches and debarking coefficient as regressors.The first models developed for Spain (Montero, 1987) used the same dependent variable.Some problems were detected when more recent researches focused on the development of cork weight models (e.g.Ribeiro and Tomé, 1992;Vázquez andPereira, 2005, 2008;Paulo and Tomé, 2010): 1) the dependent variable should be cork biomass (dry weight); 2) serious problems of heterocedasticity and non-normality usually occur in cork weight models; 3) data are usually obtained as an hierarchical data structure (with trees sampled inside plots located in regions); 4) cork age should be taken into account.A very detailed review, including a critical appraisal of all these issues can be seen in Vázquez and Pereira (2008).
Since the late 90's large effort has been put on the development of improved models for cork weight pre- diction (Ribeiro and Tomé, 2002;Tomé, 2004;Vázquez and Pereira, 2005;Paulo and Tomé, 2010).With the exception of Ribeiro and Tomé (2002), that used air dried cork weight, all the other recent models used fresh cork biomass as the dependent variable.The problem of heterocedasticity was overcome by computing the logarithm of the allometric model (Ribeiro and Tomé, 2002;Vázquez and Pereira, 2005) or by fitting the model using non-linear weighted regression (Tomé, 2004;Paulo and Tomé, 2010).The non-normality of the errors was taken into account by Tomé (2004) through robust estimation applying the Huber function to reduce the influence of data points presenting high fitting errors in the analysis (Myers, 1986).The hierarchical structure of the data implied the inclusion of random effects in the more recent cork weight models (Vázquez and Pereira, 2005;Paulo and Tomé, 2010).
Cork age was taken into account by Vázquez and Pereira (2005) through a dummy variable that indicates if cork age is 9 or 10 years.This method does not allow the simulation of cork weight for corks with less than 9 or more than 10 years, which is very important for optimization of forest management and planning of cork extraction.Tomé (2004) presented a methodology that allows for the prediction of cork weight for different cork ages.This methodology, applied and improved by Paulo and Tomé (2010), is supported by the homogeneity of the density of the cork tissue between the inner and outer cork rings (Barbato, 2004) as opposed to cork back that has about three times the density of cork: 150-225 kg m -3 in cork layers and 495-608 kg m -3 in cork back (Pereira, 2007).The method for prediction of cork weight for a cork with t years of growth (wcm t ) is based on the estimation of cork biomass if the cork had 9 years of growth (wcm 9 ).It includes four steps: 1. Estimate tree cork biomass for 9 years of age (wcm 9 ).
This estimation implies the use of a diameter growth model (e.g.Sánchez-González et al., 2005, 2006 or Tomé et al., 2006) to estimate diameter when cork will have 9 years.If the cork weight model includes cork thickness as regressor, a cork growth model (e.g.Almeida and Tomé, 2008;Almeida and Tomé, submitted, or Sánchez-González et al., 2008) is needed in order to estimate cork thickness at 9 years of age (ct9).
2. Estimate the cork back weight proportion at nine years of age (cbp9).Estimate cork biomass at 9 years free from cork back (wcm 9-b ): 3. Estimate cork biomass for t years of growth (wcmt): where wcm t is cork biomass with t years of growth, wcm t-b is cork biomass with t years free from cork back, ct t is cork thickness with t years of growth and cbp is cork back proportion.
The model developed by Paulo and Tomé (2010) allows, jointly with a dbh growth model (e.g.Tomé et al., 2006) and a cork growth model (e.g.Sanchez-González et al., 2008, or Almeida and Tomé, submitted), to predict the evolution of cork weight growth as is shown on Figure 2.
Cork growth modelling focuses on complete years of growth.Cork growth has been modelled with growth functions formulated as difference equations, using a methodology similar to the one used for dominant   height growth (Tomé, 2004;Sánchez-González et al., 2008, Fig. 3).Tomé (2004) and Almeida and Tomé (2008) also modelled the cork caliber using the cork thickness (complete years) as regressor and vice-versa.
In order to express tree cork growth rate, the authors proposed the use of a cork growth index defined as the cork thickness (complete years) of a cork with 9 years which corresponds to 8 complete growth years.Sánchez-González et al. (2007b) studied the variables influencing cork thickness at the end of the debarking period in corks with 9 years.As all the cork samples studied had the same age, cork caliber is equivalent to the cork growth index proposed by Tomé (2004) and Almeida and Tomé (submitted).This study was achieved by modelling cork caliber with mixed-models.The authors were not able to find covariates to explain cork caliber and proposed the simulation of mean cork caliber by calibration of the plot random component using a sample of 4 trees.This result expresses the high variability that can be found in cork growth index, most of which corresponds to inter-tree variability.This inter-tree variability was approached by Montes et al. (2005) through geoestatistical modelling, using ordinary kriging.
All the models described above refer to mature cork.There are just a few models for the virgin cork (the first cork produced by the first periderm).Tomé (2004) presented an equation to predict diameter over bark, essential to compute stand basal area that, for cork oak, is computed under bark.Sánchez-González et al. (2007c) developed a model to predict cork thickness at different heights using a taper equation and a linear relationship between virgin cork thickness at a given height and the corresponding diameter over cork.Both equations were simultaneously fitted using the seemingly unrelated estimation method.Tomé (2004) also presented an equation to predict virgin cork weight in juvenile trees.

Stone pine (Pinus pinea L.) nut
Stone pine nut is the main edible fruit collected in Mediterranean Forests, having been consumed by humans since ancient times (López-García, 1980).Cones are nowadays manual (climbing the trees) or mechanized (vibrating harvesters) collected mainly from natural stands, though an important effort is made to install and manage clonal grafted high productive plantations (Mutke et al., 2008).Since the 80's, large effort has been made on modelling stone pine (Pinus pinea L.) cone and nut production in Spain, Portugal and Italy.Concerning spatial variability, it has been mainly explained by tree size factors, as tree basal section or crown diameter (Cañadas, 2000;Calama and Montero, 2007;Freire, 2009), pointing to larger productions in larger trees; tree-level competition (Calama et al., 2008;Freire, 2009), with dominant trees reaching large cone crops; stocking attributes, indicating that larger tree fruit productions are attained in stands with low density (Cañadas, 2000;Calama et al., 2008;Freire, 2009); and site characteristics such as site index or soil attributes (Calama et al., 2008).While first efforts were carried on to develop stand level models (Castellani, 1989;García-Güemes, 1999), tree level models have performed better to model fruit production (Cañadas, 2000;Calama and Montero, 2007;Freire, 2009), since these models allow taking into account both between and within-tree variability in production.Nonnormality and heterocedasticity of data has been taken into account by using logarithmic transformation of the studied variable, robust regression, non-linear weighted regression together with mixed modelling techniques.As an example of a pure spatial model, Calama et al. (2008)   soils) or water presence due to a high water table.The model showed low bias and a modelling efficiency close to 40% in predicting pure between-tree variation.This model has been included as the cone production module in the general tree-level model for the multifunctional management of stone pine stands PINEA2 (Calama et al., 2007).Pure spatial models approached masting habit by def ining as response variable the average fruit production for a series of years (five years usually).In other studies masting has been modelled as a function of climate factors, showing that rainfall and temperature in the previous spring and autumn affect female flower induction in stone pine (Mutke et al., 2005;Freire 2009;Calama et al., 2010), indicating that a dependence on meteorological phenomena occurred up to four years before cone maturation.A three-year delayed negative influence on cone production was detected in Pinus pinea crops (Mutke et al., 2005;Calama et al., 2010) associated with resource depletion and exhaustion of reserves after a bumper crop.Based on this result and using a regional data series of 43 years, Mutke et al. (2005) developed a model for predicting the annual average production of cones (Wc, kg ha -1 ) for a given year in the Spanish Northern Plateau: Wc = 9.740 -0.487 log (Wc -3 ) + 0.0025 P 5,-4 + + 0.0076 P oct-4 + 0.0053 P 5, -3 -0.277 T jj-3 + [2] + 0.0024 P -1 where Wc -3 is the average regional production three years before, P 5,-4 and P 5,--3 are the rainfall from January to May four and three years before; P oct-4 is the rainfall for October four years before; T jj-3 is the average temperature for June and July three years before; P -1 indicates the annual rainfall between September and August before maturation.These climate events are in close temporal relation with important phenological processes such as floral bud initiation and differentiation, female flower emergence, post-pollination floral survival and cone ripening.Figure 4 shows the accuracy of the regional model in predicting interannual variability at regional scale.
Based on the previous works, Calama et al. (2010) have developed a model for predicting annual cone crop at tree level in the Spanish Northern Plateau.The main problem of temporal predictions at tree level is related with the abundance of zeroes in fruiting events, which can reach more than 50% of the observations, and prevent modellers to use classical statistical approaches.To overcome this problem, the proposed model was formulated as a zero inflated one, including a logistic model for defining the probability of presence/absence of cones together with a log-normal function for modelling the intensity of fruiting process (weight of cones) conditional to the occurrence.Hence, fruiting was considered as a two-phase process, one ruling the presence of cones, mainly related with tree maturity and climate events occurring before female floral emergence, and other controlling the abundance, which is also affected by post-floral conditions.The proposed model explained more than 35%, 50% and 75% of expected interannual variability at tree, plot and regional level.In that sense, the model is a largely useful tool for early detection of good and bad crop years, helping forest managers and nut industries in planning strategies and decision making.
The abovementioned models have focused on predicting the weight of cones, while the main commercial product is the pine nut without shell (white nuts), which is traditionally marketed in terms of weight.Nut industry major requirements are associated with cones attaining larger nut yield (defined as the ratio between nuts weight and cone weight) and larger nuts.Little effect of silvicultural practices or site attributes on those variables have been detected (Montero et al., 2004), while both variables show large interannual variability (Calama and Montero, 2007) and are closely related with the average size and weight of the single cone.Subsequently, attributes defining cone content quality are mainly ruled by the weather conditions during the last year of cone maturation (Calama et al., 2007).Based on this, different models which relate nut yield and single nut weight with the weight of cones (Fig. 5  . have been developed (Calama et al., 2007;Morales, 2009).These models have a great potential to support the nut industry, since they allow them to make a prior classification of cones depending on its expected nut yield and quality.

Mushrooms
Wild edible and medicinal mushrooms represent an important non-wood forest product world-wide (Boa, 2004) to the extent that the commercial value of forest fungi may equal or even surpass the value of timber (Bonet et al., 2008;Alexander et al., 2002;Arnolds, 1995; Oria de Rueda and Martínez de Azagra, 1991).In addition, mushroom picking has become one of the most important forest recreational activities in many parts of Europe (Lund et al., 1998;Mogas et al., 2005;Martínez de Aragón, 2005;Martínez de Aragón et al., 2007).Consequently, there is growing interest on the part of forest owners and managers to inventory, predict, and manage forests to improve the production of marketed mushrooms (Pilz and Molina, 2002).The inclusion of mushroom yields as an explicit objective in forest management and planning requires models for assessing, in a quantitative way, the production of mushrooms in different forest stands and management schedules.Such models can be developed with statistical methods, using empirical data on mushroom production and forest stand characteristics (Bonet et al., 2008(Bonet et al., , 2010)).Once predictive models are developed, they can be included in forest simulators to provide quantitative information on mushroom production and its economic effects in alternative forest management schedules (Palahí et al., 2009).
Site and growing stock characteristics are the most reasonable predictors when developing an empirical mushroom production model for forest planning because the site and stand characteristics are known factors within forest simulators.In addition, stand characteristics can be modified through forest management.However, mushroom productions also depend on weather conditions such as timing and quantity of rainfall, which are not equally useful in forest planning because they cannot be accurately predicted beyond a few weeks.
Since annual mushroom production varies considerably in response to weather conditions, the construction of reliable models for predicting mushroom yields requires collecting large quantities of empirical data over several years because there are multiple factors responsible for high temporal variation in mushroom productions: precipitation, temperature, frost, evapotranspiration, relative humidity, and water def icits (Wilkins and Harris, 1946;Ohenoja, 1993;O'Dell et al., 1999;Straatsma et al., 2001;Martínez de Aragón et al., 2007).
Empirical models for predicting the yield of wild mushrooms are presented in Bonet et al. (2008) to predict the total, edible, marketed mushroom yield in P. sylvestris plantations in North-east of Spain and in Bonet et al. (2010) to predict the mushroom yield in three species of pines (including plantations and natural forests) in North-East Spain.The model of Bonet et al. (2008) predicting edible mushroom yield was as follows: ln(y ij ) = 0.981 + 2.483ln(G) -0.128G + [3] + 0.934cos(Asp) -0.0135Slo 1.5   where y ij is the mushroom production of plot i in year j, G is stand basal area (m 2 ha -1 ), Asp is aspect (radians), Slo is slope (%, i.e. 45 degrees is equal to 100%).
According to such model the highest production occurs when stand basal area is within the range 10-20 m 2 ha -1 .Aspect and slope are other factors which strongly affect the predicted production.Northern aspects have the highest productions and southern ones the lowest while increasing slopes decreases the total mushroom production (Fig. 6).Bonet et al. (2008) proposed that such pattern may be explained by the fact that Catalonian pine stands also reach their peak in basal area growth rate at the basal area of 10-20 m 2 ha -1 .This relationship between peak mushroom production and peak stand growth rate is in accordance with the work of Kuikka et al. (2003) who demonstrated that formation of mycorrhizal sporocarps was strongly correlated with the growth and photosynthetic rate of the host trees.
In the model by Bonet et al. (2010) elevation also strongly affects the predicted mushroom production.When elevation increases from the lowest (500 m) to the highest level (1,500 m) mushroom productivity increases by more than ten times.The sharp increase in productivity for higher elevations and for north-facing forests is likely to be related to the availability of water.Water is the most limiting growth factor in the forests of southern Central Pyrenees where light is abundant and temperatures are mild all year round.
In addition, the model by Bonet et al. (2010) predicts that the stand dominated by P. sylvestris produce higher (more than double) mushrooms yields than the ones dominated by P. nigra or P. halepensis.In addition, mushroom yields were greater in naturally regenerated P. sylvestris stands than planted P. sylvestris stands (Fig. 7).However, such results require more research to better separate the site effects from the tree species (and plantation versus natural) effects.To do this, more balanced data are needed in which different species grow in similar sites.
In Bonet et al. (2008Bonet et al. ( , 2010)), the modelling data were characterised by multiple measurements for each individual sampling unit (several observations for the same plot) and several sampling units (plots) measured in the same year.Therefore, a multilevel linear mixed model approach with both fixed and random components was used (Fox et al., 2001).The random factor associated with the measurement year was best ex-  .Total mushroom production as a function of stand basal area, dominant tree species and, for P. sylvestris differentiated by natural or plantation-established stands.Elevation was 1,000 m, slope was 20% and aspect was east (Bonet et al., 2008).
plained by total autumn rainfall.Figure 8 shows the effect of total autumn rainfall on the production of mushrooms when the mushroom yield model is used together with the random year factor model (which has as independent variable total autumn rainfall).The models predict five times greater mushrooms yields when autumn rainfall is doubled (Fig. 8).
The study by Palahí et al. (2009) demonstrated how the mushroom yield models developed by Bonet et al. (2010) could be linked to a forest simulation-optimization system in order to optimise the joint production of timber and mushrooms.The results of such study showed that mushrooms produce more prof it than timber in most P. sylvestris and P. nigra stands of central Catalonia (north-east of Spain).In P. sylvestris stands with good site conditions for mushrooms the soil expectation value of mushroom harvesting was often 10 times higher than the soil expectation value of timber production.

Berries
Berry picking has become one of the most important forest recreational activities in many parts of Northern Europe (Saastamoinen et al., 2000;Ihalainen and Pukkala, 2001).Main berries species actually collected in Northern countries are bilberry (Vaccinium myrtillus) and cowberry (Vaccinium vitis-idaea).Similarly to nuts and mushrooms, berries show a large inter-annual variability in response to weather conditions, especially frost and precipitation events defining pollination success together with previous year crop.In Finland, bilberry produces the most abundant yields in medium shadow on medium and rather poor sites; the crown density of 10-50% is the most favourable for a good bilberry yield (Raatikainen et al., 1984).
The lack of empirical measurements on berry yields has been overcome with expert modelling so that experienced berry pickers and foresters have rated different stands according to their berry production (Ihalainen and Pukkala, 2001;Ihalainen et al., 2002Ihalainen et al., , 2005)).Ihalainen et al. (2002) employed photographs of different stands, whereas Ihalainen and Pukkala (2001) estimated visually the bilberry productivity during the field inventory of forest stands.Bilberry yields based on expertise have been modelled in an ordinal scale by using site and stand characteristics as predictors.Furthermore, Turtiainen et al. (2005) have calibrated the expert models of Ihalainen et al. (2005) with a set of measured yield data from various sources.The expert models have showed their accuracy in defining sites which are worthy of berry collection, as well as estimating national and regional yields of bilberry and cowberry in Finland.
The model by Ihalainen et al. (2003) was the first in using field measurements from permanent plots in order to estimate bilberry and cowberry yields as a function of site and stand attributes which are known in forest planning.According to the fixed part of the mixed linear model for bilberry yield, the highest productions are found in mature stands growing on medium fertility sites: ln(y + 1) = 0.0830 + 0.0103 t + 0.9904 D1 + [4] + 0.4997 D2 where y refers to bilberry yield (kg/ha); t is stand age (years); D1 is a dummy variable equalling 1 if forest site is medium, and zero otherwise; D2 is a dummy variable equalling 1 if forest site is rather poor and zero otherwise.When the logarithmic predictions for bilberry yields are back-transformed, the multiplier 2.4507 should be used for bias correction.The first empirical model identified a large amount of the variability associated with non-considered factors acting at stand and year level.In addition, the annual variation in bilberry yields was not included in the previous model.Therefore, Miina et al. (2009) proposed a new empirical bilberry yield model that first predict the coverage of bilberry as a function of e.g.dominant tree species, site fertility, regeneration method, stand age and stand basal area.Then, the annual yield .Total mushroom production as a function of total autumn rainfall in a north and south aspect of natural P. sylvestris stands using the models by Bonet et al. (2010).
Stand basal area was 15 m 2 ha -1 , elevation was 1,000 m and slope was 20%.
of bilberry is predicted as a function of bilberry coverage, dominant tree species and stand basal area.The bilberry coverage was predicted using a multilevel binomial mixed model: where UnripeBerries is the annual number of unripe berries per m 2 .No weather covariates were included in the models, but interannual variability in observed number of unripe berries was modelled by normally distributed random year effects u with zero mean and variance equal to 0.2907 for pine stands and 0.4504 for spruce stands.In stochastic simulations, the bilberry yield is predicted several times for each year by drawing the random year effects (u) from the normal distribution, and the annual bilberry yield is computed as the mean of the outcomes.
In the simulations of Miina et al. (2009Miina et al. ( , 2010)), 80% of the unripe berries were assumed to become ripe.The number of ripe berries was converted into the bilberry yield (kg ha -1 ) by multiplying it by the mean fresh weight (0.35 g) of one ripe berry.The recently developed empirical model of Miina et al. (2009) considers also the effect of stand density -and consequently thinnings-on the bilberry yields.In Finland, optimization of timber production has been used to evaluate the silvicultural recommendations for forestry practice.However, due to the increasing importance of nonwood forest products (including wild berries), there was a need for optimizing the joint production of timber and bilberries (Miina et al., 2010).The empirical bilberry yield model by Miina et al. (2009) was included in a stand growth simulator for forest planning purposes (Fig. 9).The results showed that the effect of bilberry production on the optimal stand management  increased with increasing bilberry price (Miina et al., 2010).In a Scots pine stand on mesic heath site, bilberry production affected optimal stand management already with a price of 2 €/kg, which corresponds to the average price of bilberries offered for sale in Finland (MARSI, 2008).Compared to timber production, joint production led to longer rotation lengths, higher thinning intensities, more frequent thinnings, and higher share of Scots pine in a mixed stand.In the pine stand the soil expectation value (SEV) of bilberry production, calculated with 3% discounting rate, exceeded the SEV of timber production when bilberry price was 4 €/kg.Because statistical models do not take into account all the variation observed in berry yields, the model predictions are typically around the mean yield.This means that the yields on sites where berries are usually collected can be assumed to be much higher than the predictions calculated by the models.Thus, the effect of bilberry production on the optimal management of stands that are known to be good berry stands may be stronger than the results by Miina et al. (2010) suggested.

Resin
Resin tapping has been an important traditional forest activity in Northern countries (Pisarenko and Strakhov, 1996) as well as in the Mediterranean basin.Russia is still one of the main producers of resin in the world, with Pinus sylvestris being the main species tapped (FAO, 1995), while on Finland there exist an increasing demand for Norway spruce (Picea abies) resin, due to its medical properties (Jokiaho, 2010).However, little information is available about the management of this resource in these northern countries.In Mediterranean countries resin tapping was a main rural activity until the 1970s when the international crisis in natural resin prices rendered this traditional labor no longer profitable.During that period, the resin from Pinus pinaster (and other Pinus species as Pinus halepensis) was one of the most important non timber products in some Mediterranean countries, namely Spain, Portugal, France and Greece.Recently, an increased demand for natural resins pushed up the prices, with a few abandoned stands being tapped again (Nanos et al., 2000).Presently Pinus pinaster Ait (maritime pine) is the species mainly tapped in Spain and Portugal, while Pinus halepensis is tapped for resin in some areas of Greece.Unlike other NWFP, resin production shows a large incompatibility with high quality timber production and tree survival.As in the case of cork extraction, management decisions concerning the number of opened faces or the length of annual face opening have a large influence on final resin production.
Genetic control and site conditions have been identified as the factors largely affecting spatial variability in Pinus pinaster resin production (Nanos et al., 2000).Large variability has been detected at both within and between stand for the Northern Plateau of Spain.Within stand variability was approached by modelling resin distribution at stand level (Nanos et al., 2000) by using the Weibull and the Chaudhry-Ahmad probability functions; while between stand variability was approached using geostatistical modelling followed by ordinary kriging (Nanos et al., 2001) which allows one to detect areas with larger resin production (Fig. 10).
Concerning empirical predictive models for resin production, a recent work by Spanos et al. (2009) has modelled annual tree-level resin yield production on Aleppo pine (Pinus halepensis), recognizing that, in spite of a large variability probably due to genetics, breast height diameter is the main factor explaining between trees variability (Fig. 11), a result also stated by Jokiaho (2010) in relation with Norway spruce resin.These authors have identified several other influence factors, as site, social position or crown attributes, together with management decisions, as is the number of opened faces, as main potential additional explanatory covariates.

Trends and future research in modelling NWFP
Besides timber production, forest management is increasingly focused on multiple objective uses of the forests.Taking into account the challenges and limitations shown in Table 2, the trend in modeling NWFP seems to go to multivariate hybrid models formulated taking into account the hierarchical structure and the large temporal variability inherent to the used data.These models should combine empirical measures with climatic and soil information to explain the spatial and temporal variability of NWFP production and to understand related physiological processes.Among the first approaches we can mention, e.g., the one carried by Calama et al. (2010) when modeling the weight of sound cones cropped from a tree in a given year.Subsequent steps should consider the complexity of the involved processes by considering carbohydrate allocation to non wood production in process-based models.In this sense, both non-linear processes and biological switches (Thornley and Johnson, 1990) should be considered when modeling the NWFP.
Data scarceness has been pointed as one of the main factors limiting models for NWFP.Modelling based on expert evaluations could be useful when forest use such like scenic beauty, recreational value, habitat suitability for wildlife and biodiversity is difficult or impossible to measure.In addition, expert modelling could be used to substitute too expensive or scarce empirical measurements.In the method, knowledge on the relationships between stand characteristics and the non-wood forest products are elicited and converted into an expert model (Pukkala, 2002b).Several sources of uncertainty (e.g.inconsistencies in opinions, varying levels of expertise) characterize expert modelling but the method is much quicker and cheaper than empirical modelling.However, Kangas (1998) has concluded that expert modelling should be considered only as a temporary solution until sufficient data for empirical models are available.
In modelling NWFPs, it would be possible to combine different modelling approaches.First, expert evaluation could help to find out the main factors affecting the production of NWFPs and to direct the field measurements to the most crucial gaps of knowledge.Finally, field measurements could be used to fit the parameters of the production function.As an example concerning mushroom modelling task, mushroom pickers, who offer mushrooms for sale, would be asked to evaluate which kinds of forests are good for mushroom picking.The imaginary forest stands would be used in evaluation.Field measurements are needed especially to determine the effect of stand density and thinnings on mushroom yields.If it is not possible to collect the yield data for several years, the annual variation in mushroom yields could be based on the statistics on the amounts of mushrooms offered for sale.Time series analyses (ARIMA modelling) can be used to determine the between-year variation in mushroom yields.
Regarding cork yield modeling, the main challenge will be to identify environmental factors explaining the high variability among trees that is found in cork thickness.This variability affects not only cork thickness modeling but also cork weight models.In this sense Sánchez-González et al. (2007b) suggest the use of soil isotopic signatures as indicators of cork thickness variation to be used in hybrid models.Related with this, identification of those possible variables explaining large patterns of non explained (by the actual models) variability between trees detected in cone and resin productions should be one of the main objectives in modeling research.
Finally, most of the existing models for NWFP focus on predicting the total amount of product, while quality is one of the main traits defining both end-use and price.Characteristics such as cork caliper, nut weight, chemical contents of resin, and berry size, colour or taste should be considered as sound response variables in future models for NWFP products. of resin production (kg/tree) and DBH (cm) for Aleppo pine in Greece (Spanos et al., 2009).

Figure 1 .
Figure 1.Schematic representation of a cork sample with 9 years, showing the 8 complete cork rings and the two half rings (adapted from Natividade, 1940).

Figure 2 .
Figure 2.Estimated values of cork biomass from 1 to 14 years of age, for the two pairs of hypothetical trees.Tree A (dash lines): diameter at breast height of 35 cm, vertical debarked height of 1.5 m and number of debarked branches equal to 1 (tree only debarked in the stem); tree B: (solid lines): diameter at breast height of 65 cm, vertical debarked height of 4.5 m and number of debarked branches equal to 3. Cork thickness after boiling for trees A1 and B1 (graphic lines identified with «o») was assumed as 25 mm (2 nd cork thickness class); for trees A2 and B2 (graphic lines identified with «x») this variable was assumed as 50 mm (5 th cork thickness class) (from Paulo and Tomé, 2010).

Figure 3 .
Figure 3. Cork growth curves for accumulated cork thickness in complete years of 15.75, 22.5, 29.25, 36 and 42.75 mm at 9 years for the selected model overlaid on the trajectories of observed values over time (Sánchez-González et al., 2008).

Figure 5 .
Figure 5. Nut yield (a) and single nut weight (b) as a function of the dry weight of the cone(Morales, 2009).

Figure 6 .
Figure6.Dependence of the total mushroom production (kg ha -1 ) on stand basal area, slope (%) and aspect (cosAspect).The dots indicate the means of the measured productions of the plots (3-year averages) and the curves are the respective predictions with Equation 3(Bonet et al., 2008).
Figure7.Total mushroom production as a function of stand basal area, dominant tree species and, for P. sylvestris differentiated by natural or plantation-established stands.Elevation was 1,000 m, slope was 20% and aspect was east(Bonet et al., 2008).
Figure8.Total mushroom production as a function of total autumn rainfall in a north and south aspect of natural P. sylvestris stands using the models byBonet et al. (2010).Stand basal area was 15 m 2 ha -1 , elevation was 1,000 m and slope was 20%.

Figure 9 .
Figure 9. Development of bilberry coverage (a) and yield (b) along the simulated development of Scots pine stand on mesic heath site.The management schedule was optimized by maximizing the soil expectation value (3%) of the joint production of timber and bilberries(Miina et al., 2010).

Table 2 .
Main challenges in modelling NWFP