Kaj Rosén

Hervé Jactel

Manuela Branco

Wojciech Grodzki

Jean-Sébastien Jacquet

Francisco Moreira

Sigrid Netherer

Mart-Jan Schelhaas

Z Tome

Report on specific risk analysis in regional forests of Europe under various FMAs

Keywords: Tree mortality, risk model, logistic regression, Eucalyptus, forest management

This report is a deliverable from the EU FP6 Integrated Project EFORWOOD -Tools for Sustainability Impact Assessment of the Forestry-Wood Chain. The main objective of EFORWOOD was to develop a tool for Sustainability Impact Assessment (SIA) of Forestry-Wood Chains (FWC) at various scales of geographic area and time perspective. A FWC is determined by economic, ecological, technical, political and social factors, and consists of a number of interconnected processes, from forest regeneration to the end-of-life scenarios of wood-based products. EFORWOOD produced

Introduction

During the 20 th century forest disturbances have been increasing in Europe [START_REF] Hanewinkel | Seventy-seven years of natural disturbances in a mountain forest area-the influence of storm, snow, and insect damage analysed with a longterm time series[END_REF]2003) and a similar trend has been observed in other regions [START_REF] Lutz | Twentieth-century decline of largediameter trees in Yosemite National Park[END_REF]). Consequently, there is an increase concern with forest damage, in particular regarding tree mortality. Forest disturbances are due to both abiotic agents (storms, fire and snow) and biotic, pests and diseases (Brasier and Scott 2008;[START_REF] Hanewinkel | Seventy-seven years of natural disturbances in a mountain forest area-the influence of storm, snow, and insect damage analysed with a longterm time series[END_REF][START_REF] Peltola | Model computations on the critical combination of snow loading and windspeed for snow damage of Scots pine, Norway spruce and Birch sp. at stand edge[END_REF][START_REF] Schelhaas | Natural disturbances in the European forests in the 19th and 20th centuries[END_REF][START_REF] Valinger | Models to assess the risk of snow and wind damage in Pine, Spruce and Birch forests in Sweden[END_REF]. Climate change and the increasing number of invasive forest pests and diseases have been pointed as possible causes of increasing forest disturbances (Brasier and Scott 2008;[START_REF] Holmes | P. 381-406 in The economics of forest disturbances: wildfires[END_REF][START_REF] Logan | Assessing the impacts of global warming on forest pest dynamics[END_REF]. However, changes in forest management and consequent forest structure, such as volume of growing stock, and stand age, have been suggested as the most likely explanations for an increase in damage from disturbances [START_REF] Schelhaas | Natural disturbances in the European forests in the 19th and 20th centuries[END_REF]. A revision work conducted by [START_REF] Jactel | The influences of forest stand management on biotic and abiotic risks[END_REF] it evidenced how distinct forest management strategies may impact forest susceptibility to both biotic and abiotic risks in multiple ways. Evaluating the risk of disturbance in regard to forest management is thus crucial to commit forest management towards sustainability. Since forests provide economic services as well as several ecological services, such as carbon sequestration, water resources, climate amenity and biodiversity, forest tree mortality is a public concern as much as a concern for the forest based sector.

This study is integrated in the European 4-year research project EFORWOOD "Sustainability impact assessment of the Forestry-wood chain". The principal aim of the EFORWOOD project is to develop a decision support tool for evaluating and developing the European forest based sector to sustainable development, from economic, environmental and social perspectives. One aim of the EFORWOOD project is to simulate the effect of forest management alternatives on selected sustainability indicators. The present study focuses on the Iberian Peninsula case, particularly regarding eucalyptus forests, one of the most economically important plantation forests in this region. Tree mortality is here selected as a damage indicator. The aim is to obtain a mortality risk model which will allow predicting tree mortality for distinct stand attributes and site conditions. Once developed, such mortality risk model will be included in the SIMPLOT regional simulator (Barreiro and Tomé), a regional simulator for Portuguese forest cases. Since both site selection and stand attributes may be changed by forest operations processes these variables may be subsequently linked to forest management strategies.

Eucalyptus species were first introduced in Europe during the 19 th century [START_REF] Pimentel | Eucalyptus globulus. Modo de vegetar, cultura, produção, etc. Typographia universal[END_REF], mostly for ornamental purposes [START_REF] Potts | Natural disturbances in the European forests in the 19th and 20th centuries[END_REF]). During the 20 th century eucalypts become extensively cultivated in the Iberian Peninsula, mainly for pulp. At present E. globulus is one of the main forest species cultivated in Portugal comprising a planted area of 717 432 ha, which represents about 1/3 of the cultivated forest in the country, with a stand volume of 38 million m 3 [START_REF] Tomé | Inventário florestal 2005-2006. Áreas, volumes e biomassas dos povoamentos florestais. Resultados Nacionais e por NUT's II e III[END_REF]). Most of the E. globulus plantations in Portugal are located near the littoral (Fig. 1).

For more than one century, eucalyptus plantations in Europe were extremely healthy forest ecosystems, exhibiting very low natural mortality as observed in other regions [START_REF] Paine | Insect pests of eucalypts in California: implications of managing invasive species[END_REF][START_REF] Wingfield | The future of exotic plantation forestry in the tropics and southern hemisphere: lessons from pitch canker[END_REF]. A pattern which is mostly attributed to its "enemy free" status and is common to other exotic plantation forests as evidenced with Pinus radiata plantations in New Zealand (Gadgil and Bain 1999;[START_REF] Paine | Insect pests of eucalypts in California: implications of managing invasive species[END_REF][START_REF] Wingfield | The future of exotic plantation forestry in the tropics and southern hemisphere: lessons from pitch canker[END_REF]. Within decades an increasing number of pests and diseases have changed this status. Presently, nine important eucalyptus pests, including two wood borers, Phoracantha semipunctata and P. recurva (Valente et al. 2008), three sap-sucking psylids, Ctenraytaina spatulata, C. eucaliptii and Glycaspis brimblecombei [START_REF] Valente | First record of the Red Gum Lerp Psyllid, Glycaspis brimblecombei Moore (Hemiptera: Psyllidae), in Europe[END_REF][START_REF] Valente | First record of the Australian psyllid Ctenarytaina spatulata Taylor (Homoptera: Psyllidae) in Europe[END_REF]), one eriophyid mite, Rhombacus eucalypti (Ferreira et al. 2006), two eulophid gall wasps [START_REF] Jactel | The influences of forest stand management on biotic and abiotic risks[END_REF][START_REF] Villar | Dos especies australianas de eulófidos, muy dañinas para Eucalyptus spp., introducidas en el nordeste ibérico (Hymenoptera: Eulophidae)[END_REF], one defoliator, Gonipterus scutellatus [START_REF] Mansilla | Presencia sobre Eucalyptus globulus Labill de Gonipterus scutellatus Gyll[END_REF], as well as pathogens, including several fungi of Mycosphaerella leaf blotch disease [START_REF] Silva | Species identification of Mycosphaerella leaf blotch disease (MLD) on Eucalyptus globulus in Portugal, Universidade de Évora[END_REF], Botryosphaeria sp. trunk canker disease and a root fungi (Valente et al. 2008) affect the health of eucalyptus plantations in the Iberian Peninsula, compromising its productivity and sustainability. In this context it is most relevant to understand how site selection and forest management may influence tree mortality as well as other types of damage.

In this study we intend to develop a tree mortality risk model for eucalyptus plantations, by using a regression modelling approach but without focusing or specifying any particular agent. Regression methods have become an integral component of data modelling describing the relationship between a response variable and one or more explanatory variables. In this work we use a logistic regression model which is a standard method of analysis whenever the response variable is binary or percentage such as the present case [START_REF] Hosmer | Applied logistic regression, 2nd Edition[END_REF].

The aim of this paper is to explore the structural and environmental factors influencing mortality of Eucalyptus caused by biotic and climatic factors. Both windstorms and fire mortality are excluded from the present analysis. Windstorms are rare in the studied area and were not recorded in the NFI data used for this analysis. On the contrary, fire events are very important and a well know mortality agent in the studied area. Due to its uniqueness, stands where fire events had occurred were excluded from the analysis and will be treated in a separate analysis. Our hypothesis is that both site variables and stand management are important factors in mortality risk, by affecting both directly and indirectly the damage incidence from these agents, such as through its influence upon plant physiological status (e.g. [START_REF] Hanks | Water relations of host trees and resistance to the phloem-boring beetle Phoracantha semipunctata F. (Coleoptera: Cerambycidae)[END_REF]. In particular, we are interested in the role of both site selection and stand attributes, as these may be controlled by stand management strategies.

Material and methods

Forest damage data

Data regarding stand structure and Eucalyptus globulus tree mortality came from different sources: the two last National Forest Inventories (NFI's) and permanent plots that were established and had been annually remeasured between 1984 and 2004 by SILVICAIMA and Celbi pulp companies.

The fitting data used is from the National Forest Inventory (NFI) of 1995-98 with field data collected during 1997 and 1998 and NFI 2005-06 with field work accomplished during 2005 and 2006. All the plots of the two NFI's where E. globulus was present were considered, a total of 567 plots from the NFI 1995-98 and 1115 plots from the NFI 2005-06. This dataset also includes one observation per plot and cutting cycle randomly chosen from each of the 7526 permanent plots, thus avoiding autocorrelation. The remeasurements of the permanent plots include more than one rotation in 0.29% of the permanent plots. The studied plots are located throughout the area of the distribution of the species in the country (Figure 1).

The data from NFI 1995-98 were gathered in concentric circular sample plots with an area of 250 m 2 , 500 m 2 and 1000 m 2 respectively in which trees were selected according to their breast height diameter ( larger than 7.5 cm, larger than 17.5 cm and larger than 27.7 cm, respectively) (Florestais 1999).

In the NFI of 2005-06 all trees inside a circular plot with 500 m 2 with a diameter at breast height larger than 5 cm were measured (Ferreira et al. 2005).

The shape and area of the permanent plots are not equal. Celbi permanent plots are rectangular and square with an area varies between 732 m 2 and 1200m 2 where all trees were considered. In Silvicaima permanent plots all trees inside a rectangular, square or circular area that goes from 100m 2 to 757m 2 were measured.

Explanatory structural and environmental variables

Stand description in these plots includes management regime (coppiced or high forest), composition (pure or mixed stands), stand structure (even-aged or uneven-aged) and stand age (only in even aged stands). Dominant height (hdom, m) was calculated by taking the mean of the 100 largest breast height diameter E. globulus trees per hectare. Dominant height was only used for the even-aged stands. All the trees in the sample plot, including trees of other species, were considered in the calculation of basal area (G, m 2 ha -1 ) and number of trees per hectare (N). The dead standing stems of E. globulus present in the plots during the field survey were counted and expressed as a proportion of the total number of E. globulus dead or alive trees in the plot (Pdead). All the predictor variables used for building the model are shown in table 1.

A dataset with 7510 plots coming from the remeasurements of the permanent plots, independent of the fitting dataset, was used for model validation. Only pure and even-aged plots were available. In table 2 is shown the characterizing of the fitting and validation datasets. standard deviation

Model building

The risk model was developed using a two-step approach [START_REF] Woollons | Even-aged stand mortality estimation through a two-step regression process[END_REF]. The probability of occurrence of at least one dead tree per hectare was first constructed using the whole dataset.

If we denote the event of dead occurrence as D, let X 1 , X 2 ,…,X k be the set of explanatory variables and let p be the probability of D=1 , eq. ( 1).

[1]

A mortality equation was then constructed, also by a logistic regression, to predict the proportion of dead stems per hectare, using just the data from the plots in which damage had been observed. The model can be written with the same basic forms as eq.

[1], where p is now the proportion of dead stems per hectare (Pdead) and the set of explanatory variables can eventually differ.

Both models were built using a stepwise selection (combining the elements of forward selection and backward elimination) of the regressors among those included in table 1. The significance level of the likelihood ratio test to allow a variable into the model was 0.05 and the same significance level was required for a variable to stay in the model. The multivariable logistic regression analysis was performed using the logistic procedure of SAS 9.1 statistical software (SAS 2004a).

For both models odds ratio were obtained for explanatory variables thought the eq. [2] where the β i 's are the estimated coefficients for the models. So the odds ratio estimate the so called relative risk, i.e. the variation in the damage risk (Dead and Pdead) by a particular factor or change in the value of a covariate. It was necessary to adjust the scale of some of the independent variables for which the change in odds ratio per unit was small, to ensure the appropriate scale for that variable [START_REF] Hosmer | Applied logistic regression, 2nd Edition[END_REF].

= Some nonparametric measures of association, that are helpful in assessing the validity of the fitted model, were also obtained with the same procedure. The outcome provides some basic calculation about trends among the data values (concordant, discordant and tied) and four rank correlation indexes estimated from the preceding ones (Somers'D, Gamma, Tau-a and C) (SAS 2004b). For all pairs of observations with different values on the response variable, a pair is concordant (c) if the observation with the lower ordered response value also has a lower predicted mean score than the observation with the higher ordered response value. A pair of observations is discordant (d) if the observation with the lower ordered response value has a higher predicted mean score than the observation with the higher ordered response value. Pairs that are neither concordant nor discordant are referred to as tied. Somer's D (eq.

[3]) is used to determine the strength and direction of relation between concordant and discordant variables. Its values range from -1.0 (all pairs disagree) to 1.0 (all pairs agree). The Goodman-Kruskal Gamma method ignores tied pairs (eq. [4]). Its values range from -1.0 (no association) to 1.0 (perfect association). Kendall's Tau-a is a modification of Somer's D (eq.

[5]). It ranges from 0 (no association) to 1 (perfect association). Tau-a = (c-d)/0.5n(n-1)

Model Validation

For each plot in the validation dataset the proportion of dead trees estimated by the second logistic model (Pdead) was reduced by the probability of dead occurrence estimated by the first logistic model (probability of D=1). This value was taken as the predictor of the expected value of the mortality in that plot. Bias of the predictors was evaluated by the mean of the prediction residuals (observed minus predicted mortality) and precision as 1 minus the mean of the absolute values of the prediction residuals.

In the validation dataset, for the probability of dead occurrence equation (Dead), the observations were divided into 10 equal sized groups according to their predicted probabilities and the relationship between the observed and expected numbers of dead occurrence (D=1) and dead absence (D=0) within each group were analyzed by a scatters-plot.

A contingency table were also constructed to analyze the relationship between the observed and expected numbers of dead occurrence (D=1) and dead absence (D=0). Sensitivity, as known as true positive (the model capacity to predict a event (D=1) correctly), false positive (proportion of predicted events that were observed as non events (D=0)), specificity, as known as true negative (the model capacity to predict a non event correctly), false negative (proportion of predicted non events that were observed as events), and correct percentage (percentage of events and non events classify correctly) were also calculated.

Results

Model selection

On average tree mortality rate was observed to be low for the Eucalyptus globulus stands (11.9%, 12.9%, 4.7% and 6.6% in Silvicaima, Celbi, NFI 2005-06 andNFI 1995-98, respectively), however it was a great variability, values ranging from 0 up to 100%.

The logistic regression procedure selected two climatic variables and seven attributes of the stand to predict the probability of occurrence of at least one dead tree per hectare (Dead) (table 3). All the significant explanatory variables, according to the Wald test, were significant with a p-value < 0.0001.

If interpreted odds ratio as relative risk the probability of dead occurrence in high forests would be 2 times more than in coppice cases, 9.6 times more in uneven aged than in even aged stands and just 0.11 less in mixed stands than in pure stands. In the same way, increases of every 2 m in dominant height increase the probability of mortality occurrence 1.14 times and increases of every 2°C in the mean temperature of the 3 hottest months leads to an increase of 2.24 times. Otherwise, increases in mean annual evaporation of every 100 mm decrease almost 0.7 times the probability of dead occurrence.

For the variable Pdead, the stepwise regression selected all the variables selected for the Dead model plus two new variables: the linear effect of mean annual precipitation and the stand density (table 3). The level of significance of all variables was p < 0.0001.

The interpretation of odds ratio in this model shows that Pdead decreases 0.99 times with the increase of each 10 trees in the stand density, 0.63 times with the increase of each 100 mm in evaporation and 0.94 times with the increase of each 2 m in the dominant height. In opposite, increases of each 100 mm in mean annual precipitation and 10 m 2 ha -1 in stand basal area also increase the proportion of dead 1.3 times and 1.2 times, respectively. Equally increases of each 2 years in the stand age and 2 ºC in mean temperature of the 3 hottest months, increases the proportion of dead 1.5 times and 1.4 times respectively. The proportion of dead in uneven aged stands would be 3.7 times more than in even aged stands, 0.86 times less in coppice forest than in high forest and 0.6 times less in mixed stands than in pure stands (table 3).

In table 4 some measures of association of predicted probabilities and observed responses for Dead and Pdead show that both models have more pairs of wells concordant, 79.9% and 64%, respectively, than discordant, 19.9% and 35.2% respectively and very few ties (0.2 and 0.8 respectively). The results of the four rank correlation indexes estimation shows that beside probability of dead model has better results the proportion of dead model has good results to. Somer's D's values shows that there is a large number of pares of observations that agree and the Gamma and c values indicates that there is a high association between the observed and the predicted pairs ranking (table 4). Figures 2 to 7 allow a visualization of how some of the selected variables influence the mortality risk as predicted by the models.

Model validation

Precision and bias mean values, near 0, show that both models have a good fit (

Discussion

A logistic regression approach was used to model the tree mortality risk of Eucalyptus globulus stands in Portugal. The discrete nature of the dependent variable well matched the use of logistic regression, a generally used and accepted method of analysis of binary outcome variables, widely applied in forestry to estimate tree and stand survival (e.g. [START_REF] Monserud | Modeling individual tree mortality for Austrian forest species[END_REF][START_REF] Yao | A generalized logistic model of individual tree mortality for aspen, white spruce, and lodgepole pine in Alberta mixedwood forests[END_REF][START_REF] Wulder | Estimating the probability of mountain pine beetle red-attack damage[END_REF].

Further, the logistic regression allows obtaining the relative risk, i.e. the increase or decrease in the damage risk by a particular factor. This approach allows quantifying the effect of forest management on the damage risk which is an advantage in using model risk as a tool in management decisions (e.g. González et al., 2006). .

A two-step approach was used and two models were obtained. Both models were validated with good results. The probability of tree mortality predicted by the model has a high precision (0.917) and low bias (0.02) indicating a good fit. The estimated proportion of tree mortality for Eucalyptus globulus stands has an accuracy of 78%. Similar values of accuracy were found in other models developed to estimate tree damage (e.g. [START_REF] Wulder | Estimating the probability of mountain pine beetle red-attack damage[END_REF].

One purpose of risk models is to estimate the mortality that will happen under specific conditions, such as site conditions and stand properties. The probability of mortality does not give us a certain answer if the damage will really happen but gives us a clue about the factors affecting mortality and how those can be modified through forest management. The parameters included in the present model represent easily measurable forest characteristics such as regeneration mode, stand structure, stand age, stand composition, stand density, stand basal area and stand dominant height. Those stand variables are under the control of the forest manager, which allow a forest property-owner to reduce mortality risk. The model enables the quantification of the effects of management on the probability of mortality and therefore allows for the development of silvicultural strategies that take the risk of mortality into account. [START_REF] Schelhaas | Natural disturbances in the European forests in the 19th and 20th centuries[END_REF] in a paper based on a literature review points out the importance of including risk into management plans and how that can reduce damage susceptibility. In the present model, stand traits such as tree age, dominant tree height, age structure and stand composition (mixed versus pure stands) were found significant variables indicating that forest managers may reduce tree mortality by interfering in these parameters.

The model also allows understanding the most relevant environmental factors driving tree mortality. According to the model, mortality risk increases with temperature and decreases with evapotranspiration and frost. This outcome indicates that drought stress might be a key factor for Eucalypts mortality in the studied area. Drought may either cause directly tree mortality or predispose trees to secondary biotic agents which may be the ultimate cause of death, such as stem borers. It is known that the Phoracantha spp. stem borers of Eucalypts synchronize their egg lay with the period of maximum water stress, when low bark moisture content increases survival of young larvae [START_REF] Hanks | Water relations of host trees and resistance to the phloem-boring beetle Phoracantha semipunctata F. (Coleoptera: Cerambycidae)[END_REF]. From a forest management point of view this highlights the relevance of site selection on stand vulnerability as pointed up by Jactel et al. (2009, in press).

Damage and mortality attributable to herbivore insects are observed to be generally less intense in mixed forests than in pure forests (Jactel and Brockerhoff, 2007). The concentration of food resources for the herbivore insects, increasing host finding probability and the decrease of natural enemies are some of he main factors justifying this outcome. A similar result is supported by the present model for the tree mortality on exotic eucalypts plantations, predicting that mixed forest with native trees species, mainly pines and oaks, decrease the eucalypts mortality risk in comparison to pure plantations.

On the contrary, age differentiation between trees did not decrease the risk of tree mortality; even aged plantations are estimated to have lower risk of mortality than uneven aged plantations.

The models predict a higher mortality in older stands as well as in high forest than in young and in coppice regime. Since the stands are under forest management with regular sanitary interventions the higher percentage of dead trees in older stands reflects an increase of tree mortality with stand age. This result can be justified by an increase number of biotic factors causing tree mortality in mature trees. This result suggest that young stands, and coppice regime may be less prone to mortality, agreeing that forest planatations may suffer from less mortality than more natural stands with older tree ages. This could be one of the reasons why losses through mortality attributable to pests, diseases, and abiotic disturbances are generally less in intensively managed, planted forests than in natural forests as observed by Gadgil & Bain (1999).

Models integration in forest simulators

To conclude a logistic regression model applied to the mortality of eucalypts tree in Portugal allowed to provide relevant insights to the main factors driving tree mortality in the studied area. The model may be a further useful tool as input for regional simulator of forest growth allowing predicting future scenarios under variable environmental conditions and forest management options

Barreiro, S., and M.

Abstract

In order to obtain information on the impact of the eucalyptus snout beetle, Gonipterus scutellatus Gyll. (Coleoptera: Curculionidae), a monitoring programme was established to evaluate the densities of the populations, the levels of damage and rates of parasitism by Anaphes nitens Girault (Hymenoptera: Mymaridae). Thirty five stands, located in the centre and north of the country, were surveyed. A positive relationship between the intensity of attack and elevation was found above 450m a.s.l.. Defoliation reached 90% at 800-900m a.s.l., being the highest values observed. Furthermore, an inverse relationship between damage and rates of parasitism of eggs by A. nitens was found. An analysis of inventory data allowed for an estimation of the damage caused by G. scutellatus. In the worst affected areas, showing 30 -75% of the upper third of the crown defoliated, a significant decline in timber production, of around 40% was predicted. For defoliation above 75%, timber production losses increased substantially, reaching 90%. Results point to the importance of developing efficient management practices for the control of the eucalyptus snout beetle.

Introduction

The eucalyptus snout beetle, Gonipterus scutellatus Gyll. (Coleoptera: Curculionidae), is an insect native from Australia that feeds upon the eucalyptus' (Eucalyptus spp.) leaves. The presence of this insect in Portugal was detected for the first time in 1995, in the north of the country (Sousa e Ferreira, 1996). The eucalyptus snout beetle consumes mostly newly formed adult leaves, thus the eucalyptus more vulnerable to attacks are those in the state of transition from the young leaves to the adult leaves, between 2 and 4 years old. The damages are caused both by adults, that leave a circular patters in the periphery of the leaves, and by the larvae, that make feeding trails on the limbo of the new leaves, moving after that to the petioles and older leaves. High densities may lead to the total defoliation of the terminal branches, resulting in a strong decrease of growth and wood return. Additionally, trees consecutively defoliated may become vulnerable to attacks from other biotic agents.

In most geographic regions, the populations of G. scutellatus are managed through biological control, which was effectively achieved by the egg parasitoid Anaphes nitens Girault (Hymenoptera: Mymaridae), also native from Australia. Effectively, this parasitoid allowed the reduction of damage caused by the pest in most of the affected regions after its introduction, namely: Argentina, Brazil, New Zealand, Italy, and South Africa, coming to substitute the use of insecticides (e.g. Hanks et al., 2000). In the Iberian Peninsula, A. nitens also established effectively in the temperate areas located at low elevations. Yet, the populations of the parasitoid couldn't settle effectively in higher elevations and in colder regions during winter time, which includes most of the centre and north of Portugal, where the eucalyptus is cultivated [START_REF] Valente | First record of the Australian psyllid Ctenarytaina spatulata Taylor (Homoptera: Psyllidae) in Europe[END_REF]. In these regions the damage caused by the eucalyptus snout beetle remains high and economic losses are most disturbing for forest owners and industry.

In this work we had two main goals: i) to estimate the population densities of G. scutellatus, the level of damage (defoliation) and its relation to the rate of parasitism by A. nitens in plantations in the centre and north of the country. ii) to estimate the economic loss due to the reduction of timber production as a result of this insect activity, and relate it to the damage observed.

Material and methods

Estimation of the density of the populations of G. scutellatus and defoliation

During the year of 2007, 35 stands of eucalyptus were monitored, located in the northern regions, counties of Arouca and Castelo de Paiva, and centre, counties of Castelo Branco, Águeda, Pampilhosa da Serra, Góis and Figueiró dos Vinhos. In most of these plantations releases of A. nitens were made in 2007.

To estimate the population density of G. scutellatus and the damages, we selected a plot for sampling per each 50ha of eucalyptus plantation, resulting in 51 plots. In each sample plot, an initial point was chosen randomly from which 3 trees were selected along a row, about 125m equidistant. In these 3 trees we collected egg masses; count the density of G. scutellatus larvae and adult insects per tree crown, and evaluated the damage expressed in percentage of tree defoliation. Additionally, the damage, i.e. tree defoliation, was evaluated in 7 more trees, 5m distant between each other, starting after the third tree. This way, 10 trees were evaluated per sampling plot. Defoliation was estimated through the evaluation of the percentage of defoliation in upper third of the canopy. We used the following scale of 3 levels of damage: 1. <30% of defoliation; 2. 30 to 75% of defoliation; 3. >75% of defoliation.

To estimate the number of adult G. scutellatus per tree crown, the trees were bended and shaken over a plastic sleeve 2m wide. Whenever it was impossible to bend the tree, the plastic was placed around the base of the tree and then the tree was hit with a cane. We categorised the relative density of the adults according to the following scale: 0. no adults; 1. 1 to 20 adults; 2. 21 to 30 adults; 3. > 30 adults, from which the number of adults per tree was estimated.

In each tree, the number of young branches was also estimated. In two branches of the tree with newly formed leafs the number of larvae was counted, making possible the estimation of the number of larvae per tree crown.

The collecting of the egg capsules and the samplings of larvae and adults were made in February 2007 (before the weevil activity reached its peak), while the evaluation of the damage was carried out in June of the same year.

Parasitism rates

To determine the parasitism rates, we collected in the 51 sampling plots, whenever possible, 10 fresh egg masses per tree on the 3 sampled trees (a maximum of 30 egg masses per sampling plot).

In laboratory, we put one egg mass per Petri dish, appropriately closed with parafilm and identified, and stored it in a climatic chamber, prepared for this purpose, with a diurnal temperature of 22 ± 0.5ºC and nocturnal of 20 ± 0.5ºC (photoperiod of 16h/8h). The relative humidity was kept close to 55%. After 22 days, for each egg mass, every emergence of A. nitens and hatches of G. scutellatus larvae were registered. Observations were repeated at 30 days after incubation.

The value of the parasitism rate was given by the percentage of egg capsules from which A. nitens emerged.

Timber production losses

For this analysis we compared the evolution of the wood volume production, before and after the settling of the weevil in plantations heavily attacked or moderately attacked by G. scutellatus, by using inventories data records from the pulp industry CELBI-Silvicaima.

Stands in plantations that currently suffer attack of G. scutellatus, i.e. more than 10 % defoliation, from now on are named A stands, whereas stands in plantations that currently exhibit irrelevant or no damage are named N stands. Altogether we used the data from the inventories of 18 not attacked stands (N) and 13 attacked stands (A), in relation to two different time periods: (1) 1995 to 1997, before the presence of G. scutellatus and (2) 2004 to 2006, after the establishment of G. scutellatus. Using inventory data about stand age, dominant tree height and volume of useable timber we estimated, for each stand, which would be the wood production at 10 years old stanss, this is the reference age, corresponding to the maximum annual growth of the eucalyptus stands. For the estimations we used a growth model developed by Margarida Tomé, especially for the Silvicaima enterprise, "Silvicaima Production Model", witch was adopted by this company. For this purpose, we used the Season Quality Index (SQI), given by the dominant height of the plantation at 10 years old. The correlation between the dominant height and the volume of useable timber was used as an indicator of the quality of the model adjustment. The positive correlation obtained, r = 0.83, P <0.05, allowed us to estimate the timber volumes of the plantations at that age through this model approach.

Using the estimations of the useable volume of wood at 10 years old stands achieved in the two periods of time (1995-97 and 2004-06), in the two types of stands (A and N), we calculated the loss of useable timber caused by the presence of G. scutellatus, in percentage from one period to the other. The production loss in each stand was further related to the level of defoliation damage level observed in those stands, in one field monitoring conducted in 2006.

Results

Density of the populations of G. scutellatus, intensity of defoliation and parasitism rates

The monitoring data, in the spring of 2007, as expected evidences an increase intensity of damage for increase density of adults, larvae and egg masses (Fig. 1). The lower density of larvae registered in this period, in comparison with adults and egg masses, is due to the fact that the evaluations were made before the peak of the larval activity, in the beginning of the weevil egg laying activity. A significant negative correlation (r= 0.70; p <0.05) was found between the rates of parasitism by A. nitens and the elevation. At 300m of elevation the rates of parasitism were observed to vary from 90 to 100% of egg masses parasitized. The parasitism decreased linearly with elevation, reaching values between 0 and 10% for elevations over 700 (Fig. 2).

Similarly, a linear increase of defoliation damage was observed with the elevation (Fig. 3).

The average damage was recorded, generally, over 50% of defoliation for elevations above 600m.

Timber production loss

From the inventory data collected for the two groups of stands defined previously, we could observe a decrease in the productivity, revealed both in the dominant height and in the volume of useable timber at 10 years, in the stands that were at present attacked by G. scutellatus (A).

On the contrary, no significant decrease in timber production was observed on the other stands (N) (Fig. 4 e 5). The estimated losses in the highly attacked plantations were about 29% for the average dominant height and about 51% in volume of useable timber at 10 year old plantations. Before the arrival of G. scutellatus dominant height and wood volume at 10 years old stands was identical for the two types of stands (Fig. 4 e 5).

We also verified that the losses of volume caused by the G. scutellatus are strictly related to the intensity of the defoliation caused by this pest. In the plantations where the trees have a level of defoliation between 30 and 75%, the estimation of the decrease of timber production was about 40%, whereas the loss in the wood production increased substantially up to 90% in useable timber, when defoliation levels were over 75% (Fig. 6).

Conclusions

This study attests that the density of the eucalyptus snout beetle, as well as the damage in the canopy caused by this insect, are reduced, with percentages of defoliation close to zero, in elevations below 450m, as shown in other papers [START_REF] Valente | First record of the Australian psyllid Ctenarytaina spatulata Taylor (Homoptera: Psyllidae) in Europe[END_REF]. Nevertheless, the defoliation increases linearly with the elevation, reaching average values over 75% in elevations above 700m. Likewise, the parasitism rates by A. nitens decrease linearly with the elevation, which indicates an increase off the difficulty of the A. nitens settlement as the elevation grows, probably related to the low temperatures registered during the Winter in these regions. On the other hand, these results attest the importance of the parasitoid in the biologic control of the eucalyptus snout beetle, visible in the low damages caused by the de G. scutellatus in the plantations located at elevations below 450m, as observed in other temperate regions where the parasitoid was introduced (e.g. Hanks et al., 2000). Nonetheless, currently about 28% of the total area of eucalyptus plantations are located in the centre and north of Portugal above 450m of elevation, from which the damages are considerable.

The results show that the defoliation caused by the snout beetle have a quite significant impact on the productivity of the eucalyptus plantations, since it can lead to significant economic losses, estimated in about 40% to 90% of losses in the volume of useable timber, for medium and severe defoliations, respectively. Actually, the loss of volume may be even bigger than the one estimated in this study, since the model used to extrapolate the data is only based on the height and diameter measured at a given moment in the tree's age. In fact, the used model assumes rates of growth that doesn't take in consideration the effect of defoliation in the future growth of the trees, and estimates the growth as if the tree had never been defoliated. However, the consecutive reduction of the leaves surface, essential for the tree's grows, will reflect negatively also in the following years. Additionally, there is the possibility that the defoliation continues in the subsequent years to the making of the inventories. Because of these two reasons the method underestimates the real losses, which will be greater than those values. The high losses of productivity indicate the need of intervention with adequate management measures. Valente, C., Vaz, A., Pina, J., Manta, A., Sequeira, A. 2004. Control

Predisposition of Norway spruce (Picea abies) dominated forests to the European spruce bark beetle, Ips typographus in Austria

Authors: Netherer S, and Pennerstorfer, J.

Objectives

The main objective of the specific risk analysis is to evaluate the effect of new stand management alternatives (according to FMAs defined by Eforwood) on the level of bark beetle induced stand mortality in a central European forest dominated by Norway spruce. The modelling is based on the hypotheses that FMAs have an effect on some forest stand characteristics (e.g. composition and structure) and that certain stand characteristics drive susceptibility to bark beetle infestation. The damages caused by bark beetles result in wood volume loss through tree mortality.

Study area

The specific risk analysis was based on data provided by Kalkalpen National park, a large protected forest area situated in the Upper Austrian limestone Alps. The forest stands are dominated by Norway spruce (Picea abies), fir (Abies abies) and beech (Fagus sylvatica).

Picture source: http://www.kalkalpen.at/

Methods

The specific risk analysis was performed by implementing a knowledge-based expert model for hazard assessment developed at our institute. Under the assumption that forest susceptibility to damaging agents is influenced by specific site and stand characteristics [START_REF] Berryman | Forest Insects. Principles and practice of population management[END_REF], key factors for the predisposition of forests to diverse biotic and abiotic agents were identified by literature research and expert discussions and summarized by way of checklists [START_REF] Netherer | Assessment of the predisposition of forest sites and stands to epidemics of the Little Spruce Sawfly Pristiphora abietina (Christ) (Hym. Tenthr.)[END_REF][START_REF] Nopp | Erarbeitung von Identifikationsschlüsseln der Prädisposition fichtenreicher Bestände gegenüber verschiedenen abiotischen und biotischen Schadauslösern[END_REF], Nopp et al. 2001, Führer & Nopp, 2001). The assessment systems have already been applied and verified in several research areas (Ips typographus: [START_REF] Netherer | Monitoring and risk assessment of the spruce bark beetle Ips typographus[END_REF][START_REF] Netherer | Modelling of bark beetle development and of site-and stand-related predisposition to Ips typographus (L.) (Coleoptera; Scolytidae). A contribution to Risk Assessment[END_REF], Netherer & Nopp-Mayr 2005[START_REF] Schopf | Risikoabschätzung von Borkenkäfer-Massenvermehrungen im Nationalpark Kalkalpen (Risk assessment of bark beetle massoutbreaks in the National Park Limestone Alps[END_REF].

For the present analysis, predisposition to the spruce bark beetle as well as to storm and snow damage was evaluated for each Eforwood-forest management alternative (FMA) from closeto-nature management to biomass production. The analysis was performed in time steps of 20 years from now to 100 years in future in order to compare susceptibility trends between the different forest management regimes depending on stand age. The original dataset of the study area, provided by the management of the national park, represents the FMA "nature reserve". According to the management strategies of the different other alternatives the parameters "proportion of Norway spruce", "rotation period", and "canopy closure" were modified (see table 1). Assessment of predisposition to Ips typographus -relevant criteria Assessment of predisposition was based on a simplified expert model comprising the parameters species composition, stand age, stand structure, predisposition to storm damage and snow breakage, which most fundamentally affect stand-related susceptibility to the spruce bark beetle (comp. Fig. 1234). The causalities are shortly summarised in the following: typographus meets ideal developmental conditions in monocultures of Norway spruce. High availability of breeding material is combined with low population densities of antagonists, e.g. due to the lack of herbal or shrub layers. The incidence of bark beetle infestations markedly declines with a proportion of deciduous or conifer trees other than Norway spruce of at least 50%.

Damage is mostly concentrated in stands older than 60 years and most likely for high age classes. Especially mature and over-mature forests of reduced vitality are disposed to dieback caused by bark beetle infestation.

Generally, structural features leading to increased solar irradiation and temperature conditions in the stand enhance bark beetle development. Forest edges, the margins of clear-cuts, small clearings or roads are favoured locations of primary infestation, as well as stands of low crown cover. Intense thinning may be followed by an abrupt increase in solar irradiation and sudden stress of the remaining trees and should be avoided. On the other hand, too dense stocking of trees again results in higher susceptibility to bark beetle attack, because of reduced stand vitality due to competition.

The sudden blow down of tree groups or entire stands in case of storm events or snow breakage of trees results in an excess supply of food and breeding material for the spruce bark beetle. Usually beetle populations build up during 1-2 years after the incident. Storm events create gaps of different size in the affected stands, where solar irradiation and thus temperature increases. Sun-exposed stems lying on the ground show significant higher bark temperature than standing trees. The remaining standing trees are unconditioned to the changed site conditions and, as a result of stress, show less resistance to bark beetle attack.

Weighting

The stand-related criteria "predisposition to storm damage" and "predisposition to snow damage" are again output of (also simplified) expert models.

Parameters considered for the evaluation of stand-related predisposition to storm damage: possible to evaluate probable wood volume loss in this study, as the necessary data were not provided by the management of the national park.

Assumed "close-to-nature" management (CN) of the forest stands resulted in the most heterogeneous pattern of stand-related predisposition due to variable species composition and canopy closure, followed by "combined objectives" (CO) (comp. Fig. 4.4-7). CN was the only alternative to show areas of low predisposition and exhibited the lowest extent of high predisposition areas over time. The most homogeneous pattern of predisposition was obtained for the management alternative "biomass production" (BIOM), as the proportion of spruce was generally set 100% and canopy closure 60%.

All management alternatives except for BIOM showed a clear dependency of stand-related predisposition on time.

Conclusions

The trend of stand-related predisposition over time is mainly determined by the parameter stand age. The shorter the rotation periods and consequently the younger the stands are in general, the less heterogeneous are the patterns of predisposition. Structural diversity is lowest under the objective "biomass production", with uniform species composition and canopy closure, and age classes only ranging between <=20 yrs and 21-40yrs. This constellation leads to a constant number of forest districts highly predisposed to bark-beetle induced mortality within the different time periods and over time (comp. Fig. 4.8).

The likelihood of stand mortality is lowest for "close-to-nature" management due to naturally mixed species composition and heterogeneous canopy closure. However, as a rotation period of 100 years is assumed for this management alternative, probabilities of biotic disturbance are increasing with time and highest in mature forest parts. This trend can also be observed for the variant "combined objectives", whereby a much larger forest area is prone to barkbeetle induced mortality due to the assumed dominance of Norway spruce.

While CN management implies a slow increase of stand-related predisposition with stand age, intensive even-aged forestry might result in an abrupt increase in areas of highly vulnerable stands at the time when the age classes susceptible to bark beetle infestation are reached. Predisposition will then remain on high levels until the end of the rotation period.

Estimation of the effects of new forest management alternatives on wood volume loss using hazard-specific risk analyses. Case: Silesia (Beskidy Mts.), Norway spruce, Ips typographus (L.)

Wojciech Grodzki, Forest Research Institute (IBL), Kraków, Poland

Materials and methods

The data on tree mortality due to Ips typographus (L.) infestation in individual forest subcompartments were used. The relative (%) mortality was calculated by the comparison of the volume of infested trees processed during a year with the volume stock of Norway spruce in a given sub-compartment (record). For each record the number of species and the share of Norway spruce and broadleaved species (in 10% intervals) were assigned using forest inventory data. The slope exposure and the mean altitude of the sub-compartment (centroid of a polygon) was derived from the forest digital map (DTM generated using izolines). For the analysis only the sub-compartments with recorded occurrence of infested trees in 2006 were selected. The effect of 5 factors: number of species, share of spruce, share of broadleaves, exposure (8 geographical directions) and altitude (4 zones: <600, 601-800, 801-1000, >1000 m a.s.l.) was assessed using Kruskal-Wallis test. The non-parametric test was used due to the lack of normal distribution of analyzed data (concentration in the lowest classes). The mean age, DBH and height was calculated using weighted averages of the spruce fractions occurring in each sub-compartment, then the 20-years age classes and 10-cm DBH classes sere used as grouping variables in the analyses.

The relative tree mortality in 2006 was analyzed in two objects: Forest District Węgierska Górka and Ujsoły, located in the Western Carpathians in the area of a large Ips typographus (L.) outbreak and related spruce decline [START_REF] Grodzki | Spatio-temporal patterns of the Norway spruce decline in the Western Beskidy mountains in Poland[END_REF][START_REF] Grodzki | Forest decline in mountain spruce stands affected by bark beetle outbreaks in Poland. Attempt to spatial characteristics[END_REF].

Object 1: Forest District Węgierska Górka -forest area 9 029 ha, pure spruce 3 913 ha (43,3%). Share of spruce (pure and dominated) stands according to Forest Inventory data in 2007: 83,0%. The forst bark beetle outbreak took place in 1997-2001 [START_REF] Grodzki | Threats to mountain Norway spruce stands in the Carpathians from the insect pests[END_REF]) and the second -from 2003 till present, with culmination in 2007 [START_REF] Grodzki | Forest decline in mountain spruce stands affected by bark beetle outbreaks in Poland. Attempt to spatial characteristics[END_REF]. The data from 527 out of 1741 sub-compartments included into the analysis.

Forest District Ujsoły -forest area 12 886 ha, pure spruce 7 643 ha (59,3%). Share of spruce (pure and dominated) stands according to Forest Inventory data in 2007: 92,5%. Stable bark beetle population lasted till 2002 [START_REF] Grodzki | Threats to mountain Norway spruce stands in the Carpathians from the insect pests[END_REF], then the outbreak, started in 2003 and still ongoing, has the culmination in 2007 [START_REF] Grodzki | Forest decline in mountain spruce stands affected by bark beetle outbreaks in Poland. Attempt to spatial characteristics[END_REF]. The data from 619 out of 1892 subcompartments included into the analysis.

Object 2. Although some factors (number of species in the stand, share of spruce, altitude, age and BHD of spruce) have a significant effect on the tree mortality (tab. 1), the relationships are not linear (fig. 12345678), thus the reliable regression equation can not be defined. In this case the effect of some factors (number of species in the stand, share of spruce and broadleaved species, altitude) on the tree mortality is significant (tab. 2), the relationships are not linear (fig. 101112131415161718), thus the regression equation can not be defined.

Results

Object

Object 2.

General remarks

-the analyzed data concern outbreak conditions -no evident patterns can be defined -every spruce can be attacked in epidemic bark beetle population level [START_REF] Christiansen | Resistance of conifers to bark beetle attack: searching for general relationships[END_REF]) -in outbreak conditions probably there are no relationships with forest management alternative -similar situation occurred in actively or passively protected stands in national parks [START_REF] Grodzki | Windthrowns and spruce bark beetles in protected areas in Polish mountains[END_REF])

Effect of Pine Processionary Moth defoliation on tree growth: comparison of empirical, dendrochronogical and meta-analytical outcomes

Jacquet, J.S., Jactel, H.

Abstract

In the context of climate change, contribution of forests to carbon sequestration is widely acknowledged. However this function is often adversely affected by biotic damage such as defoliation by pest insects and the risk is likely to increase with damaging agents that can benefit from warmer temperatures. The Pine Processionary Moth (PPM) is the most damaging pine defoliator in southern Europe and it is now spreading to higher latitudes and elevations due to the effects of global warming. To better predict the effect of future PPM outbreaks on forest productivity we aimed at evaluating the impact of PPM defoliation on tree growth.

We used three methodological approaches. First we correlated loss in radial growth to level of PPM defoliation at the tree level. Second we performed a dendrochronological analysis to relate individual ring increments to past records of PPM infestation at the stand level. Lastly we used the meta-analysis method to combine outcomes of published studies on the effect of Thaumetopoea spp. on Pinus and Cedrus tree growth.

The three methods provided consistent results. They showed that, overall, defoliation by processionary caterpillars resulted in a reduction of tree growth that linearly increased with the rate of defoliation. However the growth loss was not significantly different from zero in the lower class of defoliation, i.e. 1 -25%. There was no significant interaction between classes of defoliation and initial tree diameter on the effect on tree growth. It took only 1 to 2 years for the trees to recover from a single defoliation.

Among the tested methods, the use of dendrochronological analyses combined with quantile regressions seems promising as it allows detecting past effects of insect defoliation on tree growth even if no individual tree survey has been conducted.

Introduction

In the current global warming context, contribution of forest to carbon sequestration and biofuel production is widely acknowledged (Laurent Saint-André, 2008). However these functions may be adversely affected by biotic agents that damage tree crown and reduce photosynthesis. For example tree defoliator can cause consistent reduction of effective leaf area. And yet, climate change, with increasing temperature, is likely to pest insects as poikilotherm organisms [START_REF] Netherer | Potential effects of climate change on insect herbivores in European forests-General aspects and the pine procassionary moth as specific example[END_REF]. There is therefore an urgent need to better evaluate the impact of insect defoliator on tree growth, mainly for those insect species models that are relevant to climate change issues.

The pine processionary moth (Thaumetopoea pityocampa Dennis & Schiff.) is one of the most harmful insects that defoliate Pinus and Cedrus species in Mediterranean countries [START_REF] Demolin | Incidence of some factors affecting the social behaviour of larvae of Thaumetopoea pityocampa Schiff. (Lepidoptera) during the period of the pupation procession. Their repercussions on the effectiveness of parasites in 'Biological Control in the Forest', at Pont-a-Mousson[END_REF][START_REF] Kanat | Thaumetopoeiadae) on annual diameter increment of Pinus brutia Ten[END_REF], Markalas 1996). It is a univoltine oligophagous moth species with the unique particularity of a winter larval development. There is usually one generation per year. Moths are present from July to September. Eggs are laid in batches, wrapped around a couple of needles. Caterpillars live in groups. Larvae go through five instars. Young caterpillars make small provisional nests. They feed on nearby needles and cause yellow-brownish discoloration. In late autumn, caterpillars spin a final silky nest (winter nest) in the most illuminated part of the crown. At the end of winter, caterpillars walk in a procession towards warm spots on the soil, where they dig themselves in and pupate. Pupae may stop their development (diapause) during several months to several years (prolonged diapause of up to three years). Outbreaks are periodic, occurring every 5 to 7 years [START_REF] Robinet | Modelling the effects of climate change on the potential feeding activity of Thaumetopoea pityocampa (Den. & Schiff.) (Lep., Notodontidae) in France[END_REF]). Unfavourable climatic conditions, starvation and natural enemies are the main factors regulating populations. The pupal diapause has been considered as a way to escape these constraints. Old larva (from the third instar) are the most destructive, they can cause complete defoliation on vast surfaces.

The pine processionary moth (PPM) Thaumetopoea pityocampa Denis & Schiffermüller (Lepidoptera: Thaumetopoeidae) is spreading towards higher latitudes and altitudes, probably due to the global warming [START_REF] Battisti | Expansion of geographic range in the pine processionary moth caused by increased winter temperatures[END_REF][START_REF] Robinet | Modelling the effects of climate change on the potential feeding activity of Thaumetopoea pityocampa (Den. & Schiff.) (Lep., Notodontidae) in France[END_REF]. PPM females use visual cues to locate host trees, flying towards dark tree silhouette against clear backgrounds [START_REF] Demolin | Incidence of some factors affecting the social behaviour of larvae of Thaumetopoea pityocampa Schiff. (Lepidoptera) during the period of the pupation procession. Their repercussions on the effectiveness of parasites in 'Biological Control in the Forest', at Pont-a-Mousson[END_REF].

PPM outbreaks constitute a health hazard to humans and domestic animals, due to urticating hairs. However the main impact of T. pityocampa is on tree growth leading to potential economic losses and ecological impacts [START_REF] Hodar | Host utilisation by moth and larval survival of pine processionary caterpillar Thaumetopoea pityocampa in relation to food quality in three Pinus species[END_REF]. Defoliations result in reduction of needles activity and/or availability for the photosynthesis. Severe and repeated defoliation may even lead to tree mortality, mainly on young trees and poor soil conditions. Several studies have been carried out to evaluate the loss in radial growth [START_REF] Laurent-Hervouët | Mesure des pertes de croissance radiale sur quelques espèces de Pinus dues à deux défoliateurs forestiers. I Cas de la processionnaire du pin en région méditerranéenne[END_REF][START_REF] Lemoine | Contribution à la mesure des pertes de production causées par la chenille processionnaire Thaumetopoea pityocampa Schiff au pin maritime dans les Landes de Gascogne[END_REF][START_REF] Babur | The effect of damage of Thaumetopoea pityocampa (Schiff.) in Calabrian pine seedling in Turkey, symposium on problems Thaumetopoea pityocampa (Schiff.) in Turkey and their solutions[END_REF], Gachi et al. 2005[START_REF] Kanat | Thaumetopoeiadae) on annual diameter increment of Pinus brutia Ten[END_REF]), volume growth [START_REF] Bouchon | Etude préliminaire sur les pertes de production des pinèdes soumises aux attaques de Thaumetopoea pityocampa Schiff[END_REF], Graph et al. 1995) and height growth (Markalas 1996) caused by T. pityocampa's defoliation.

However, several key questions remain to solve:

-Is there any threshold value of defoliation under which no effect on tree growth would be detected? -Does the effect of defoliation on tree growth also vary with cofactors such as tree size (or tree vigour)? -How long would last the effect of a single event of defoliation (or how fast can trees recover from one defoliation)? -Most of the methods used to evaluate the effect on tree growth are individually -based (survey of individual trees before and after defoliation): can we develop a retrospective and stand -based method that could be more widely applied? -Using very different methodologies makes it difficult to compare their results: are these results consistent?

In the present study, we present the results from three different methodologies: 1) an empirical approach: it is based on data collected by Lemoine in 1970-71 in maritime pine (Pinus pinaster) stands that we re-analysed with modern statistics; 2) a dendrochronological approach: The French Forest Health Department has developed a permanent survey of Pine Processionary Moth infestations that provide information about defoliation and number of winter nests in many pine stands for the last twenty years or so. However the effect on tree growth is not evaluated in this network. In six of its plots, we collected cores and retrospectively calculated individual tree growth. Then we tried to correlate annual growth of individual trees in a given year (N values) with the density of PPM nests (one value) in the same plot and year; 3) a meta-analysis approach: Meta-analysis is a powerful method to combine results from several independent studies into a quantitative and objective review. We used this tool to compare tree growth between defoliated trees and control trees and to evaluate the magnitude of the effect on tree growth for different levels of PPM defoliation. We used 20 cases studies reported in the literature about Thaumetopoea pityocampa (17 cases), T. bonjeani (1 case) and T. wilkinsoni (2 cases) on Pinus pinaster, P. brutia and Cedrus atlantica.

Each of these methods and associated outcomes will be presented successively. Then a synthesis will be drawn to compare results, provide new insights in the effect of PPM defoliation on pine tree growth and propose further studies to clarify some specific issues.

Empirical analysis based on individual tree data

Re-analysis of data collected by Lemoine in 1970 on two maritime pine stands

The study of [START_REF] Lemoine | Contribution à la mesure des pertes de production causées par la chenille processionnaire Thaumetopoea pityocampa Schiff au pin maritime dans les Landes de Gascogne[END_REF] was conducted in the "Landes de Gascogne" area, a one million hectares plantation forest of Pinus pinaster. The objective of this study was to estimate the effect of severe PPM defoliation of individual tree growth. The main attack took place in 1969 with a defoliation rate ranging from 0 to 100% of crown defoliation at the tree level, followed by a second year of lower defoliation in 1970. A total of 773 trees were surveyed in 4 different plots. To get a balanced sample of trees among four defoliation classes (0 -5%, 6 -25%, 26 -50%, 51 -100% defoliation) we restricted our analysis to only two representative plots: "Parentis" (226 sampled trees) and "Leon". (141 sampled trees). The Parentis stand was 20 years old in 1969 with an average diameter of 11,0cm (and average height of 69,7m) whereas the Leon stand was 45 years old with an average diameter of 32,2cm (and average height of 21,16m).

For each sampled tree we used the original data sheets to enter individual diameter and defoliation values, from 1970 to 1975, in a new access file. Then we recalculated the radial growth area as the main response variable to PPM defoliation.

The radial growth area (RGA n+1 ) in year n+1 was calculated as the difference between the area in the year n+1 and the area in the year n:

RGA n+1 = π × R² n+1 -π × R² n
With R= CIR/2 π (CIR =circumference at Breast Height 1.3m)

Statistical Analyses

For the data analysis we used a two steps approach.

We started with data concerning the effect of the 1969 outbreak on tree growth in the first year following this attack. We made three tree circumference classes of equal size (same number of trees per class within each stand) in order to test the effect of PPM defoliation (4 classes), initial tree circumference (in 1969) and their interaction on radial tree growth (RGA n+1 ) with a two-ways ANOVA. Beforehand RGA data were transformed with a Box-Cox transformation (Lambda = 0.01) to normalize the residuals. Then we analysed separately the data from the three circumferences classes and used a Dunett's unilateral test to compare the mean tree growth (RGA) within each defoliation class (6 -25%, 26 -50%, 51 -100% defoliation) with the growth in control trees (0 -5% defoliation).

In a second step we tried to estimate the recovery time from a single event of defoliation. We selected trees that were attacked only once (in 1969) and not the following years (e.g. 1970) in order to avoid any confounding effects of repeated defoliations. In these trees we used the RGA calculated for each of the four years following the initial PPM attack as response variables (i.e. RGA n+1 , RGA n+2 , RGA n+3 , RGA n+4 ). We made two-ways ANOVAs to test the effect of circumference class, defoliation class and their interactions on RGA n+i . Then we pooled the data from the three circumferences classes and used a Dunett's unilateral test to compare the mean tree growth (RGA) within each defoliation class (6 -25%, 26 -50%, 51 -100% defoliation) with the growth in control trees (0 -5% defoliation) for each year following the defoliation.

Results and discussion

Effect of PPM defoliation on tree growth in the first year following the defoliation

As shown in Table 1, both in initial tree circumference and PPM class of defoliation had a significant effect on radial tree growth. However the interaction between these variables was not statistically significant. PPM defoliation had therefore similar effect on tree growth irrespective to the initial circumference of attacked trees. Then we calculated a Relative Growth Loss (RGL) as following:

ctl ctl def RGA RGA RAG RGL - =
with RAG def the radial growth area of defoliated trees and RAG ctl the radial growth area of control trees (0-5% defoliation).

For both plots, we observed the same gradual effect of defoliation by Thaumetopoea pityocampa on tree growth irrespective to initial tree circumference (Fig. 1). Defoliation resulted in a reduction of tree growth and the growth loss always increased with defoliation intensity. However the growth in defoliated trees was significantly lower than the growth in control trees only for a defoliation rate higher than 50% (with one exception of large trees in the Parentis stand).

Tree growth was significantly higher in large than in medium and in small trees (data not shown). This result indicates that control trees should be selected within the same class of diameter than defoliated trees to make proper comparison of growths. This recommendation is not that trivial because trees attacked by T. pityocampa are on average higher, and then larger, than unattacked trees due to the selection behaviour of female moth that preferably lay eggs on taller trees (Démolin 1969). Asterisked bar indicates that the growth in defoliated trees was significantly different from the growth in control trees, according to a Dunett's unilateral test (or that the relative growth loss was significantly different from zero).

The relative growth loss differed between the two sampled stands. In highly defoliated trees (51 -100% defoliation) the growth loss was almost twice as important (70%) in the Leon stand than in the Parentis stand (40%). The difference between the two stands may be due to the age of trees which are older in the plot "Leon". Older trees have a slower growth and may not have enough reserves to compensate the reduction in carbohydrates synthesis resulting from the loss in needle biomass [START_REF] Carus | Impact of defoliation by the pine processionary moth (Thaumetopoea pityocampa) on radial, height and volume growth of Carlabrian pine (Pinus brutia) trees in Turkey[END_REF].

To summarise, we found that the PPM can cause, on average, a relative loss in tree growth of 25% in Parentis and 50% in "Leon". These results are consistent with the findings of [START_REF] Lemoine | Contribution à la mesure des pertes de production causées par la chenille processionnaire Thaumetopoea pityocampa Schiff au pin maritime dans les Landes de Gascogne[END_REF] who estimated a mean reduction of 30% in radial growth.

Effect of a single PPM defoliation on tree growth in the four years following the defoliation.

As shown in Table 2, for the trees attacked only once (in 1969), both circumference and percentage of defoliation classes still had a significant effect on tree growth the first year following the attack and the interaction of these variables was not statistically significant. In the "Parentis" stand, the defoliation classes had no significant effect on tree growth the following years. For "Leon" stand, the defoliation classes had still a significant effect on tree growth the second year following the attack, but the circumference classes had not and the interaction had no significant effect as well. Then we pooled the data from the three circumferences classes. In the "Parentis" stand, the growth loss was significantly lower in defoliated trees than in control trees only the first year following the attack. This means that the trees managed to recover already from the second year following the attack. In the "Leon" stand, the effect of defoliation on tree growth lasted two years in trees that experienced an initial defoliation higher than 25%.

These results are therefore different from the original conclusions of Lemoine in 1977 who found a three years recovery time with his own method. It is probably because we used a different methodology with more recent statistic tools. We restricted our analysis to only two representative plots to get a balanced sample of trees between defoliation classes and we selected trees that were attacked only once (in 1970) to avoid any confounding effects of repeated defoliation. Asterisked bar indicates that the growth in defoliated trees was significantly different from the growth in control trees, according to a Dunett's unilateral test (or that the relative growth loss was significantly different from zero).

3 Dendrochronological, retrospective analysis to estimate the impact of PPM defoliation on tree growth.

Dendrochronology was used to test the relationship between annual ring area growth at the tree level and defoliation by the pine processionary moth (PPM) at the stand level. We had a record of mean defoliation rate at the stand level (by Forest Health department) for the last 24 years but not at the individual tree level. We had ring area growth data for a sample of individual trees, each of the 24 years, thanks to tree cores analysis. We could then plot and correlate individual ring area growth of each tree for a given year against the mean number of nest of the corresponding year in the same stand.

We performed this analysis with data from Laurent-Hervouet (1986) on Austrian pine (Pinus nigra Arn. ssp. nigricans Host) in the Mont Ventoux region and our own data collected in 2007 -2009, in two maritime pine stands of the "Landes de Gascogne" forest.

We used the quantile regression method to test whether the PPM nests density at stand level placed an upper limit on ring area growth at tree level. To quantify this effect, and to make the results comparable with the first empirical part, the number of PPM nests was converted into percentage of defoliation using data from the Forest Health Department in maritime and Austrian pine.

Stand description and core collection

• Forest Health Department (DSF) plots in the Landes de Gascogne

Cores were collected in two plots in the Landes de Gascogne: the "Biscarosse" plot and the "Lit et Mixe" plot. The "Biscarosse" plot was a pure stand of maritime pine in the Landes district, close to Biscarosse (latitude: 315.7 and longitude: 194.7). It is located in a flat area; the slope is null and the altitude is 26m. This stand was planted in 1978 and was 29 years old at the core collection time. We sampled 60 trees at the plot edge and 112 trees within the plot. We collected trees on the edge because PPM females are attracted by the silhouette of tree crown and then often show a preference to border trees [START_REF] Demolin | Incidence of some factors affecting the social behaviour of larvae of Thaumetopoea pityocampa Schiff. (Lepidoptera) during the period of the pupation procession. Their repercussions on the effectiveness of parasites in 'Biological Control in the Forest', at Pont-a-Mousson[END_REF]. We then expected the trees on the edge to be more attacked. The cores collection was realised in October 2007.

The "Lit et Mixe" plot was a pure stand of maritime pine planted in 1983 (and was 26 years old at the core collection time) in a flat area. In June 2009, we collected cores in 100 trees within the stand and 30 trees on the edge. Interior and edge trees were pooled for this plot because of the small edge sample size.

We collected the cores manually at 1.3m above the ground level with a 5 mm Pressler auger. We collected one core per tree at the first core collection in 2007 and two cores per tree at the second collection.

After collection the cores were flattened and scanned. Pictures were analysed with Windendro, a specific software designed for dendrochronological analyses. It assists users in estimating ring width for each year of growth. The ring growth area (RGA n ) in year n was calculated as the difference between the area in year n and the area in year n-1:

RGA n+1 = π × R² n+1 -π × R² n
The ring growth area was calculated here automatically through an access request.

In each plot the number of pine processionary moth nests per ha was counted yearly since 1986 by the Forest Health Department giving data for the last 24 years.

Picture 1. Tree core ready to be analysed in Windendro. The core was flattened and scanned.

• Plots of Laurent-Hervouët (1986) at the Mont Ventoux

The study area was near the "Mont Ventoux", in the state owned forest of Malaucène. Three plots have been chosen where the forest services have estimated the mean number of nests per tree at the stand level, each spring from 1958 to 1981. The first two plots were in Pinus nigra unevenaged stands of ca. 90 years at an altitude of 700 m. They were regularly attacked by the pine processionary moth. The third plot was at an altitude of 560 meters, in a more Pinus nigra even-aged stand. This plot was on average less attacked. The tree cores had been taken in October 1983 with a 5 mm Pressler auger. For the three plots, we had 12, 16 and 13 tree cores available respectively. We had access to the archive data of the core analysis performed by Laurent-Hervouët (annual radial growth (mm) measures from 1958 to 81). We used annual radial growth (mm) as growth response variable for this stand.

Statistical analyses

For each stand, we plotted the ring growth area (or the radial growth for the Ventoux stand) of each trees in a given year (N values) against the total number of nests per ha (one value) in the same year (so that the X axis did not correspond to a time scale but to a gradient of increasing PPM infestation). This graphical analysis revealed a triangular relationship between ring growth area and the number of nests (Fig. 3). When PPM nests were abundant, the variability in ring growth area was then consistently reduced, while when the nests were rare, variability in ring growth area was higher. This means that the number of nests did not fully control tree growth but was an active limiting constraint (according to the limiting factor concept) and induced a pattern of unequal variation. The objective of our analysis was then to test whether the nest abundance placed an upper limit on individual ring growth area. Due to the effect of other, ignored driving factors, changes near the maxima of nest abundance were better estimates of the relationship than at the distribution center [START_REF] Cade | Estimating effects of limiting factors with regression quantiles[END_REF]. Ordinary least square regressions were therefore not appropriate for this kind of data (Johnson and Van Der Wal, 2009). Indeed least square regression assumes that the variance of the Y variable (ring growth area) is uniform with respect to the value of the X variable (nest abundance) witch was not the case in our study. The pertinent statistical analysis would be to test the significance of the constraint line that marks the upper boundary of the scatter diagram. The quantile regression method has been designed for such pattern analysis [START_REF] Cade | Estimating effects of limiting factors with regression quantiles[END_REF]. It was used to fit quantile value of the ring growth area as a function of the measured variable, the nest abundance. Several methods for constructing confidence intervals of slopes estimates (CI) have been developed. We used the rank inversion method (using the quantreg library in the R statistical package version 2.8.1, http://www.R-project.org) which provides, for each quantile, a 99% bootstrap confidence interval. The evolution of the magnitude of the slope and its confidence interval shows the changing effect of the number of nests per ha as a limiting factor on individual ring growth area.

Once the effect of PPM infestation on individual tree growth was ascertained, we calculated relative losses in ring growth area in order to make the results comparable with the empirical approach. We used a two steps approach: -First, the number of nests per stand was correlated with mean % of defoliation using the Forest Health Department data on maritime pine and Austrian pine plots. Using an ANOVA or a Friedman non-parametric test, we determined if this relation was significant.

We calculated a mean number of nests and its 95% Normal confidence interval for three defoliation classes: 1-5%, 6-25% and 26-60% defoliation. -Then the effect of the defoliation was quantified using the least square regressions on mean values.

Then we could compare relative loss in growth area from the empirical and the dendrochronological approaches for two defoliation classes: 6-25% and 26-60%.

Results

Biscarosse Plot

For the DSF Biscarosse plot, plotting ring growth area as a function of number of nests per ha showed a triangular pattern (Fig. 4). Least square regressions (Fig. 4) were significant (P < 0,0001 for the interior and P < 0,0001 for the edge) but explained just a very small part of variability (R² = 0.01 for the interior and 0.03 for the edge). Least square regressions slopes were respectively -1.21 and -2.50 for the interior and the edge of the plot.

The results of the different quantile regressions for the interior of the plot are shown in Figure 4a. The slope was always negative indicating a reduction in tree growth when PPM infestation increased. The magnitude of the slope, which is a measure of the size of the effect for each quantile, increased as increasingly higher quantiles data were used, up to the 0,995 th quantile (Fig. 4a). The slope estimated for the 0.995 th quantile was -8.66 and had a 99% CI (confidence interval) that did not bracket zero, indicating a slope value significantly different from zero.

For quantiles higher than 0.995 th , the regression slopes was not significantly different from zero probably due to a low number of data within these quantiles (as also shown by larger confidence intervals).

We observed the same general pattern at the edge of the plot (Figure 3b and4b) but the slopes were steeper: for example, for the 0.995 th quantile the slope was -14.46 for edge trees instead of -8.66 for the interior trees. The Figure 6 shows 0,9 th quantile regression curves for interior and edge trees of the Biscarosse plot. The regression lines are crossing because the slopes are different. The effect of defoliation seems to be higher at the edge of the plot. a.

b. Two explanations were possible.

First, trees are often more vigorous at the edge of pine stands because of higher light interception and defoliation may not have the same effect on trees of different vigours: we demonstrated that this hypothesis was not correct with the empirical approach (cf. 2.3.1, Table 1).

Second, the number of nests at the edge is probably higher than the average estimate provided by the Forest Health Department because PPM females are attracted by the silhouette of tree crown [START_REF] Demolin | Incidence of some factors affecting the social behaviour of larvae of Thaumetopoea pityocampa Schiff. (Lepidoptera) during the period of the pupation procession. Their repercussions on the effectiveness of parasites in 'Biological Control in the Forest', at Pont-a-Mousson[END_REF]) and then are preferentially laying eggs on dominant trees [START_REF] Geri | Etude sur les populations de processionnaire du pin Thaumetopoea pityocampa[END_REF]Geri 1983). To correct the underestimation of the number of nests at the edge, we calculated a correction coefficient using the following model: Using this linear model, we adjusted the edge number of nests. The Figure 7 shows the example of the 0.9 th quantile regressions after correcting the number of nest in edge trees. The effect of nests number on ring growth area is now similar at the edge and within the plot (parallel quantile regression lines). Using the equation of the least square regression on mean values, we were able to quantify the ring growth area for a determined number of nests per ha. Using data from the Forest Health Department in maritime pine plots, a correlation between the number of nests per ha and the mean percentage of defoliation by the pine processionary moth was established (Figure 7). The ANOVA analysis showed a significant effect of the class of defoliation on the number of nests per ha indicating a significant difference in nests number between classes of defoliation (no overlapping between classes). A mean number of nests and its 95%confidence interval were then calculated for three defoliation classes: 1-5% defoliation corresponded to 119 ± 39 nests per ha 6-25% defoliation corresponded to 310 ± 72 nests per ha 26-60% defoliation corresponded to 700 ± 189 nests per ha

Using the equation of the least square regression on mean values (Figure 4, red regression), corresponding ring growth area values for each of these defoliation classes were calculated. We used the least square regression on mean values, which was significant, instead of quantile regression equation because there was no objective reason to use one particular quantile or another. We could then calculate a relative growth loss (using the mean value of ring growth area in unattacked trees as reference) for each of the three defoliation classes (Figure 8). The Relative Growth Loss was calculated as following:

RGA def is the mean ring growth area of defoliated trees and RGA ctl is the mean ring growth area of controlled trees (0-5% defoliation).

For the interior of the Biscarosse plot, we got an average reduction of 25% in relative ring growth for the 6-25% defoliation class and an average reduction of 57% for the 26-60% defoliation class. As shown in the Figure 8, these results are consistent with the results from the empirical approach. The "Leon" pattern (Figure 1.b) shows an average of 19% relative ring area growth loss for the 6-25% defoliation class and of 43% ring area growth loss for the 26-60% defoliation by the PPM. The "Parentis" pattern's growth reduction is lower: 14% in average for 6-25% of defoliation and 27% for 26-60% of defoliation.

The Biscarosse plot edge relative growth losses were calculated with the least square regression after the nest number correction. We get an average reduction of 19% in relative ring area growth for the 6-25% defoliation class and an average reduction of 43% for the 26-60% defoliation class. These results are really consistent with results from the empirical approach (particularly "Leon results") Comparison with results from the empirical approach in "Parentis" and "Leon" stands.

3.3.2"Lit et Mixe" Plot

For the DSF plot in "Lit et Mixe", plotting ring area growth as a function of number of nests showed a triangular pattern as well (Fig. 9). Least square regressions (Fig. 9) was not significant (P = 0.1468) and explained just a very small part of variability (R² = 0.0008). The least square regression slope was -0.57.

The results of the different quantile regressions for the "Lit et Mixe" plot are shown in figure 9. The slope was always negative indicating a reduction in tree growth when the PPM infestation increased. The magnitude of the slope increased as higher quantiles data were used, up to the 0,995 th quantile (Fig. 9). The slope estimated for the 0.995 th quantile was -13.2 and had a 99% CI (confidence interval) that did not bracket zero, indicating a slope value significantly different from zero. For quantiles higher than 0.995 th , the regression slopes was also significantly different from zero.

Slope

For the "Lit et Mixe" plot, we got an average reduction of 12% in relative growth for the 6-25% defoliation class and an average reduction of 27% for the 26-60% defoliation class. As shown in Fig. 11, these results are consistent with the results from the empirical approach. Comparison with results from the empirical approach in "Parentis" and "Leon" stands.

Mont Ventoux plots

For the three plots from the Mt Ventoux area, data also revealed a triangular relationship between the ring growth area and the mean number of nest per tree.

Least square regressions (Fig. 12) were significant (P < 0.0001 for the three plots), but explained just a small part of variability (R² = 0.16, 0.23 and 0.20, respectively). Least square regressions slopes were -2.26, -3.27 and -2.4, respectively (Fig. 12).

The results of the different quantile regressions for plot 1 from the Mt Ventoux area are shown in Figure12. The magnitudes of slope remained negative and stable (Figure13) with higher quantiles. The slope estimated for the 0.9 th quantile was -2.55 and had a 99% CI (confidence interval) that did not bracket zero.

For quantiles higher than 0.9 th , the regression slopes was not significantly different from zero probably due to a low number of data within these quantiles.

Plot 2 and 3 from the Ventoux area showed similar patterns. Using a non parametric test of Friedman with the Forest Health Department data on Austrian pine, we found three significantly different numbers of PPM nests per tree for the three Slope defoliation classes (P < 0,0001). We could then calculate a mean number of nest/tree for each class of PPM defoliation in black pines: 1-5% defoliation corresponded to 0.20 ± 0,10 nests per tree 6-25% defoliation corresponded to 0.72± 0.27 nests per tree 26-60% defoliation corresponded to 1.76 ± 0.59 nests per tree Using the least square regressions on mean values, ring growth area values were calculated for each of these classes of PPM defoliation. We compared the reduction of ring growth area for these three plots with the results from the empirical approach. For the three plots, the growth loss averaged -8.5% for the 6-25% class, and -21 % for the 26-60% class.

As shown in the Figure 14, these results are consistent with the results from the empirical approach. The results from the plot 2 and 3 are very similar to the "Parentis" outcomes. Figure 14: Effect of PPM defoliation on relative growth loss compared to unattacked trees using the least square regression on mean values for "Ventoux" plots. Comparison with results from the empirical approach in "Parentis" and "Leon" stands.

Meta-analysis

Introduction

Here we present the results from a meta-analysis of 20 different case studies that compared a growth measurement between Thaumetopea sp. defoliated or control conifer trees. Metaanalysis is a powerful method to conduct quantitative, objective reviews of several independent studies (Jactel, 2007). Our objectives were (i) to determine if there was a significant effect of defoliation on tree growth, (ii) whether this effect significantly differed with the intensity of defoliation, (iii) and to test if the same relationship could be observed with other species of Thaumetopoea, namely T. bonjeani and wilkinsoni.

Method

• Data collection

We searched for suitable studies using on-line bibliographic databases (Web of Knowledge) and cited references in relevant publications. Our research yielded 20 case studies provided by 6 different publications published from 1977 to 2009. They involved three tree species (Pinus pinaster, Pinus brutia and Cedrus atlantica) and three insect species (Thaumetopoea pityocampa, bonjeani and wilkinsoni). Thaumetopoea bonjeani is present in North Africa and mainly feed on Cedrus. Thaumetopoea wilkinsoni is feeding on Pinus brutia in the eastern part of the Mediterranean Basin (e.g. Turkish region).

The criteria for the studies to be included in the meta-analysis were:

(i) tree growth, estimated with circumference variables (16 cases), volume variables (2 cases) or height variables (two cases), was compared between defoliated and unaffected (control) trees (ii) the mean of the growth response variable, a measure of its variance and the sample size was mentioned in the text or available via graphic digitization.

To avoid autocorrelation, results from different sampling dates were not use as replicates unless they were collected in different, independent plots. Data concerning independent groups of trees with different defoliation rates in the same study were compared to the same control group and considered as individual cases (i.e. a study comparing a control group of unattacked trees with two groups of attacked trees with defoliation rate of 25% and 100% respectively generated two individual case studies: 0 vs. 25% and 0 vs. 100%).

• Meta-analyses

Analyses were carried out using METAWIN 2.0 software. To compare the magnitude of an effect (here the effect of defoliation) across different independent studies, we calculated a standardised effect size. In our study we used the Response Ratio (eqn.1.1):

(1.1)

The Response Ratio was calculated as the logarithm of the ratio between mean measure of growth in defoliated trees ( DEF X , defoliated trees) and in control trees ( CTL X , unaffected trees). If DEF X is lower than CTL X , the ratio is lower than 1, and the log ratio is negative.

For each of this Response Ratio, a variance can be calculated. It depends on the sample size ( CTL N for the controlled trees and DEF N for the defoliated trees) and the standard deviation in the original study (respectively CTL s and DEF s for the control and the defoliated trees).

(1.2)

Then we combined all the studies to calculate a grand mean effect size as a weighted mean of E i = ln R i .. The weight attributed to each individual study was calculated as the inverse of

2 2 2 2 ln ) ( ) ( ) ( ) ( CTL CTL CTL DEF DEF DEF R X N s X N s v + = ) ln( ) ln( ) ln( ln CTL DEF CTL DEF X X X X R - = = variance ( i w = 1/ (vi)
) and therefore depended on the sample size and on the variance of initial growth measures, namely of the quality of the individual study. A lower weight was given to studies that were of lower experimental quality.

The equation of the grand mean effect size was as following: (eqn 1.3) (1.3) A bias -corrected bootstrap (9999 iterations) confidence interval was estimated for the mean. The mean was considered as statistically significant if its confidence interval did not bracket 0.

Results

Considering first the 17 case studies concerning Thaumetopoea pityocampa, our metaanalysis revealed a negative and significantly different from zero grand mean effect size of -0.32 (CI = -0,50 to -0,16), indicating a significant loss in tree growth Then we tested the effect of defoliation class:

-it was significant (P = 0.04) -it showed a response gradient (Table 4): the mean effect size was negative and significantly different from zero for the three higher defoliation classes [25-50%] [50-75%] and [75-100%] and its magnitude increased with the rate of defoliation (Fig 15). The growth loss was not significantly different from zero when the defoliation rate was lower than 25%. The relative growth loss was estimated by back transforming the mean effect size (lnR): Then we performed a second meta-analysis with all case studies (N= 20) including those concerning other Thaumetopoea species: Thaumetopoea wilkinsoni and bonjeani.

Relative Growth Loss = 1 -R= =1-℮ lnR = 1 -℮ Ei . ∑ ∑ = = + + = n i i n i i i w E w E
We obtained a grand mean effect size of -0.33 (CI = -0.50 to -0.16) then very similar to the value obtained with only T. pityocampa data.

In almost all of the cases, the individual effect sizes were negative (Fig. 16), indicating an overall growth loss in defoliated trees.

Only two positive effect sizes were observed which could be due 1) to imprecision in the measurement or 2) to the fact that dominant trees are more prone to T. pityocampa attacks; being more vigorous they could display higher growth than unattacked, smaller trees under low infestation levels.

The magnitude of growth loss due to T. bonjeani and wilkinsoni defoliation in Pinus and Cedrus were of similar magnitude than those observed in T. pityocampa defoliated trees. Green bars correspond to Thaumetopoea pityocampa, blue bars to Thaumetopoea wilkinsoni and the yellow bar to Thaumetopoea bonjeani.

General discussion

5.1. The three methodological approaches consistently showed that, overall, defoliation by processionary caterpillars resulted in a reduction of tree growth irrespective to the variable used to evaluate growth (diameter, height or volume). The Relative Growth Loss estimated with the three approaches was of the same order of magnitude, e.g. 25 -60% of loss for 26-60% of defoliation (Table 5). Similarly Gachi found that growth losses caused by Thaumetopoea bonjeani on Cedrus atlantica averaged 34% (M. Gachi and al., 2005). It was reported by Turkish researchers that following the damage by PPM (Thaumetopoea wilinsoni) on Pinus brutia, relative radial growth losses were of 12% [START_REF] Carus | Impact of defoliation by the pine processionary moth (Thaumetopoea pityocampa) on radial, height and volume growth of Carlabrian pine (Pinus brutia) trees in Turkey[END_REF]) and 21% [START_REF] Kanat | Thaumetopoeiadae) on annual diameter increment of Pinus brutia Ten[END_REF]. Cadahja and Insua (1970) found out relative losses in wood volume increment of 14% to 33% in young reforested area of P. radiata defoliated by Thaumetopoea pityocampa. [START_REF] Bouchon | Etude préliminaire sur les pertes de production des pinèdes soumises aux attaques de Thaumetopoea pityocampa Schiff[END_REF] reported that the Mt Ventoux forest of P. nigra periodically subjected to heavy attacks of T. pityocampa showed ca 45%-volume decrease during a 50 years period.

5.2.

There was a gradient of response with higher growth loss in more defoliated trees.

Through the three methodologies, we never found any asymptotic pattern. The growth loss still decreased until a complete defoliation.

The same gradient pattern was found in several papers. [START_REF] Markalas | Biomass Production of Pinus pinaster after defoliation by the Pine Processionary Moth (Thaumetopoea pityocampa) Other references[END_REF] found that one year after PPM defoliation on Pinus pinaster, losses in total shoot biomass were 41-50% in moderately defoliated trees and 54-64% in completely defoliated trees. [START_REF] Graf | L'infestation du cèdre de l'Atlas Cedrus atlantica Manetti par la chenille processionnaire du pin au Maroc et son impact sur la productivité[END_REF] observed losses of wood increment between 17 and 25% for the medium infestation rate (30 to 50% of defoliation) and 38 to 57% for high infestation rate (70 to 100% of defoliation) in Cedrus.

5.3. However the growth loss was not significantly different from zero in the lower class of defoliation, i.e. 1 -25%.

We found through the empirical approach and the meta-analysis that the relative growth loss was significant since a minimum threshold value.

In the empirical study, the growth in defoliated trees was significantly lower than the growth in control trees only for a defoliation rate higher than 25%. Based on the meta-analysis results, the growth loss was not different from zero for the 6 -25% defoliation class.

This threshold effect has been already reported in the literature. Only defoliation exceeding 25% caused significant losses in height and DBH (Diameter at Breast Height) of the pines (pinus pinaster), according to Barrento and al. (2008). Laurent-Hervouët (1986) also observed in pines (pinus nigra) that PPM defoliation in Corsica in 1974-1976 would not have any effect since the level of defoliation observed during that period was always lower than 50%.

There was no significant interaction between classes of defoliation classes and initial

tree diameter on the effect on tree growth.

With both the empirical and the dendrochronological approaches, we found parallel response patterns to PPM defoliation between large and small trees or between edge and interior trees. Even if the absolute growth loss was more severe in smaller trees than in larger trees, the relative growth loss was similar in the two types of trees.

Nevertheless the initial circumference of the trees had a significant effect on further tree growth. Moreover taller trees are known to be more prone to PPM infestation than smaller, dominated trees. This means that samples should involve trees of the same size and location in order to properly compare growth in defoliated and un-defoliated trees.

It took only 1-2 years for the trees to recover from one single defoliation.

We found that it took an average of one to two years after the main attack for the trees to recover using the empirical approach. Old trees from "Leon" stand may need more time to recover from the defoliation that the younger trees from "Parentis". [START_REF] Lemoine | Contribution à la mesure des pertes de production causées par la chenille processionnaire Thaumetopoea pityocampa Schiff au pin maritime dans les Landes de Gascogne[END_REF] estimated a recovery time of three years. Results from Gachi and al. (2005) showed that PPM defoliation slowed down tree growth during the following five years. Laurent-Hervouët (1986) found a recovery time of 1 -2 years. The recovery time may depend on specific pine species susceptibility as the three examples cited above concerned Pinus nigra, Pinus pinaster and Cedrus atlantica. However other factors may account for these discrepancies such as the recurrence of defoliation and local climatic and site conditions.

Future prospects:

The use of the dendrochronological methodology seems promising as it allows detecting past effects of insect defoliation on tree growth even if no individual tree survey has been conducted. However it would need further development to calibrate this approach and particularly to indicate which quantile should be selected in order to accurately estimate relative growth loss. Another avenue of improvement would be to reduce variance via incorporation of climatic data in multiple quantile regressions.

The pine processionary moth is spreading towards northern latitudes and higher altitudes due to the global warming. It would be then of great interest to develop a process-based model, based on knowledge on the pine processionary moth life cycle and nutrition on the hand, and on pine growth on the other hand. This eco-physiological model would then provide data for ecological studies of the relationships between PPM outbreaks, climate events, and stand management. It could be validated with the empirical, dendrochronological and metaanalytical data that we collected in this study. With such a model, we would be better able to predict and anticipate trends in defoliation by Thaumetopoea with the climate change.

The wind stability of different silvicultural systems for Douglas fir in the Netherlands: a model based approach

Mart-Jan Schelhaas, Centre for Ecosystem Studies, Alterra, Wageningen UR

Summary

The aim was to evaluate different silvicultural systems for Douglas fir (Pseudotsuga menziesii (Mirb.) Franco) in the Netherlands in terms of timber production and wind stability over a full rotation. This was done using the ForGEM-W model, which combines a distant dependent tree growth simulator with a mechanical-empirical wind damage module. Six different silvicultural systems were evaluated: normal yield table management, free thinning from above in a monospecies and a mixed stand (50% mixture of beech (Fagus sylvatica L.)), uneven-aged system, no thinning regime with low initial stand density of Douglas fir with and without admixture of beech. Silvicultural systems leading to low height-diameter (h/d) ratios were most successful in avoiding damage. Low h/d ratios were obtained in the system with low stand density and no thinning, and in the uneven-aged system by systematically removing trees with the highest ratios during thinning. Especially the uneven-aged system combined a high timber production with low risk. The use of Douglas-fir-beech mixtures changed the competition pressure on Douglas fir, and thus the h/d ratio and the wind risk. Results from this study indicate that the current trend towards more nature-oriented management could lead to lower wind risks and even to an increase in overall productivity.

Introduction

Wind damage is an important risk factor in forest management, especially in the temperate and boreal zone. In North-western Europe, the importance of wind is illustrated by a series of severe wind events during the last decades, interspersed with numerous smaller events [START_REF] Schelhaas | Natural disturbances in the European forests in the 19th and 20th centuries[END_REF]. In addition to wind speed, the vulnerability of trees and stands to wind damage are determined by many different factors such as tree and stand characteristics. The height and crown size of a tree largely determine how much force the wind exerts on the tree: the larger the tree and the crown, the higher the wind drag will be. Wind force will cause the stem to bend and exerts a force on the root-soil system. The resistance to stem breakage can be determined by the strength properties of the wood and by the diameter of the stem [START_REF] Peltola | A mechanistic model for assessing the risk of wind and snow damage to single trees and stands of Scots pine, Norway spruce, and birch[END_REF]. Correspondingly, the resistance to overturning depends on the shape and size of the root system, the soil type and on soil conditions (e.g. [START_REF] Coutts | Components of tree stability in Sitka Spruce on peaty gley soil[END_REF]). Stem diameter or stem mass has been found to correlate well with root mass and anchorage [START_REF] Nicoll | Anchorage of coniferous trees in relation to species, soil type, and rooting depth[END_REF][START_REF] Schelhaas | Introducing tree interactions in wind damage simulation[END_REF].

Trees having high height to stem diameter ratio at breast height (h/d ratio) for a fixed tree height require less wind speed to be damaged [START_REF] Peltola | A mechanistic model for assessing the risk of wind and snow damage to single trees and stands of Scots pine, Norway spruce, and birch[END_REF][START_REF] Wilson | Stability and density management in Douglas-fir plantations[END_REF][START_REF] Cameron | Importance of early selective thinning in the development of long-term stand stability and improved log quality: a review[END_REF]. With increasing height, h/d ratios should be kept lower to maintain the same risk level [START_REF] Cremer | Effects of stocking and thinning on wind damage in plantations[END_REF][START_REF] Becquey | L'existence de zones de stabilité des peuplements[END_REF][START_REF] Ruel | Understanding windthrow: silvicultural implications[END_REF]. However, forest stands with relatively high h/d values can still be stable if stand density is high enough, because of mutual support and sheltering amongst the individual trees [START_REF] Schelhaas | Introducing tree interactions in wind damage simulation[END_REF].

Management decisions and silvicultural treatments influence the state of the forest and thus its susceptibility to windthrow. For example, clear fellings suddenly expose adjacent stands to the wind [START_REF] Zeng | Influence of clear-cutting on the risk of wind damage at forest edges: A gis-based integrated models approach[END_REF]. Thinnings are known to reduce the stand stability for several years [START_REF] Cremer | Effects of stocking and thinning on wind damage in plantations[END_REF][START_REF] Savill | Silviculture in windy climates[END_REF][START_REF] Ruel | Understanding windthrow: silvicultural implications[END_REF][START_REF] Cameron | Importance of early selective thinning in the development of long-term stand stability and improved log quality: a review[END_REF], but can lead in the longer term to lower average h/d ratios. In highly exposed sites it is recommended to use a relatively sparse initial stand density and not to thin [START_REF] Quine | Forests and wind: Management to minimize damage[END_REF][START_REF] Ruel | Understanding windthrow: silvicultural implications[END_REF]. As low initial stand densities will generally negatively affect the wood quality [START_REF] Deans | Effects of respacing on young Sitka spruce crops[END_REF][START_REF] Fahlvik | Influence of precommercial thinning grade on branch diameter and crown ratio in Pinus sylvestris in southern Sweden[END_REF], the use of a self-thinning mixture is proposed as alternative [START_REF] Quine | Forests and wind: Management to minimize damage[END_REF][START_REF] Ruel | Understanding windthrow: silvicultural implications[END_REF].

Tree species choice is also an important management decision with regard to vulnerability to wind damage. Species differ in wind vulnerability among others due to differences in wood properties, rooting strategy and crown shape. Broadleaves are generally considered windfirm [START_REF] Holmsgaard | Historical development of wind damage in conifers in Denmark In: Minimizing wind damage to coniferous stands[END_REF][START_REF] Quine | Forests and wind: Management to minimize damage[END_REF] whereas among the conifers especially Norway spruce (Picea abies (L.) Karst.) and Sitka spruce (Picea sitchensis (Bong.) Carr.) are seen as unstable species [START_REF] Slodicák | Thinning regime in stands of Norway spruce subjected to snow and wind damage[END_REF]. Several authors [START_REF] Eriksson | Windthrow damage in forests in relation to stand treatment-present state of knowledge in Sweden[END_REF][START_REF] Slodicák | Thinning regime in stands of Norway spruce subjected to snow and wind damage[END_REF][START_REF] Schütz | Vulnerability of spruce (Picea abies) and beech (Fagus sylvatica) forest stands to storms and consequences for silviculture[END_REF] advocated that introducing more windfirm species into a stand of less windfirm species would yield a higher stability of the total stand (see [START_REF] Lüpke | Aspects of stability, growth and natural regeneration in mixed Norway spruce-beech stands as a basis of silvcultural decisions[END_REF] for a literature overview), although others suggest that this will result only in the loss of the unstable component [START_REF] Lüpke | Aspekte der Stabilität und des Wachstums von Mischbeständen aus Fichte und Buche als Grundlage für waldbauliche Entscheidunge[END_REF] or that the combination could become even more unstable [START_REF] Savill | Silviculture in windy climates[END_REF]. As stabilising species, usually broadleaves are recommended [START_REF] Eriksson | Windthrow damage in forests in relation to stand treatment-present state of knowledge in Sweden[END_REF].

Despite a vast body of literature on the relation between vulnerability of trees to wind, forest characteristics and management, forest managers still express the need for more support. [START_REF] Blennow | Risk perception among non-industrial private forest owners[END_REF] showed that forest owners in Southern Sweden rank the risk of wind damage highly, but they generally do not know how to change their forest management to reduce the risk of damage. Mathematical models provide an opportunity for objectively calculating these risks and can show forest managers the long-term implications of their actions [START_REF] Gardiner | Management of forests to reduce the risk of abiotic damage -a review with particular reference to the effects of strong winds[END_REF].

In recent years, substantial progress has been made in the modelling of windthrow. [START_REF] Peltola | A mechanistic model for assessing the risk of wind and snow damage to single trees and stands of Scots pine, Norway spruce, and birch[END_REF] and Gardiner et al. (2000) have developed mechanistic models to estimate, for a given tree, the critical wind speed needed for breakage or uprooting. In these models, the stand under consideration is assumed to be even-aged and comprised of a single species. Furthermore, it is assumed to consist of an evenly-spaced array of trees with the same height and diameter (the socalled "average tree"). [START_REF] Ancelin | Development of an individual tree-based mechanical model to predict wind damage within forest stands[END_REF] extended this approach and evaluated all trees in a particular stand. [START_REF] Schelhaas | Introducing tree interactions in wind damage simulation[END_REF] introduced a similar mechanistic module into a single tree model and improved shelter and support mechanisms between trees in Douglas fir (Pseudotsuga menziesii (Mirb.) Franco). This integrated model, ForGEM-W, is able to simulate a wide variety of silvicultural practices under a specified wind climate.

One of the countries in North-western Europe that is affected regularly by storms is the Netherlands. Major wind damage to the forest occurred in 1972/73 and 1990 [START_REF] Schelhaas | Database on Forest Disturbances in Europe (DFDE)[END_REF], with wind-felled timber between 0.6 and 1.5 million m3. Lesser damages occurred in 1976 and 2007 with respectively 0.15 and 0.25 million m3 [START_REF] Schelhaas | Database on Forest Disturbances in Europe (DFDE)[END_REF]. Douglas fir was one of the species reported to be vulnerable in these and other events ( [START_REF] Soest | Stormschade aan Douglas[END_REF][START_REF] Faber | Stabiliteit van bos ten opzichte van wind. Een theoretisch gezichtspunt[END_REF][START_REF] Sissingh | Stabiliteit van bos ten opzichte van wind en storm gezien vanuit de praktjik[END_REF].

Although the share of Douglas fir in the total forest area is only about 6% [START_REF] Dirkse | Meetnet Functievervulling Bos 2001-2005[END_REF], it is considered an important species due to its high productivity and good wood quality. Moreover, it regenerates well naturally and is expected to increase its share in future.

Traditionally it was managed in monocultures with thinning from below. Nowadays the future tree method and thinning from above are common ways of managing Douglas fir. In the future tree method, about 60-100 trees are selected that represent the most valuable trees and that are expected to form the final crop. Selection usually takes place at an age of about 40. During the further life of the stand, these trees are consistently favoured in thinnings. As Douglas fir is a rather shade-tolerant species, it can also be grown in uneven-aged and selection systems. In regeneration groups it can be mixed with fast growing species as larch (Larix sp.) and birch (Betula sp.), but eventually these species will be overgrown by the Douglas fir. Similarly, Douglas fir is known to out compete also beech (Fagus sylvatica L.) in individual mixtures, and therefore, [START_REF] Hekhuis | De menging douglas -beuk : onderzoek naar de individuele gelijkjarige menging van douglas en beuk, met douglas als hoofdboomsoort, op holtpodzolen in Nederland[END_REF] recommended mixing both species in strips of 10-15 m or in groups. Despite the importance of Douglas fir and its reported wind vulnerability, not much is known yet about the consequences of different silvicultural systems on its vulnerability to wind damage, for example, in the Netherlands.

The aim of this study was to evaluate the effectiveness of a set of silvicultural systems with regard to timber production and wind stability by means of ForGEM-W simulations. Silvicultural systems were selected to represent the whole range of possible systems that would likely be applied to Douglas-fir stands in the Netherlands. They were evaluated over a full rotation period to enable a good comparison. The focus was on sheltered stands. However, although trees are known to adapt to more exposed conditions in permanent stand edges [START_REF] Cucchi | Root anchorage of inner and edge trees in stands of Maritime pine (Pinus pinaster Ait.) growing in different podzolic soil conditions[END_REF][START_REF] Gardiner | The stability of different silvicultural systems: a wind-tunnel investigation[END_REF], underlying mechanisms are still insufficiently understood and not included in the ForGEM-W model.

Methods

Outline of ForGEM-W

ForGEM

ForGEM (Forest Genetics, Ecology and Management) is a physiologically-based model that simulates the growth and development of individual trees, on a scale up to several hectares [START_REF] Kramer | Effects of silvicultural regimes on dynamics of genetic and ecological diversity of European beech forests; impact assessment and recommendations for sustainable forestry[END_REF]. Several processes can be simulated with different levels of detail according to the needs of the user. In the current application, the simplest simulation methods were used. This means simulation of light interception based on the gap-type approach and no simulation of water and nitrogen balances. Here only the basics of the model are explained, a more detailed description can be found in [START_REF] Kramer | Effects of silvicultural regimes on dynamics of genetic and ecological diversity of European beech forests; impact assessment and recommendations for sustainable forestry[END_REF] and [START_REF] Kramer | Bridging ecophysiological and genetic knowledge to assess the adaptive potential of European beech[END_REF]. Most approaches are essentially similar to those used in SORTIE [START_REF] Pacala | Forest models defined by field measurements: 1. The design of a northeastern forest simulator[END_REF][START_REF] Pacala | Forest models defined by field measurements: II. Estimation, error analysis and dynamics[END_REF], except in this case the gap-type approach for light interception.

ForGEM simulates individual trees that have known coordinates. Tree growth is driven by light interception. In the gap-type approach, light is assumed to come straight from above. The light interception in each 20x20 m grid cell is divided amongst the trees in the cell, according to their foliage mass and its vertical distribution (see [START_REF] Bugmann | A review of forest gap models[END_REF]. Intercepted radiation is converted into photosynthates via the radiation use efficiency (RUE) parameter [START_REF] Landsberg | A generalised model of forest productivity using simplified concepts of radiation-use efficiency, carbon balance and partitioning[END_REF]. The RUE must be calibrated against independent data, such as growth and yield tables.

Photosynthates are allocated in such a way that specific ratios between different tree components (stem mass, branches, root mass, etcetera) are maintained (see [START_REF] Kramer | Bridging ecophysiological and genetic knowledge to assess the adaptive potential of European beech[END_REF]. These ratios depend on the total aboveground biomass. Photosynthates allocated to the stem are converted into height and stem diameter increment. Maximum height increment follows a Richards' growth curve [START_REF] Richards | A flexible growth function for empirical use[END_REF]Jansen et al., 1996). Actual height increment can be lower due to insufficient resources. The relationship between tree volume, height and stem diameter at breast height (DBH) is described by species-specific allometric functions, derived from yield table data (Jansen et al., 1996). Crown expansion is influenced by overlap with other tree crowns [START_REF] Schelhaas | Introducing tree interactions in wind damage simulation[END_REF].

Competition-induced mortality is assumed to occur if the crown of a tree is completely suppressed by those of its neighbours (see [START_REF] Schelhaas | Introducing tree interactions in wind damage simulation[END_REF] for details). Age-related mortality is assumed to follow a 2-parameter Weibull distribution. Mortality is calculated at the end of each month. Dead trees are removed from the tree list and are assumed to be left in the forest. Regeneration is explicitly simulated by seed production, dispersal and germination.

Seedlings are treated as cohorts in 5x5 m grid cells. Self thinning is simulated in these cohorts using the -3/2 power law (e.g. [START_REF] Reineke | Perfecting a stand density index for even-aged forests[END_REF][START_REF] Drew | Stand density management: an alternative approach and its application to Douglas-fir plantations[END_REF]. Upon reaching 2 m height they are further treated as individuals and receive x,y-coordinates. Several processes have a stochastic character (mortality, weather), so simulations should be repeated several times.

A wide range of management options is available in ForGEM. These include tree planting, tending, thinning, cutting gaps, clear-cut, shelterwood, a future tree selection system, and thinning aimed at a target diameter distribution. Thinnings can be specified according to agedensity schemes, removal of a certain fraction of the basal area increment or according to the available crown space. Thinnings can be done from below, from above or randomly.

Wind damage module

The wind damage module uses a static mechanistic approach, largely following the principles of the HWIND model [START_REF] Peltola | A mechanistic model for assessing the risk of wind and snow damage to single trees and stands of Scots pine, Norway spruce, and birch[END_REF]Gardiner et al., 2000). Based on characteristics of the tree (height, DBH, crown size), but also its surroundings, it calculates which trees will be broken or uprooted for a given mean hourly wind speed. Wind drag on a tree is calculated according to [START_REF] Monteith | Principles of environmental physics[END_REF], where the crown and stem shape determine the effective sail area of a tree, adjusted for streamlining of the crown (see [START_REF] Peltola | A mechanistic model for calculating windthrow and stem breakage of Scots pines at stand edge[END_REF]. The weight of the displaced stem and crown add to the total turning moment at stem base.

However, unlike in HWIND, this turning moment is modified according to the shelter and support received by surrounding trees. Support is assumed to be linearly related to the crown overlap of trees. For each tree a shelter zone is defined, consisting of a strip equal to the crown width and reaching until the first upwind edge. This can be either the forest edge or a gap inside the forest. A gap is defined as an area free of trees higher than the crown base of the subject tree, with a minimum gap length of half the tree height of the subject tree. Sheltering effects depend on the stand density in this zone, modified for relative heights of trees and the presence of foliage. The relationship between stand density and sheltering was adapted from wind tunnel studies in even-aged stands [START_REF] Gardiner | Field and wind tunnel assessments of the implications of respacing and thinning for tree stability[END_REF], as described in [START_REF] Peltola | A mechanistic model for assessing the risk of wind and snow damage to single trees and stands of Scots pine, Norway spruce, and birch[END_REF] and Gardiner et al. (2000).

Trees are assumed to break or uproot if respectively the maximum stem resistance or the maximum anchorage resistance is exceeded. Stem resistance is a function of MOR and DBH to the power 3 (see [START_REF] Peltola | A mechanistic model for assessing the risk of wind and snow damage to single trees and stands of Scots pine, Norway spruce, and birch[END_REF].

Resistance to uprooting depends also on DBH according to:

2 1 a RES DBH a ROOT × =
where a 1 and a 2 are species specific parameters, estimated from tree pulling tests [START_REF] Schooten | Fysische aspecten van windworp bij de douglasspar, Vakgroep Bosteelt[END_REF]. Trees can experience additional loading if hit by fallen trees (the domino effect). Trees smaller than 5 m cannot be uprooted or broken, but they can be destroyed by falling trees. Trees are removed from the tree list at the end of each monthly time step. More details on the wind damage module can be found in [START_REF] Schelhaas | Introducing tree interactions in wind damage simulation[END_REF].

Model computations

Selected silvicultural systems

The following six silvicultural systems were selected to be tested, providing a broad range of likely systems:

1. Normal yield table system. Every five years, the number of stems was reduced to the values indicated in the yield table for normal stand density from Jansen et al. (1996). A thinning from below was applied, i.e. trees to be removed were those with the lowest stem diameters, compatible with Jansen et al. (1996). In figures, this system is further referred to as "yield table". 2. Free thinning from above. The implementation of this silvicultural system was based on the management guidelines for Douglas fir forests from the Dutch State Forest Service. Every five years the best trees were selected and liberated from their immediate competitors.

Competitors were defined as trees that have an overlapping crown projection and that reached up to at least half the crown of the tree of interest. The number of trees to be liberated decreased stepwise from 200 at age 20 to 60 at age 80. Preferred maximum h/d ratios for these trees decrease from 100 at age 20 to 70 at ages above 40. This system is further referred to as "free thin mono". 3. Free thinning from above with mixture of more windfirm beech. The same silvicultural system as described for the "high mono" was used, but 50% of the selected trees to be liberated had to be beech. This system is further referred to as "free thin mix". 4. Uneven-aged system. This system aims to maintain an exponentially decreasing diameter distribution:

q j j N N 1 + =
where N j is the number of trees in diameter class j and q a constant. Q was 1.3 and the target diameter 60 cm (De [START_REF] Klein | Planning en controle in ongelijkjarig bos[END_REF]. Thinning starts from the highest diameter class. The required number of trees was selected with preferably low h/d ratios and liberated as described before. Excessive trees within the diameter class were removed. Subsequently the lower diameter classes were evaluated for thinning, taking into account the already removed trees. This system is further referred to as "uneven". 5. No thinning. A no-thin system with an initial stand density of 200 trees per hectare was selected, which is the final density reached by the normal yield table management. This system is further referred to as "no thin mono". 6. Self-thinning mixture. A mixture of 200 Douglas firs and 3800 beech trees per hectare were used, where the beech was expected to be outcompeted by the Douglas fir. This system is further referred to as "self-thin mix".

Simulations and comparisons of silvicultural systems

All simulated systems were studied on a square forest plot of 1 ha. On the north, west and south side a buffer strip of 50 m was assumed to provide shelter to the experimental plot, since only sheltered stands were dealt with. This buffer strip consisted of exactly the same forest as the simulated stand and received the same treatment. However, no wind damage was allowed in the buffer. A buffer strip to the east was deemed unnecessary because high wind speeds are always connected to south to north-westerly directions. Still, in the simulations, 5-10% of the volume damaged was due to easterly winds. However, omitting the easterly buffer strip considerably reduced the computing time needed. The whole forest plot was assumed to be surrounded by an open strip of 100 m, needed to calculate the wind profile outside the forest. The site class was assumed to be 16 for Douglas fir and 8 for beech (Jansen et al. (1996), referring to the maximum mean annual volume increment over a rotation that can be reached for that species on that site).

A young 3 ha stand was generated as initial situation for all even-aged mono-species systems. Initial conditions were taken from Jansen et al. (1996) for a stand with normal density (4000 stems/ha) at age 10, assuming a regular spacing. For the free thinning from above with mixture of beech (system 3) the same initial situation was used, but assuming that all trees in every second 12.5 m strip were beech. Strips were located north-south, providing most shelter to the predominant westerly winds. For the self thinning mixture (system 6) all trees were assumed to be beech, except 200 trees per ha located in a regular grid. To create an uneven-aged initial stand, gaps were cut randomly in the regular stand every 10 years for a period of 100 years.

Thereafter uneven-aged management was applied until the diameter distribution was more or less in balance. For the even-aged simulations, no regeneration was allowed to occur in order to save computing time.

Model parameters for Douglas fir were the same as those used by [START_REF] Schelhaas | Introducing tree interactions in wind damage simulation[END_REF], except that the maximum crown radius was decreased from 7.7 to 6.5 m. This was done because of excessive mortality at older ages compared with yield table values (Jansen et al., 1996). Consequently the RUE had to be re-calibrated, as well as the parameters that govern the decrease of RUE over age (T and c, Table 1, see [START_REF] Schelhaas | Introducing tree interactions in wind damage simulation[END_REF]. Beech was chosen as an example of wind stable species for admixture. For purposes of demonstration and not to complicate the analysis further, this species was assumed not to be damaged by wind at all, so not all parameters were needed for this species. Relevant parameters are shown in Table 1, others are equal to those in [START_REF] Kramer | Bridging ecophysiological and genetic knowledge to assess the adaptive potential of European beech[END_REF]. Furthermore, a 10% variation around the mean was introduced in the height growth parameter for individual trees. This was done to ensure that trees emerging from the same cohort would develop in slightly different ways, allowing for processes as competition and selection to work well. Meteorological variables such as daily radiation, temperature and precipitation are generated by a weather generator for the growth simulations. A separate generator was added to provide daily maximum mean hourly wind speed data (further referred to as wind speed). This generator is based on the Rijkoort-Weibull model for mean daily wind speed [START_REF] Rijkoort | A compound Weibull model for description of surface wind velocity distributions[END_REF] and can be parameterised using measured wind speed series. For a simulation, a year was divided in six seasons of two months duration, starting with January/February. Each season has its own wind direction frequencies. Each of the combinations of the eight wind directions and six seasons has its own Weibull distribution. In the simulation, a wind direction is randomly selected each day according to the frequency distribution for that season, and a wind speed is also randomly selected based on the corresponding Weibull distribution.

Wind speed data for the period April 1961 -December 2005 were obtained from the Dutch Royal Meteorological Institute (KNMI, 2007) for the measurement station Leeuwarden (53°13' N, 5°46' E). Measurements were already corrected for roughness of the surrounding terrain. Figure 1 gives an impression of the wind climate with regard to wind direction and wind speed. The highest wind speed measured at Leeuwarden was 28.1 m/s. Based on this data, the daily maximum mean hourly wind speed and the corresponding direction were extracted. Data were grouped per two-month season and eight wind directions (north between 337.5 and 22.5 degrees, etc). For each of the 48 combinations of season and wind direction a 2-parameter Weibull distribution was fitted, using the maximum likelihood method. Furthermore, seasonal direction frequencies were derived. Using the weather generator, 10 time series of daily weather were generated, each covering a period of 100 years. All silvicultural systems were tested for these 10 weather scenarios and all results presented are the means of these 10 replicates. Each of the silvicultural systems had its own dynamics with respect to the development of average diameter, standing volume, harvest and other stand characteristics. Considering the differences in diameter development, it was not realistic to compare the systems over the same time span of 100 years. Therefore in all cases a target diameter of 60 cm was assumed, equal to the maximum diameter in the uneven-aged system (De [START_REF] Klein | Planning en controle in ongelijkjarig bos[END_REF]). In the evenaged cases, the analysis was limited to the period needed to obtain an average (quadratic) diameter of 60 cm of the standing Douglas fir trees. In case this diameter was not reached during the simulation period, the full 100 year period was analysed. In the uneven-aged case a 100 years period was analysed, since no information was available on the average time required for a single tree to reach the target diameter.

All systems were compared visually for the development of h/d ratio against height, where for the uneven-aged system the average height and h/d ratio per 2.5 m height class for a randomly selected year was used as comparison. Even-aged systems were also visually compared with regard to development of stand density, DBH, height and h/d ratio over time.

For all silvicultural systems, the gross annual volume increment (GAI) over the period under consideration was assessed, calculated as following:

i T T t it T t it T t it i T S S H M W GAI i i i i / )) ( ( 0 0 0 0 - + + + = ∑ ∑ ∑ = = =
where GAI i is the gross annual increment (m 3 ha -1 yr -1 ) for system i, W it , M it , H it respectively the wind damage, the other mortality and the harvest in system i in year t (all in m 3 ha -1 ), T i the evaluated timespan for system i and S Ti and S 0 the standing wood volume at the end and the start of the simulation. It was assumed that all dead timber is left in the forest, both due to mortality and wind damage. Total timber production is then calculated as the components harvest and standing wood volume increase in above formula. All silvicultural systems were also visually compared in terms of gross annual increment and the distribution of its components.

Results

Stand dynamics development

The development of the mean diameter over time differed between systems (Figure 2) and was related to the differences in stand density development. The normal yield table system showed a very regular development of diameter, caused by a very regular decrease in stand density over time. After 100 years of simulation a diameter of 58.7 cm was reached. The system free thinning from above mono maintained a higher stand density in the first half of the simulation, leading to lower mean diameters. Later stand densities were lower, leading to a more rapid increase in diameter. After 100 years a diameter of 59.9 cm was reached.

In the system free thinning from above mixed with beech, total stand density equalled that of the free thinning from above mono system and was higher later on in the simulation. However, the share of Douglas fir in the stand density decreased quickly to around 20%. The mean diameter of the Douglas fir increased rapidly, and reached the 60 cm target already after 69 years. Diameter growth was the fastest in the system no thinning mono, but in the second half of the simulation the speed of growth decreased. On the contrary, the system self thinning mixture showed a slower diameter growth to start with, but a higher growth rate in the second half of the simulation. As a result, both systems reached the 60 cm diameter at about the same time: after 85 years for the no thinning mono and 82 years for the self thinning mixture system. Development of average height was almost equal in all systems, except for the system free thinning from above mono where average height was 4-6 m lower during most of the simulation.

Also the development of the h/d ratio showed large differences between systems. In the normal yield table system the h/d ratio varied between 80 and 90 most of the time and dropped only below 80 when trees were over 35 m high. In the system free thinning from above mono, the h/d ratio increased to over 120 at a tree height of 15 m, but decreased steadily to 60. The system free thinning from above mixed followed a similar pattern, but reached at maximum only a ratio of 114. In the free thinning from above mono system, the h/d ratio varied mostly between 50 and 60, whereas the system self thinning mixture showed a gradual decline from 80 to 60. In the uneven-aged system, average h/d ratios of small trees were high, with values close to 120 for trees until 15 m. However, h/d ratios decreased rapidly for higher trees, with a value of 55 for the tallest trees in the stand.

Timber production and intensity of wind damage

Systems differed substantially in development of growing stock, harvest, mortality, wind damage and gross annual increment (Figure 3). In the normal yield table system, wind damage started to be substantial after about 60 years, leading to a decrease in growing stock level from 470 m 3 ha -1 to 200 m 3 ha -1 in the end. Also after 50 years some mortality started to occur. The total gross annual increment was 11.8 m 3 ha -1 yr -1 , with standing volume and harvest accounting for just over half of it (Figure 4). In the system free thinning from above mono, wind damage started to occur already after about 40 years, but increased less quickly over time. Mortality occurred especially between 20 and 40 years. The gross annual increment was 12.5 m 3 ha -1 yr -1 . Developments in the system free thinning from above mixed were largely comparable, but the target diameter was reached much earlier. The gross annual increment was 13.8 m 3 ha -1 yr -1 , with harvest and standing volume accounting for about 75%. The system no thinning mono showed hardly any wind damage and only little mortality. However, gross annual increment was substantially lower with 8.5 m 3 ha -1 yr -1 . In the system self thinning mixture the gross annual increment was 8.9 m 3 ha -1 yr -1 , with about half of it in harvest and standing volume. Wind damage started to occur after 55 years. The uneven-aged system had the highest gross annual increment with 18.0 m 3 ha -1 yr -1 . This system showed hardly any wind damage and mortality.

There was no clear relationship between gross annual increment (Figure 5) and the amount of wind damage. Low level of damage was found both in the no thinning mono and uneven-aged system, but gross annual increment was nearly double in the uneven-aged case. Similarly, the systems free thinning from above mono, free thinning from above mixed and normal yield table had comparable gross annual increment, but the amount of wind damage was more than twice as high in the normal yield table system as compared to the system free thinning from above mixed. Within the 10 replicates per each system, gross annual increment showed not much variation, whereas wind damage was much more variable (Figure 5). Systems with higher average wind damage levels showed more variation than systems with low average wind damage levels, too. Furthermore, the systems differed in the lowest wind speed that caused damage. In the uneven-aged system, damage occurred already at wind speeds of 7 m/s. In the system free thinning from above mono, wind damage started at 9 m/s, while in the systems free thinning from above mix and normal yield table this was at approximately 11 m/s. The systems no thinning mono and self-82 thinning mixture showed the highest wind speed needed to cause damage with 16 m/s. The highest wind speed present in the simulation was 26.9 m/s.

Discussion and conclusions

4.1 Evaluation of the findings Differences between systems with regard to stand dynamics can be explained by inter and intra species competition. The system free thinning from above mono had the highest stand density early in the simulation and, thus, the highest competition. This resulted in the lowest diameter growth of all systems, in retarded height growth and in considerable mortality early in the simulation. The normal yield table system had a lower stand density than the system free thinning from above mono until 40 years and thus a higher diameter growth and virtually no mortality. At higher ages the stand density was higher than the system free thinning from above mono, leading to slower diameter growth and more mortality. The system no thinning mono had the lowest stand density, the highest diameter growth and no mortality up to about 50 years. After that, competition increased, causing lower diameter growth and start of mortality. In the system self thinning mixture, 3800 beech trees per hectare were planted in between the sparse (200 stems ha -1 ) Douglas-fir stand. This resulted in increased competition, as evidenced by the lower diameter growth of Douglas fir as compared to the system no thinning mono. However, at a later stage diameter growth recovered as the Douglas fir outcompeted the beech. In the system free thinning from above mixed, half of the Douglas fir were replaced by beech. As Douglas fir is a better competitor, its diameter growth increased as compared to the system free thinning from above mono. Both mixed cases showed that Douglas fir is out competing beech, as expected. Furthermore, reactions of stand characteristics to changes in stand density and tree species composition showed realistic and consistent patterns.

Since height growth was more or less equal among the even-aged systems (except the system free thinning from above mono), the development of the h/d ratio depended on the diameter development and, thus, on the stand density development. Management only had an indirect influence on the h/d ratio through manipulation of stand density. However, in the uneven-aged case management also had a direct influence on the h/d ratio, because h/d ratio was the main selection criterion when the stand density had to be decreased in a certain diameter class. By removing the trees with the highest h/d ratio, the average ratio in a diameter class was decreased after each thinning. Since the number of trees to be retained decreases with increasing 84 diameter class, there was a strong selection on trees with a low h/d ratio. This explains the rapid decrease of h/d ratio with increasing height (Figure 2). Moreover, remaining trees profit from decreased competition and will grow quickly towards the target diameter.

Not all silvicultural systems utilised the growing space equally efficient. Gross annual increment ranged from 8.5 in the system no thinning mono to 18.0 m 3 ha -1 yr -1 in the uneven-aged case (Figure 4). Lower increments appeared when stand density deviated from the optimal stand density. Too dense stands, as for example in the first half of the system free thinning from above mono, causes high competition, with relatively little investment in height and diameter growth of the stem. Too open stands, as for example the start of the system no thinning mono, do not utilise all incoming light. Furthermore, part of the difference between even-aged and uneven-aged stands was caused by the fact that no regeneration was allowed in the evenaged systems. Especially in the systems with low stand densities a second tree layer could considerably increase the overall productivity.

Differences in wind damage level between systems can be linked to differences in development of the h/d ratio. The longer a tree remains in a situation with relatively high h/d ratios, the larger the chance that it actually will be damaged. At tree heights above 25 m, the system normal yield table showed the highest h/d ratios of all silvicultural systems (Figure 2). Moreover, this situation lasted quite some time, increasing the risk that critical wind speeds will be reached. This system also showed the highest damage amount (Figure 4). The system free thinning from above mono showed the highest absolute h/d ratios reached in any system and had a relatively high damage level as well. In the system free thinning from above mixed, h/d ratios were somewhat lower, and so was also the damage level. However, this was probably also related to a shorter simulation period, and thus the time trees were exposed to risk of wind damage. The system no thinning mono had hardly any wind damage, connected to very low h/d ratios throughout the rotation and a short time needed to attain the target diameter. Admixture of beech (the system self thinning mixture) resulted in higher h/d ratios and a higher damage level in Douglas fir, although the time to reach the target diameter was about the same.

However, the combination of h/d ratio and height alone are not sufficient to predict the damage level. The systems self thinning mixture and free thinning from above mono had about the same level of damage, but h/d ratios in the high mono system were much higher. Higher h/d ratios could, to a certain extent, be compensated for by more shelter and support of 85 surrounding trees as was suggested also previously by [START_REF] Ruel | Understanding windthrow: silvicultural implications[END_REF], for example. This was also demonstrated by the uneven-aged case where initial h/d ratios were rather high, but overall damage was small. This is also caused by the fact that overstory trees will grow quickly to the desired diameter due to the large growing space. Thus, the exposure time is relatively limited.

Several authors state that mixed stands are more stable than mono species stands (e.g. Slodicak, 1995;[START_REF] Schütz | Vulnerability of spruce (Picea abies) and beech (Fagus sylvatica) forest stands to storms and consequences for silviculture[END_REF]. In the case where half of the Douglas fir was replaced by beech (system free thinning from above mixed vs. free thinning from above mono), damage was indeed reduced by 50%, in line with findings of [START_REF] Lüpke | Aspekte der Stabilität und des Wachstums von Mischbeständen aus Fichte und Buche als Grundlage für waldbauliche Entscheidunge[END_REF]. However, this contradicts the study of [START_REF] Schütz | Vulnerability of spruce (Picea abies) and beech (Fagus sylvatica) forest stands to storms and consequences for silviculture[END_REF] who found that an admixture of 10% wind-firm tree species in Norway spruce stands lead to a decrease in damage by a factor of 3.4. These differences might perhaps be explained by differences in stand history, or preference of mixed stands on certain sites, like sites with lower fertility. Furthermore, use of beech as mixture in the no thinning mono system lead to a considerable increase in damage. This was caused by the effect beech has on the h/d ratio of Douglas fir through increased competition. This possible effect has also been hypothesised by [START_REF] Cameron | Importance of early selective thinning in the development of long-term stand stability and improved log quality: a review[END_REF]. In the current simulations, the effect of mixtures was limited to aboveground competition and sheltering effects. However, belowground competition is an important factor as it can affect the development of deeper or shallower root systems [START_REF] Hendriks | Root density and root biomass in pure and mixed forest stands of Douglas-fir and Beech[END_REF][START_REF] Schmid | Vertical distribution and radial growth of coarse roots in pure and mixed stands of Fagus sylvatica and Picea abies[END_REF], but may also result in loss of interlocking between trees of the same species [START_REF] Elie | Windthrow hazard modelling in boreal forests of black spruce and jack pine[END_REF].

The uneven-aged system was found to suffer less damage in comparison with normally practised even-aged systems. This was in line with observations [START_REF] Dvorak | Sturmschäden in ungleichförmigen Beständen[END_REF][START_REF] Schütz | Vulnerability of spruce (Picea abies) and beech (Fagus sylvatica) forest stands to storms and consequences for silviculture[END_REF]. Many authors attribute the greater stability of uneven-aged systems to the fact that overstorey trees will become accustomed to the wind [START_REF] Cameron | Importance of early selective thinning in the development of long-term stand stability and improved log quality: a review[END_REF][START_REF] Mason | Are irregular stands more windfirm?[END_REF]. However, such a mechanism was not present in the current model. Instead, the greater stability of the larger trees was caused by high diameter growth of these trees (due to a higher availability of light, nutrients and water) and a strong selection on low h/d ratios. The observed trend for larger trees in this system to have more favourable h/d ratios was confirmed by observations from [START_REF] Kenk | Management of transformation in central Europe[END_REF], but this might at least be partly caused by acclimation to wind. 86

Model aspects

Especially DBH was strongly influenced in simulations by the silvicultural system. But, still developments of DBH and thus, h/d ratios showed realistic and consistent patterns between the systems. However, it was not possible to validate the absolute values, other than for the normal yield table system.

Although the model was calibrated on the normal yield table, differences occurred. At age 100, simulated height and DBH were 38.3 m and 53.0 cm, respectively, compared to yield table values of 37.1 m and 51.4 cm. Differences were most likely caused by the inclusion of wind damage events in the simulations, leading to lower stand densities (i.e. 95 simulated trees per hectare against 206 according to the yield table). However, only the last 10-15 years of the simulation densities were considerably lower than the yield table, so the differences in height and DBH are still rather small.

A larger problem is the sensitivity of wind damage module to anchorage parameters. The underlying sample size was rather small and for trees larger than 40 cm DBH anchorage had to be extrapolated [START_REF] Schelhaas | Introducing tree interactions in wind damage simulation[END_REF].

Although other tree pulling studies found no changes of relationships between tree size and anchorage at diameters above 40 cm [START_REF] Papesch | Mechanical stability of Pinus radiata trees at Eyrewell forest investigated using static tests[END_REF][START_REF] Lundström | Anchorage of mature conifers: resistive turning moment, root-soil plate geometry and root growth orientation[END_REF], more experimental work is needed to confirm this. The shelter function of the wind damage module was also adapted from wind tunnel studies in even-aged stands [START_REF] Gardiner | Field and wind tunnel assessments of the implications of respacing and thinning for tree stability[END_REF]. In the future, these functions should be refined using wind tunnel studies in irregular stands, as was done for example by [START_REF] Gardiner | The stability of different silvicultural systems: a wind-tunnel investigation[END_REF]. Another useful approach would be to replace the shelter functions by models that simulate airflow through and above the stand directly. However, this would considerably increase the computing time needed.

At the moment, the ForGEM-W model is not able to take into account the acclimation of trees to wind stress, which would be particularly important in edges and stands with sparse canopies. However, due to the simulation setup the simulation of wind prone permanent edges was avoided. On the other hand, trees in different silvicultural systems were growing under different degrees of shelter, and the taller trees in the no thinning mono and unevenaged systems, in particular, would be liable to show acclimation to wind loads. However, that would most likely only decrease the already low levels of damage found for these systems and would not influence the outcomes of this study. Currently, the wind damage module also simulates a storm as a short event in which trees are either damaged or not. In reality, storms can last for hours or even days, causing the wind to penetrate further into the stand as trees start to collapse. However, it is not likely that this would 87 change the order of susceptibilities of the silvicultural systems, since all systems were treated in the same way.

Conclusions

The wind damage module was found to be most sensitive to tree height, DBH and anchorage parameters, and especially to h/d ratio of trees [START_REF] Schelhaas | Introducing tree interactions in wind damage simulation[END_REF]. As a result, it was found in this work that the silvicultural systems leading to low h/d ratios suffered the least amount of damage. Low h/d ratios were attained through low stand densities. Furthermore, low h/d ratios imply a high diameter growth, and thus a shorter rotation, which in turn lowers the time stands are at risk. Especially the uneven-aged system combined a high timber production with low risks.

From the even-aged systems the system free thinning from above mixed had the highest timber production, but included some risk. The no thinning system had the lowest risk, but on account of a considerably lower timber production.

However, low stand densities will also be unfavourable for wood quality, leading to wide growth rings, strongly tapered stems and thick branches. Increasing wood quality through higher stand densities will lead to higher risks of wind damage through increased h/d ratios. Forest owners have to decide individually on this trade-off and how much risk they want to take. This decision will strongly depend on the actual location of a site. The more sheltered a site is, the more focus could be put on aims other than risk consideration. To aid them in this decision, future studies should take into account wood quality and economic parameters, too.

The wind data used in this study are derived from one of the most exposed meteorological stations and risk of wind damage will probably be lower for most of the Douglas fir forests in the Netherlands than as shown here. It is therefore unlikely that foresters will massively adopt a no-thin strategy with low densities, also because of the consequences for wood quality. Currently there is a tendency towards more nature-oriented management. This implies both a move from traditional silvicultural systems like the normal yield table system towards more free thinning and uneven-aged systems, as well as increased interest in admixture of other species. Results from this study indicate that these trends could lead to lower wind risks and even to an increase in overall productivity.

  .................................................................................................................................... 4.1 INTRODUCTION ................................................................................................................................... 4.2 METHOD ............................................................................................................................................. 4.3 RESULTS ............................................................................................................................................. ............................................................................................................................................... 2.1 OUTLINE OF FORGEM-W ................................................................................................................... 2.2 MODEL COMPUTATIONS ...................................................................................................................... 3 RESULTS ................................................................................................................................................... 3.1 STAND DYNAMICS DEVELOPMENT ...................................................................................................... 3.2 TIMBER PRODUCTION AND INTENSITY OF WIND DAMAGE .................................................................... FINDINGS ........................................................................................................... 4.2 MODEL ASPECTS .................................................................................................................................

Figure 1 .

 1 Figure 1. Location of the sample plots.

  [3] Somer's D = (c-d)/t [4] Gamma = (c-d)/(c+d) [5]

Figure 2 .

 2 Figure 2. Predicted values for the probability of damage and amount of damage in function of stand dominant height and stand composition: pure and mixed forest.

Figure 3 .

 3 Figure 3. Predicted values for the probability of damage and amount of damage in function of stand dominant height and stand regeneration mode: high forest and coppice.

Figure 4 .

 4 Figure 4. Predicted values for the probability of damage and amount of damage in function of stand age and stand composition: pure and mixed forest.

Figure 5 .

 5 Figure 5. Predicted values for the probability of damage and amount of damage in function of stand age and stand regeneration mode: high forest and coppice.

Figure 6 .

 6 Figure 6. Predicted values for the probability of damage and amount of damage in function of the mean annual evaporation.

Figure 7 .

 7 Figure 7. Predicted values for the probability of damage and amount of damage in function of the mean temperature of the 3 hottest months .

Figure 8 .

 8 Figure 8. Relationship between the observed and expected numbers of dead occurrence (D=1) and dead absence (D=0) in 10 equal sized groups according to their predicted probabilities.

Figure 1 .

 1 Figure 1. Average number per tree (± standard error) of egg masses, larvae and adult G. scutellatus for 3 levels of damage: 1-<30% of defoliation; 2-30 a 75% of defoliation; 3> 75% of defoliation. Spring 2007. Northern and central regions of Portugal.

Figure 2 .Figure 3 .

 23 Figure 2. Relation between the rate of parasitism (percentage of G. scutellatus egg masses parasitized by A. nitens) and the elevation. Spring 2007. Northern and central regions of Portugal.

Figure 4 .

 4 Figure 4. Dominant height estimated at 10 years old in stands currently under attack (A) and not attacked by G. scutellatus (N), respectively in two time periods: 1995-98 before G. scutellatus arrival (A1 and N1); 2004-06 after G. scutellatus arrival (A2 and N2). For the same stands, bars with different letters mean significant differences (P <0.05). Northern and central regions of Portugal.

Figure 5 .

 5 Figure 5. Timber production estimated at 10 years old in stands currently under attack (A) (n=13) and not attacked by G. scutellatus (N) (n=18), respectively in two time periods: 1995-98 before G. scutellatus arrival (A1 and N1); 2004-06 after G. scutellatus arrival (A2 and N2). For the same stands, bars with different letters mean significant differences (P <0.05). Northern and central regions of Portugal.

Figura 6 .

 6 Figura 6. Relation between the level of damage caused by the G. scutellatus and the percentage of loss in useable timber volume in 13 eucalyptus stands exhibiting moderate to intense defoliation level. Northern and central regions ofPortugal, 2007.

Figures 3

 3 Figures 3.1-4: Stand related parameters affecting the predisposition of forests to infestation by Ips typographus and relative importance (weighting) of the indicators in the hazard model Proportion of Norway spruce

Fig. 1 .

 1 Fig. 1. Relative tree mortality due to I. typographus infestation related to the share of Norway spruce in Węgierska Górka (n=525)

Fig. 10 .

 10 Fig. 10. Relative tree mortality due to I. typographus infestation related to the share of Norway spruce in Ujsoły (n=618)

Figure 1 :

 1 Figure 1: Relative growth loss (radial area) compared to control trees in maritime pine trees of different circumferences, according to classes of Thaumetopoea pityocampa defoliation.

Figure 2 :

 2 Figure 2: Relative growth loss (radial area) compared to control trees in maritime pine trees, according to classes of Thaumetopoea pityocampa defoliation, in the four years following the attack.

Figure 3 .

 3 Figure 3. Transformation of ring growth area from a function of time to a function of number of PPM nests.

Figure 4 :Figure 6 :

 46 Figure 4: Quantile regressions of ring growth area as a function of number of winter nests per ha. 4a. Interior of the Biscarosse plot -4b . At the edge of the Biscarosse plot.

  Number of nests at the edge (EN) = a × Number of nests within the stand (IN) + b Using data collected in 2005 by the INRA Laboratory of Forest Entomology in pure maritime pine stands of the same forest area, we found a significant linear model (N =80 ; P<0.0001 R²= 0.52): EN =1.99 × IN +84.31

Figure 7 :

 7 Figure 7: 90% quantile regression curves for edge and interior tree in the Biscarosse Plot after correction of nests number at the edge.

Figure 7 :

 7 Figure 7: Effect of defoliation classes on the nest number per hectare. The three defoliation classes are [1-5%], [6-25%] and [26-60%].

Figure 8 :

 8 Figure8: Effect of PPM defoliation on relative growth loss using the least square regression on mean values for the interior and the edge of Biscarosse plot. Comparison with results from the empirical approach in "Parentis" and "Leon" stands.

Figure 9 :

 9 Figure 9: Quantile regressions plotting ring growth area as a function of number of winter nests per ha for the "Lit et Mixe" plot

Figure 10 :

 10 Figure 10: Slope of linear quantile regressions as a function of quantiles. Grey areas represent 99% CI. Solid red lines represent least square regression value and dashed red lines the corresponding 99% CI.

Figure 11 :

 11 Figure11: Effect of PPM defoliation on relative growth loss using the least square regression on mean values for the "Lit et Mixe" plot. Comparison with results from the empirical approach in "Parentis" and "Leon" stands.

Figure 12 :Figure 13 :

 1213 Figure 12: Quantile regressions plotting the ring growth area as a function of the mean number of winter nests per tree, for the plot 1 from the "Ventoux". (The plot 1 is shown as an example. The patterns for the plot 2 and 3 are similar.)

Table 4 :

 4 Meta-analyses of the influence of Thaumetopoea pityocampa defoliation rate on tree growth.Defoliation Class(%) Number of studies Mean effect sizeBias corrected

Figure 15 :

 15 Figure 15: Effect of defoliation rate by Thaumetopoea pityocampa on the relative growth loss in pine species.

Figure 16 :

 16 Figure 16: Relative growth loss in individual studies on the response of tree growth to defoliation by Thaumetopoea species. The negative effect sizes for most of the studies indicate that defoliation is almost always responsible for growth loss.Green bars correspond to Thaumetopoea pityocampa, blue bars to Thaumetopoea wilkinsoni and the yellow bar to Thaumetopoea bonjeani.

Figure 1 .

 1 Figure 1. Relative number of observed daily maximum hourly mean wind speeds per wind speed class and wind direction at measurement station Leeuwarden (in % of total number of observations).

Figure 2 .

 2 Figure 2. Development of a) mean (quadratic) diameter, b) stand density, c) tree height of Douglas fir in the even-aged silvicultural systems and d) development of h/d ratio in relation to average height for all systems. All figures show the Douglas fir part only and are averages over 10 replicates.

Figure 4 .

 4 Figure 4. Gross annual increment and its components, averaged over the period needed to get an average diameter of 60 cm or over a period of 100 years. Average of 10 replicates.

Figure 5 .

 5 Figure 5. Wind damage in relation to timber production (volume development plus cumulative harvest), including one standard deviation.

ANALYSIS TO ESTIMATE THE IMPACT OF PPM DEFOLIATION ON TREE GROWTH

  ANALYSIS OF DATA COLLECTED BY LEMOINE IN 1970 ON TWO MARITIME PINE STANDS ............... 2.2 STATISTICAL ANALYSES ..................................................................................................................... 2.3 RESULTS AND DISCUSSION ..................................................................................................................

	2.1 RE-3 DENDROCHRONOLOGICAL, RETROSPECTIVE

2 EMPIRICAL ANALYSIS BASED ON INDIVIDUAL TREE DATA .................................................. . ........................................................................................... 3.1 STAND DESCRIPTION AND CORE COLLECTION ..................................................................................... 3.2 STATISTICAL ANALYSES ..................................................................................................................... 3.3 RESULTS .............................................................................................................................................

Table 1 .

 1 Predictor variables used for building the model.

	Unites	Description

a cont-continuous variable; disc-discrete variable.

Table 2 .

 2 Characterization of the fitting and validation dataset

			s.
	Variables	Fitting data	Validation data
		mean (s.d.) a	mean (s.d.) a
	N (nº trees.ha -1 ) 998 (445.03)	1053 (407.69)
	G(m 2 .ha -1 )	10.53 (6.35)	10.42 (5.53)
	hdom (m)	15.23 (6.01)	16.31 (4.60)
	t (years)	7.17 (3.55)	7.42 (2.71)
	T (˚C)	22.12 (1.23)	22.42 (0.95)
	P (mm)	783.17 (176.16) 746.51 129.13)
	Frost (days)	23.17 (12.52)	23.54 (12.07)
	Evap (mm)	553.35 (67.15) 541.26 (52.29)
	Alt (m)	224.22 (131.01) 220.52 (119.55)

a.

Table 3 .

 3 Estimated logistic regression coefficients (β The order of the variables is the same of the outcome of the stepwise procedure. Only significant variables are shown.

	Model			β i	S.E.	Wald χ 2	Pr > ChiSq	O.R.	C.I.	Type
	Pdead
	Intercept		-9.5117	0.8561	123.44	<.0001
	hdom (2 m)	0.0670	0.0110	37.34	<.0001	1.143	1.095 -1.194 cont
	T (2 °C)			0.4022	0.0283	201.86	<.0001	2.2350	2,001 -2,499 cont
	r			0.6991	0.0947	54.47	<.0001	2.0191	1.964 -2.445 dis
	Comp			-1.1043	0.1025	116.09	<.0001	0.1100	0.074 -0.164 dis
	RI			1.1305	0.1091	107.29	<.0001	9.5930	6.254 -14.714 0/1
	Evap (100 mm)	-0.0036	0.0006	37.09	<.0001	0.6970	0.620 -0.782 cont
	G (2 m	2 . ha	-1 )	0.0420	0.0072	33.54	<.0001	1.2330	1.149 -1.324 cont
	t (2 anos)		0.3791	0.0379	100.02	<.0001	2.1340	1.846 -2,485 cont
	t	2			-0.0175	0.0019	83.05	<.0001	0.9830	0.979 -0.986 cont
	Dead
	Intercept		-1.0850	0.1523	50.7656	<.0001
	N (10)			-0.0012	0.0000 2764.0767	<.0001	0.989	0.988 -0.989 cont
	t	2			-0.0048	0.0003	273.6704	<.0001	3.146	2.763 -3.582 cont
	RI			0.5731	0.0331	299.3290	<.0001	1.559	1.516 -1.604 0/1
	t (2 anos)		0.2220	0.0072	944.3351	<.0001	0.907	0.897 -0.917 cont
	hdom (2 m)	-0.0488	0.0029	283.9868	<.0001	0.389	0.311 -0.482 cont
	Comp			-0.4724	0.0559	71.3358	<.0001	1.185	1.135 -1.238 dis
	G (10 m 2 . ha	-1 )	0.0170	0.0022	58.7124	<.0001	0.956	0.931 -0.981 cont
	r			0.1917	0.0379	25.6319	<.0001	0.514	0.411 -0.637 dis
	Frost (10 days)	-0.0026	0.0005	26.9473	<.0001	0.975	0.965 -0.984 dis
	T (2 °C)			-0.0145	0.0064	5.1764	0.0229	0.971	0.948 -0.996 cont

i ) and respective standard errors (S. E.), Wald test, p-value, odds ratio (O. R.) and odds ratio confidence intervals (C. I.) for the logistic regression models of probability of dead (Dead) and proportion of dead (Pdead) from the whole data a . a

Table 4 .

 4 Non parametric measures of association of predicted probabilities and observed responses for the probability of dead (Dead) and proportion of dead (Pdead) logistic models.

	Models	Concordant Discordant Somers'D Gamma	Tau-a	c
	Pdead	79.9%	19.9%	0.601	0.602	0.275	0.800
	Dead	64.0%	35.2%	0.288	0.290	0.072	0.644

Table 4 .

 4 Mean bias and precision, and its standard deviation, for the predicted values of the probability of tree mortality. Predicted values are estimated by combining sequentially the probability of dead occurrence (Dead) and of the proportion of dead trees (Pdead).

		N	Mean	Std deviation
	Bias	7510	0.0209143	0.1369942
	1-Precision	7510	0.0875716	0.1074014

Table 5 .

 5 Contingency table and indicators for the accuracy of the proportion of dead trees model (Pdead) using the validation dataset.

	Observed data
	damage occurrence damage absence	Total

Wood loss impact due to the eucalyptus snout beetle, Gonipterus scutellatus, in eucalyptus stands, relation to defoliation and population Authors: Branco, M.

	Tomé. SIMPLOT: Simulating the impacts of fire severity on sustainability
	of eucalyptus forests in Portugal. Ecological Indicators In Press, Corrected Proof.
	Bozic, M., O. Antonic, R. Pernar, S.D. Jelaska, J. Krizan, J. Cavlovic, and V. Kusan. 2004.
	Modelling the damage status of silver fir trees (Abies alba Mill.) on the basis of
	geomorphological, climatic and stand factors. P. 202-208 in 4th International Workshop on
	Environmental Applications of Machine Learning (EAML), Bled, SLOVENIA.
	Branco, M., C. Boavida, N. Durand, J.C. Franco, and Z. Mendel. 2009. Presence of the
	Eucalyptus gall wasp, Ophelimus maskelli (Ashmead) and its parasitoid Closterocerus
	chamaeleon Girault, in Portugal: first record, geographic distribution and host preference. P.
	51-54 in Phytoparasitica.
	Brasier, C.M., and J.K. Scott. 2008. European oak declines and global warming: a theoretical
	assessment with special reference to the activity of Phytophthora cinnamomi. P. 221 -232 in
	EPPO Bulletin.
	Ferreira, C., J. Carvalho, and M. Baptista. 2005. Instruções para o trabalho de campo do
	Inventário Florestal Nacional -IFN 2005/2006. Florestal, D.N.d.G. (ed.). Autoridade Florestal
	Nacional, Lisboa.
	Ferreira, M., A. Manta, and C. Valente. 2006. Primeiro registo de um ácaro eriofídeo do
	eucalipto em Portugal, Rhombacus eucalypti Ghosh & Chakrabarti. P. 227-229 in Agronomia
	Lusitana
	Florestais, D.d.I.e.E. 1999. Manual de Instruções para o Trabalho de Campo do Inventário
	Florestal Nacional. Florestas, D.G.d. (ed.), Lisboa.
	Gadgil, P.D., and J. Bain. 1999. Vulnerability of planted forests to biotic and abiotic
	disturbances. P. 227-238 in New Forests.
	Gonzalez, J.R., M. Palahi, T. Pukkala, and A. Trasobares. 2006. Modeling the risk of forest fires
	in Catalonia (North-East Spain). P. 85-91 in Seminar on Scientific Tools and Research Needs

Table 1 :

 1 Forest Management Alternatives (FMAs):

		Nature reserve	Close-to-Nature	Combined-objectives Intensive even-aged	Wood-Biomass
	Stand Composition	mixed spruce forest, uneven-aged	mixed spruce forest, uneven-aged	pure -uneven-aged	pure, even-aged	pure, even-aged
	Proportion of Norway Spruce	variable	variable	100%	100%	100%
	Regeneration Type	natural	natural	planting and natural	planting (and natural)	planting
	Thinning Regime	no	selective	several moderate thinning operations in the course of the rotation period	several moderate period thinning operations in the course of the rotation	1-2 heavy thinnings
	Harvesting Regime	no	selective	strip and femel system	clear cut at 80	clear cut
	Rotation Period	no	100 years	100 years	80 years	40 years
	Canopy Closure	variable	variable	variable	homogenous 80%	homogenous 60%

Area of high to very high predisposition to Ips typographus

  Areas of high to very high susceptibility to infestation by Ips typographus increased with stand age. The most abrupt change from less-disposed young stands to high-disposed older stands was observed for the variant "intensive even aged" (INT). The assumption of rotation periods shorter than 100 years (INT and BIOM) anticipated over-maturation of forest stands (comp. Fig.4.1-3) and thus levels of very high stand-related predisposition to Ips typographus. However, given the alternative BIOM, predisposition levels were constantly high due to low values of canopy closure and the dominance of Norway spruce. Proportion of predisposition classes (low, medium, high, very high predisposition to Ips typographus) for the different management alternatives in time steps of 20 years from presence to 100 years in future

					BIOM
		50.000
		45.000
		40.000
		35.000
	area (ha)	20.000 25.000 30.000					<= 20 yrs 21-40 yrs
		15.000
		10.000
		5.000
		0
		presence	+20 yrs	+40 yrs	+60 yrs	+80 yrs	+100 yrs
	Close-to-nature presence +20 yrs +40 yrs +60 yrs Combined objections +80 yrs +100 yrs presence +20 yrs +40 yrs +60 yrs +80 yrs +100 yrs +20 yrs +40 yrs +60 yrs Figures 4.4-7: Intensive even aged 0% 10% 20% 30% 40% 50% 60% 70% 80% 90% 100% 0% 10% 20% 30% 40% 50% 60% 70% 80% 90% 100% 0 5.000 10.000 15.000 35.000 40.000 45.000 50.000 presence 25.000 30.000 area (ha) 20.000	very high high medium low very high high medium low +80 yrs	+100 yrs	CN CO INT BIOM
		100%
		90%
		80%
		70%
		60%
							very high
		50%					high
		40%					medium
		30%					low
		20%
		10%
		0%
			presence +20 yrs +40 yrs	+60 yrs	+80 yrs +100 yrs
						end of rotation period

Figures 4.1-3: Proportion of age classes assuming different rotation periods (CN, CO: 100yrs; INT: 80yrs; BIOM: 40yrs) in time steps of 20 yrs from presence to 100 years in future Figure 4.8: Forested area of high to very high predisposition to I. typographus -comparison of trends between the management alternatives close-to-nature (CN), combined objectives (CO), intensive even-aged (INT) and biomass production (BIOM)

1 Table 1 .

 11 Relative tree mortality due to I. typographus infestation in Węgierska Górka -effect of selected factors estimated using Kruskal-Wallis test

	Factor	sample size (n)	df	H	p
	number of species (1-6)	527	5	17,73081	0,0033
	share of Norway spruce (10-100%)	520	9	29,42919	0,0006
	share of broadleaved (0-100%)	526	9	13,47323	0,1424
	exposure (flat excluded)	525	7	12,08510	0,0978
	altitude (4 zones)	527	3	13,27106	0,0041
	age variability	523	3	1,24409	0,7425
	age (20-years classes)	523	9	54,42876	0,0000
	BHD (10-cm classes)	523	6	36,97471	0,0000

Table 2 .

 2 Relative tree mortality due to I. typographus infestation in Ujsoły -effect of selected factors estimated using Kruskal-Wallis test

	Factor	sample size (n)	df	H	p
	number of species (1-4)	619	3	10,02266	0,0184
	share of Norway spruce (10-100%)	618	9	28,10486	0,0009
	share of broadleaved (0-100%)	617	9	20,09340	0,0174
	exposure (flat excluded)	595	7	12,82806	0,0765
	altitude (4 zones)	618	3	13,35762	0,0039
	age variability	616	3	4,40583	0,2209
	age (20-years classes)	616	6	12,32952	0,0550
	BHD (10-cm classes)	616	5	3,96178	0,5549

Table 1 : Results of ANOVA: impact of both circumference and PPM defoliation classes on maritime pine radial tree growth one year following the main attack. Parentis stand

 1

	Source	Degrees of freedom	F	Pr > F
	circumference class in 1969	2	138,7	< 0,0001
	defoliation class in 1970	3	20,1	< 0,0001
	circumference class × defoliation class	6	0,6	0,693
	Leon stand
	Source	Degrees of freedom	F	Pr > F
	circumference class in 1969	3	14,9	< 0,0001
	defoliation class in 1970	2	10,8	< 0,0001
	circumference class × defoliation class	6	0,6	0,702

Table 2 : Results of ANOVA: impact of both circumference and PPM defoliation classes on maritime pine tree growth in the four years following the main attack. Parentis stand

 2

	Year			N+1		N+2		N+3		N+4
	Source	Df	F	Pr > F	F	Pr > F	F	Pr > F	F	Pr > F
	defoliation	3	14,3 < 0,0001	1,6	0,197	1,4	0,263	2,4	0,073
	circumference	2	111,6 < 0,0001 125,9 < 0,0001 145,1 < 0,0001 143,0 < 0,0001
	defoliation × circumference	6	0,7	0,665	0,9	0,455	1,2	0,287	0,7	0,663
	Leon stand
	Year			N+1		N+2		N+3		N+4
	Source	Df	F	Pr > F	F	Pr > F	F	Pr > F	F	Pr > F
	defoliation	3	11,7	< 0,0001	11,3	< 0,0001	3,5	0,018	1,3	0,288
	circumference	2	8,2	0,000	2,2	0,116	2,3	0,107	0,5	0,600
	defoliation*circumference	6	0,4	0,860	2,4	0,032	1,6	0,139	1,6	0,142

Table 3 : Results from ANOVA: effect of the mean percentage of defoliation on the total number of nests per ha.

 3

	Source	Degrees of freedom	F	Pr > F
	Class of % defoliation	2	231,04	< 0,0001

Table 5

 5 Mean Relative Growth Loss (%) for two PPM defoliation classes according to three analytical methods.

	Method			Dendrochronological			Empirical	Meta-analysis
	Defoliation classes	Biscarosse interior	Biscarosse edge	Ventoux (Plot 1)	Ventoux (Plot 2)	Ventoux (Plot 3)	Lit et Mixe	Parentis	Leon	17 cases
	6-25%	-25,2	-18,9	-6,0	-9,6	-10,0	-11,8	-14,0	-18,7	0
	26-60%	-57,0	-42,8	-14,6	-23,5	-24,4	-26,6	-26,6	-43,4	-33

Table 1 .

 1 Relevant parameters used in ForGEM-W.

	Parameter	Douglas fir Beech	Unit	Reference
	Height of crown centre,				Dutch Forest Reserves
	relative to crown length	0.33	0.41	-	Database
	CR MAX (maximum crown				Dutch Forest Reserves
	radius)	6.5	14	m	Database
	RUE	2.56E-09	1.78E-09	kg DM MJ -1	calibrated
	T (onset of age effect on
	RUE)	30	30	year	calibrated
	c (impact of age effect on
	RUE)	0.2	0	-	calibrated
	b 2 (speed of crown				Estimated from Jansen
	expansion)	0.1137	0.0722	-	et al., 1996

, as an output, a tool, which allows for analysis of sustainability impacts of existing and future FWCs.

The European Forest Institute (EFI) kindly offered the EFORWOOD project consortium to publish relevant deliverables from the project in EFI Technical Reports.

The reports published here are project deliverables/results produced over time during

88

To conclude, one should still be cautious to transfer the results of this study directly to other wind climates, tree species and soils. More or less severe wind climates might change the specific order of silvicultural systems as obtained in this study because underlying relationships are non-linear [START_REF] Schelhaas | Introducing tree interactions in wind damage simulation[END_REF]. Other tree species will differ in growth patterns and react differently to inter-and intra-tree species competition, too. However, conclusions on the effects of stand density and h/d ratios on wind stability are likely to be generally valid.