Modelling damage occurrence by snow and wind in forest ecosystems

Snow and wind damages are one of the major abiotic disturbances playing a major role in forest ecosystems and affecting both stand dynamics and forest management decisions. This study analyses the occurrence of wind and snow damage on Norwegian forests, based on data from four consecutive forest inventories (1995–2014). The methodological approach is based on boosted regression trees, a machine learning method aiming to demonstrate the effects of different variables on damage probability and their interactions as well as to spatialize damage occurrence to make predictions. In total, 313 models are fitted to detect trends, interactions and effects among the variables. The main variables associated with damage occurrence are consistent across all the models and include: latitude, altitude and slope (related to site and location), and tree density, mean diameter and height (related to forest characteristics). The results show that stand dominant height is a key variable in explaining damage probability, whereas stand slenderness has a limited effect. More heterogeneous forest structures make birch dominated stands more resistant to damage. Finally, the models are translated into occurrence maps, to provide landscape-level information on snow and wind damage hazard. Further application of the models can be oriented towards assessing the probability of damage for alternate stand management scenarios.


Introduction
Snow and wind damage in forests can result in large economic losses and changes in stand dynamics, thus their consideration is vital if we want to improve the chances of obtaining the specific management objectives (Gadow Von, 2000). Such a consideration is even more relevant under changing conditions, which can happen due to climatic variations or shifts in the forest state. It is not precisely clear how climate is going to evolve in Northern Europe (Feser et al., 2015;Mölter et al., 2016), but climatic changes leading to higher frequency and intensity of the disturbances agents may increase snow and wind damage in forests. Changes in the state of forests can occur, for example, by an increase in the forest age and/or variations in the structure. If forests become more vulnerable, thereby increasing their probability of experiencing snow and wind damage, we are to expect also modifications in biodiversity (Lain et al., 2008;Thom and Seidl, 2015), carbons stocks (Schelhaas et al., 2003;Thom and Seidl, 2015) or the provision of goods and services, such as timber (Zubizarreta-Gerendiain et al., 2017). Thus, it is important to investigate why, and which, forests are vulnerable to snow and wind damage, if we want to adopt management actions that reduce the uncertainty in the expected outcomes.
Several studies, using different approaches and tools, have assessed the forest responses to snow and wind damages through modelling (e.g. Seidl et al., 2011). Some of these studies have focused on the effects following high-intensity storms that lead to catastrophic damages events, and only include the affected areas (Dobbertin, 2002;Valinger and Fridman, 2011;Hale et al., 2015). Approaches that only focus on catastrophic events and affected areas obscure the reasons behind why certain stands are more vulnerable than others, because such events have devastating effects on all stands, regardless of the variations in the forest stability or the resistance traits of certain trees. In this regard, instead of targeted post-event inventories, we required large, reliable and comprehensive datasets covering larger temporal and spatial scales. Long-term observations provide a better understanding of the overall forest susceptibility to damage, and relate such damage to forests and site variables (Martín-Alcón et al., 2010;Albrecht et al., 2010;Valinger and Fridman, 2011), while large spatial scales capture the variability of disturbance regimes and forest characteristics.
Models that aim to predict snow and wind damages in forests can be hybrid-mechanistic or empirical. Hybrid-mechanistic models com bine knowledge and processes information, facilitating the analysis of potential different management scenarios. Examples of hybrid-mechanistic models have been listed by Gardiner et al. (2008). Mechanistic models have the advantage of being able to include tree failure in a physical manner, for example by focusing on the critical wind speeds needed to uproot and break the stem of trees (Gardiner et al., 2000;Schelhaas et al., 2007). On the other hand, and although with exceptions (Nicoll et al., 2006), mechanistic models were developed based on a limited number of observations, and only dealing with uniform stands of single species, leaving opportunities to improve the representation of wind and snow damage processes on more heterogeneous stands. Empirical models can account for all of the ecological complexity reflected in the collected dataset, and help to identify the key factors involved. Studies using empirical approaches often diverge in the relative importance of the variables, due to the obvious variations in climate, location and forest conditions among studies from different areas (Mitchell, 2013). Predictors used in these empirical studies have included: variables concerning the topographic situation (Mitchell et al., 2001;Albrecht et al., 2012;Pasztor et al., 2014;Hanewinkel et al., 2014); variables referring to the soil and site conditions (Mitchell et al., 2001); and variables describing stand-and tree-level characteristics such as mean basal area (Albrecht et al., 2010;Valinger and Fridman, 2011;Pasztor et al., 2014), density (Pukkala et al., 2016), tree size (height, diameter and slenderness) (Martín-Alcón et al., 2010;Albrecht et al., 2010;Pukkala et al., 2016), stand structure (Hanewinkel et al., 2014;Pukkala et al., 2016), dominant species and species diversity (Albrecht et al., 2010;Valinger and Fridman, 2011;Albrecht et al., 2012;Pasztor et al., 2014).
One of the disadvantages of the traditional empirical modelling, such as logistic models (Fridman and Valinger, 1998;Jalkanen and Mattila, 2000) is the difficulty in understanding the complex damage process from the resulting model. Identifying how variables are involved in, and affect, forest vulnerability to snow and wind is complex. There are many correlations among the variables involved, making it difficult to interpret the influence of certain of these without considering the correlated variables. New methods based, on machine learning, such as boosted regression trees, can consider a large range of variables and, at the same time, fit complex non-linear relationships, accommodate missing values, handle variable interactions and still deliver good predictive models (Elith et al., 2008). Although machine-learning models are difficult to interpret, they can be summarised in a way that provides powerful insights, and represent a valuable tool for use in making informed decisions. Several studies have proved that machine-learning approaches are powerful tools for identifying characteristics related to damage phenomena, and providing a clear understanding of the meaningful variables Albrecht et al., 2012).
In this study we focused on the probability of the occurrence of snow and wind damage in forest stands. The analysis relies on empirical data, collected from several consecutive National Forest Inventories from Norway, including 20 years' worth of data. Based on an adequate methodological framework that addresses the complexity of the data structure, the specific objectives of this contribution are to (i) assess which variables are relevant to explain the damage occurrence, evaluating their relationships and interactions (ii) identify areas susceptible to damage spatializing the probability of occurrence.

Data sources
The data used was from the Norwegian National Forest Inventory (NFI), collected during the period 1995-2014, corresponding to the 7 th , 8 th , 9 th and 10 th inventories. The Norwegian NFI is a systematic inventory of permanent circular plots of 250 m 2 in a 3 X 3 km grid represent ing the entire distribution and diversity of Norwegian forest types. For plots divided by a stand border (e.g., a plot is part forest, part water, or encompasses two forest types with very different productivity parameters), only the most forested part of the plot was considered. Plots were also divided into four different forest types, based on their dominant species and according to their basal area presence: spruce (spruce > 70%); pine (pine > 70%); birch (birch > 70%); and mixed forests (when no single species accounted for more than 70% of the basal area).
Most of the observations were categorized as mixed forests and birch forest types, 62 and 10% respectively; conifer dominated plots represented the 12 and 16% for pine and spruce, respectively. The occurrence of forest damage was identified and measured from an area of 1000 m 2 around the centre of each permanent plot; we considered those plots where such damage was recorded in that area to be damaged plots. Following the inventory protocols (Skoglandskap, 2007), stand damage occurrence was recorded when: it was believed to have taken place during the five years prior to measurement; it had an effect on the future economic development of the stand; and it compromised regeneration, or decreased volume production and wood quality. We considered snow and wind damage together, as climatic conditions in Nordic areas make the combination of both agents the cause of damage (Valinger and Fridman, 1999).
The 7th NFI (1995)(1996)(1997)(1998)(1999) was used as a reference inventory to provide information concerning initial stand characteristics, and the estimated variables that described the forest stand and site were used as predictors ( Table 1). The predicted variable was presence/absence of snow and wind damage, and this was defined and evaluated for those plots included in the 8th-10th inventories (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014); damage presence was assigned to those plots that were recorded, at least once, as Gini coefficient [0,1] where N is the number of trees in the plot, i, j two different trees in the plot, g i the basal area of the tree i, g j the basal area of the tree j and the mean basal area of the plot. Shannon coefficient Species diversity index, calculated by the basal area, for Spruce, Birch, Pine, other conifer and other broadleaves, before the damage. Higher values indicate more diverse plots, lower values indicate less diverse plots. Shannon index = [0, ln S] where p i is the basal area probability: , being g i the basal area of the i specie and G the total basal area of the plot. i, is the number of species ranging from 1…S Ecological Modelling xxx (xxxx) xxx-xxx damaged, whereas damage absence was assigned to plots with no damage recorded in the three inventories (15 years). In total, 7097 plots were considered for modelling purposes, of which 6486 plots were not classed as damaged.

Statistical methods
The statistical analysis was based on boosted regression trees (BRT). This approach combines machine learning and statistical techniques, and consists of adding regression trees fitted in a forward stepwise process (Friedman et al., 2000;Elith et al., 2008). The stepwise process starts with the creation of a single tree that minimises the loss function; then a new tree is fitted, using the residuals of the first one. By these means, the first developed model is updated to incorporate the information from the new tree. The residuals of the latest updated model are recursively used to fit a new tree, and so forth. The existing trees do not change, as only the fitted values of each observation are re-estimated each time a new tree is added. The final model reflects the contributions of all of the developed trees, providing more stable and accurate values than would a single tree.
The BRT model calibration is defined by four parametersthe learning rate (or shrinkage parameter), the bag fraction, the tree complexity, and the number of trees required for optimal prediction. The learning rate determines the contribution of each new tree to the growing model, and it is always substantially lower than 1, higher values being related to faster learning. The bag fraction provides information on which fraction of the entire data should be drawn randomly to fit the new tree. This parameter includes a random probabilistic component, making each run model different, and is aimed at improving model accuracy, speed of model creation, and the reduction of overfitting (Friedman, 2002). Tree complexity controls the number of fitted interactions among variables, and determines the number of splits in each tree; for example, a value of 1 will present only one split, meaning that the model does not consider interactions; a value of 2 will result in two splits, and two interactions. The optimal number of trees is selected based on the three previous parameters. The values fitted by the final model are computed as the sum of all of the predictions of the trees, multiplied by their respective learning rates.
We built different models, based on alternative values for learning rate (0.005, 0.001, 0.01, 0.05, 0.1), bag fraction (0.2, 0.5, 0.8) and tree complexity (every unit between 1 and 10), and 10-fold cross-validation was performed to estimate the optimal number of trees in each case, following Elith et al. (2008). The models were fitted assuming a Bernoulli distribution, since the response variable was the presence/absence of snow and wind damage, and were designed for each of the forest types and for all of the forest types combined, resulting in a total of 750 models (5 forest types + combined, 5 levels of learning rate, 3 bag fractions and 10 levels of tree complexity). The statistical analyses were developed in R version 3.3.2 (R Core Team, 2013), and the BRT models were based on the 'dismo' extension of the 'gbm' package (Elith, 2017).
The deviance value of the models was used to select the model with the best predictive performance for each forest type. Additionally, the performance and behaviour of all developed models was analysed. For each model, the cross-validation deviance was calculated, and analysed as a function of the different model parameters (tree complexity, bag fraction and learning rate). The relative importance of the predictor variables was calculated by considering the number of times each variable was used for splitting, weighted by improvement to the model as a result of the split, and averaged across all of the trees in the models, following Elith et al. (2008). Model behaviour was evaluated through partial dependence plots, showing the effect of each of the variables on the response by accounting for the average effects of all other predictors in the model. In order to identify the fitted interaction effects, we tested the partial dependencies between the response and predictor variables. The model predictions were obtained for each pair of predictor variables, setting all other predictors to their means. We plotted the most significant pairwise interactionsthose variables with higher relative significance. Finally, the models with the best predictive performance were used to create a predictive map of the estimated probability of damage occurrence. The map was created by using observations from the period 2010-14, divided by forest type, as the data input. The predicted values were presented in the map as a spatially-weighted average, calculated through Gaussian kernel smoothing.

Results
Some models could not be computed for the parameters considered due to the lack of sufficient data. In total, 92 models were developed for all of the forest types combined, 63 models were developed for spruce-dominated stands, 68 for birch stands and 90 for mixed stands; due to the reduced number of damaged plots, pine only produced one model, which was not included in the study. We found that, as we increased the BRT parameters, we obtained better performance levels, thereby increasing the randomness in the models, the learning rate and the tree complexity (Fig. 1). Variation in the BRT parameter values also changed the number of trees that were fitted per model; increasing learning rate and tree complexity decreased the number of trees. Variation of the BRT parameter values also provided different model performance values. Relatively few models were fitted with the fastest feasible learning rates, and none with the highest learning rate considered (0.1). In all of the developed models, reducing tree complexity decreased model deviance, although such variation was limited. In models predicting damage occurrence for all species together, increasing bag fraction reduced model deviance. For all of the models, increasing learning rate and tree complexity reduced the number of optimal trees required for minimum error.
Variable importance of the best predictive models (Fig. 2), and across all models, were consistent. Among the site-related variables, latitude, altitude and slope were the most significant. Density, diameter and height were the three most significant variables related to forest characteristics, although the significance of diameter varied for each forest type. On the other hand, the dominating species had a low contribution to the models predicting damage occurrence for all species together.
The marginal effect of each of the variables was consistent across all models (Figs. 3 and 4). Steeper slopes made all forest types more susceptible to being affected by wind and snow. Similarly, northernmost forests at higher elevations had a higher probability of damage occurrence, particularly on mixed forests, or when the data was analysed for all forest types without previous segregation. How forest characteristics affected the probability of damage occurrence varied depending on the dominant species. Increasing stand dominant height above 20 m caused a sharp increase in the probability of damage occurrence on spruce and mixed stands. Increasing mean diameter in the initial stages of development for stands with small mean diameters reduced the probability of damage on mixed stands, while relatively small increments, for diameters below 40 cm, increased the probability of damage occurrence in birch-and spruce-dominated stands. Increasing density up to a 2000 tree ha −1 threshold increased the damage occurrence probability in all of the forest types. Increasing slenderness did not influence the damage occurrence probability on any of the forest types. Increasing stand heterogeneity, once a threshold level of 0.5 was reached on the Gini coefficient, had a negative impact on the damage occurrence probability, except for in birch stands, where highly heterogeneous stands appeared to be less susceptible to being affected by snow and wind. Finally, species diversity, as explained by the Shannon index, did not have an effect, except for high values in mixed Ecological Modelling xxx (xxxx) xxx-xxx Fig. 1. Models deviance as a function of the tree complexity (number of interactions), bag fraction (stochasticity) and learning rate (shrinkage). The coloured dots correspond with the selected model in each forest type with minimum deviance resulting from the best combination of bag, tree complexity and learning rate values. stands, where highly diverse stands were more prone to suffer damage episodes.
The nature and magnitude of the fitted interactions varied across the selected models (Fig. 5). The model for all species together had the highest interaction levels between altitude and latitude, stand dominant height and latitude, and density and altitude. In the selected model for spruce-dominated stands, the most important pairwise interactions of continuous variables were between density and stand dominant height, Gini coefficient and altitude, but their interaction values were much lower compared to the model for all species together. In the selected model for birch-dominated stands, the strongest interaction was between the Gini coefficient and stand mean diameter. In the selected model for mixed stands, the most important pairwise variables, due to their strongest interactions, were between latitude and altitude, Shannon index and latitude, slenderness and altitude, and dominant height and diameter. In the selected model for mixed stands, the interaction between latitude and altitude had the same pattern as the one presented for all species together, and therefore was not plotted again in Fig. 5.
The most important interactions of the models presented different effects. In the model predicting damage occurrence for all species to gether, the highest fitted values occurred when the stands were at high altitude regardless of latitude values, when the stands were at high altitude regardless of stand density, and when the stands had dominant heights above 20 m for all latitude values. In the model predicting damage occurrence for spruce-dominated stands, the highest fitted values occurred when the dominant height was above 25 m for all density values. On birch-dominated stands, the highest damage probability values occurred when the stand was older and even-agedhigh stand mean diameter and homogeneous structure (low Gini coefficient values). In the model predicting damage occurrence for mixed stands, the highest fitted values occurred when the stand was less diverse (lower Shannon index values) at high latitudes, at high altitude regardless of the slenderness, and when the mean stand diameter was low, and the dominant height was high.
Finally, using the best predictive models for each forest type, a spatial prediction was developed in order to identify vulnerable areas (Fig.  6). According to the model predictions using the most up-to-date plot data available, the most vulnerable areas for spruce-and birch-dominated were in the south-west of the country, although birch-dominated vulnerable areas extended across the whole country. Mixed stands were most vulnerable in the north of the country. Díaz-Yáñez et al. Ecological Modelling xxx (xxxx)

Discussion
Among the different disturbances resulting in forest damage, snow and wind play an significant role in Norway (Díaz-Yáñez et al., 2016). Analysis of damage occurrence, however, presents serious challenges. On one hand it requires a long period of time to effectively reveal patterns and trends, and on the other hand, the potential exploratory variables are complex, often interact with each other, and do not necessarily present a linear relationship with the probability of damage occurrence. The present study addressed this first challenge by covering an interval of 20 years, which resulted in accessible information on dam age occurrence of 15 years from the measurements, making it possible to compare the characteristics of damage and undamaged stands, and to use that information to identify which variables and values make a stand more vulnerable to experiencing damage. The second challenge, related to the complexity of potential exploratory variables, was addressed by using a BRT-based approach. BRT is a machine-learning method that provides an accurate description of variable relationships and interactions, having the advantage of being flexible, providing useful insights into complex ecological datasets, making it possible to identify complex interactions among variables, and facilitating interpretation of model behaviour with a high level of accuracy.
Unfortunately, one of the limitations of this approach relates to the applicability of the models to other areas or conditions; whereas predictions resulting from BRT models have a good level of accuracy compared to other modelling approaches , the original models are not easily transferable, as is the case for regression models, for example, where the estimates of the parameters would suffice for replication in other conditions or locations. Another limitation concerns the risk of overfitting. Although BRT models may allow overfitting, this is appropriately addressed by using different mechanisms (e.g., introducing stochasticity into the fitting process); also, it has been proved that overfitting does not compromise the predictions, even providing better performance than other methods, such as logistic models (Elith et al., , 2008Leathwick et al., 2006). Arguably, the higher predictive power, as well as the capacity to reveal complex interactions, compensates for the overfitting limitation, as reflected in its increasing popularity. Particularly in forest sciences, the use of BRT has increased in recent years. It has been used to investigate storm resistance between species (Albrecht et al., 2012), to map site indices of different species (Aertsen et al., 2010), to estimate yield biomass (Van Meerbeek et al., 2014;Mola-Yudego et al., 2016) and forest productivity (Verkerk et al., 2015), to account for the geographical variation in soil organic carbon stocks (Nanko et al., 2017), and to estimate species or disease distribution and richness, related to ecological variables (Elith et al., , 2008Leathwick et al., 2006;Randall and van Woesik, 2015), among others. In this study, we assessed the effects of different values for the different model parameters, and analysed their effects on the predictive performance of the models. As discussed in previous studies (Leathwick et al., 2006), this is a necessary step for Fig. 3. Marginal dependency plots for the three site variables with highest variable importance and their comparative effect among the three forest types (Spruce, Birch and Mixed) and all of them together (All). The solid lines represent the selected models and the coloured areas the variability around the selected models from all the developed models. Y axes are centred to have zero mean over the data distribution. Díaz-Yáñez et al. Ecological Modelling xxx (xxxx) xxx-xxx Fig. 4. Marginal dependency plots of the six forest variables with highest variable importance across all the models and their comparative effect among the three forest types (Spruce, Birch and Mixed) and all the forest types together (All). The solid lines represent the selected models and the coloured areas the variability around the selected models from all the developed models. Y axes are centred to have zero mean over the data distribution.
each dataset, as preselecting the values of the model parameters could lead to models with higher predictive errors. Although the response variable was measured on a bigger area than the plot and, in consequence, may relate to other trees than those described in the plot, the plot-level information is a valid representation of the stand, and therefore highly correlated to the damage measurement. The portfolio of variables considered in predicting the damage occurrence was in agreement with previous studies (Albrecht et al., 2012;Mitchell, 2013;Pasztor et al., 2014;Hanewinkel et al., 2014;Pukkala et al., 2016), and the results showed that their role across all the developed models was consistent. Concerning variables related to location, latitude and altitude had a contribution in our models, as the spatial distribution of the stand plays an important role in how exposed it is to damaging agents; typically, stands located at higher altitudes and more northerly latitudes are more exposed to snow and wind. Considering the levels of the damaging agents in the models would have helped to better understand the effect of the exposure to damaging agents, however climate information was not available for the NFI plots and times of measurement. Previous studies have also considered topographic information, by using a variable proxy for topographic exposure (Scott and Mitchell, 2005), whilst in this study we used directly measured variables instead, as we want to understand their behaviour individually. Our results showed, however, that the combined effects of altitude and latitude did not necessarily increase the probability of damage occurrence, and did not equally affect all species. This can be explained by the growing conditions of the trees located in northern latitudes that make them more adapted to extreme conditions (Nicoll et al., 2008), and even more resistant than stands in more sheltered areas. Although other studies have found negative correlations between altitude and damage (Lanquaye-Opoku and Mitchell, 2005;Albrecht et al., 2012), they have also concluded that there are many other factors that might be relevant when considering topographic information; for example, local conditions such as stand borders, which are very relevant in relation to wind damage, and the effects of local climatic conditions, as in this study, with many coastal areas and high elevations. From our results, there was also a clear positive correlation between slope and damage occurrence, which contradicts previous wind studies (Hanewinkel et al., 2004;Lanquaye-Opoku and Mitchell, 2005). Díaz-Yáñez et al. Ecological Modelling xxx (xxxx) xxx-xxx Concerning variables related to the forest conditions, the results underlined the importance of tree's height, this being consistently the most significant variable in all of the models, which is acknowledged as being a primary variable for estimating stand vulnerability to wind damage (Schmidt et al., 2010). Diameter was the second most significant variable, also indicating that increasing diameter increases vulnerability, which contradicts many previous studies, but agrees with Canham et al. (2001) and Pukkala et al. (2016), probably indicating that the average stand height was also increasing due to the high correlation between diameter and height (Pukkala et al., 2016).
In addition to height and diameter, independently, previous studies have also discussed the importance of their ratio in relation to damage vulnerability (Fridman and Valinger, 1998;Peltola et al., 1999a). Our results agree with the general idea that increasing trees slenderness, or height, results in higher stand vulnerability; however, the results showed that height explained a more substantial part of the probability of damage occurrence, whereas slenderness had a very limited effect. The fact that slenderness provides information on tree-level stability and not stand level could explain its lack of influence in our models. There are recommendations that suggest that only looking at slenderness should be avoided, as it could lead to a false conclusion (Albrecht et al., 2010); for example, lower ratios have been found to represent more stability against snow breakage damage (Peltola et al., 1999b), but lower values could also represent trees with larger crowns that easily store snow, thus making them more vulnerable (Weiskittel et al., 2009). Our study's conclusions are in the same direction, indicating that mean stand slenderness by itself is not a good predictor for determining stand-level resistance to snow and wind damage occurrence, although it has proven to be a relevant variable in determining damage intensity and type in spruce and mixed stands (Díaz-Yáñez et al., 2017).
Although the individual form of the trees in a stand is highly relevant in determining their vulnerability, it can be greatly influenced by stand density. Previous studies have highlighted that denser stands, with taller trees, are less stable, while low-density stands result in the development of more resistant trees (Petty and Worrell, 1981;Mitchell O. Díaz-Yáñez et al. Ecological Modelling xxx (xxxx) xxx-xxx On the other hand, other studies have not found a relationship between density and a decrease in stand vulnerability (Hanewinkel et al., 2014). Our results showed that dominant height and mean diameter were highly correlated with stand density in spruce-dominated stands, indicating that denser stands, with higher dominant heights and larger mean diameters, had higher probabilities of damage occurrence. It is assumed that stand structure has an influence on stand stability, making more heterogeneous structures more resistant to damage from snow and wind (Griess and Knoke, 2011;Pukkala et al., 2016). In this study, the effects of stand structure were analysed through the Gini index (Valbuena et al., 2012), which indicated that more heterogeneous stand structures reduced the probability of damage occurrence in birch-dominated stands, but not in spruce or mixed stands, or when considering all species together. It has been also discussed that irregular structures tend to have individual trees with lower slenderness, making the trees less vulnerable (Kenk and Guehne, 2001;Mason, 2002). Although silvicultural treatments oriented towards more irregular stands can result in more stable stands, it would be problematic to change from an even-to uneven-aged structure in a stand, as that would leave individual dominant trees that were not adapted to increased wind loading exposed (Mason, 2002;Martín-Alcón et al., 2010). It is also relevant to consider the site-specific biomechanical adaptation of trees to agent exposure, as acclimatising trees to a certain level of exposure will fail if the exposure levels are exceeded (Coutts and Grace, 1995;Bonnesoeur et al., 2016). Therefore, silvicultural treatments aimed at modifying stand structure, such as thinning (Pukkala et al., 2016), should be applied with caution, and other related parameters, such as tree anchoring through soil depth and the thigmomorphogenetical acclimation of trees, should be considered. Moreover, silvicultural treatments that change forest structures should aim for slow transformation, making changes in the early stages, when stand heights are low.
In the models for all species combined, the variable defining the species had a mean contribution below 5%, indicating that location and stand-specific variables had a greater impact on damage probability than the dominant species, which agrees with previous studies, in which variable species composition was not significant (e.g., Hanewinkel et al., 2014). Nevertheless, there were differences between a stand's dominant species and the behaviour of certain other variables, such as diameter or height. It is important to note that species-specific characteristics, such as rooting or crown shape, were not included as specific variables in our models, and were assumed to be implicitly represented by the variable that differentiated plot species dominance (Albrecht et al., 2010) as a necessary simplification due to lack of more detailed records. This induces us to consider that the seemingly lack of importance of this variable does not mean that there are not differences, but rather that the species-based differences related to tree's architecture should be more specifically considered, whenever available data make it possible.
Finally, there are at least four potential applications for the models developed in this study: 1) the production of landscape-level hazard maps for snow and wind damage; 2) identification of the key factors controlling damage occurrence; 3) evaluation of risk for individual stands; and 4) the simulation of stand vulnerability to damage occurrence under variable conditions. Combined with the probability of stand damage proportion (Díaz-Yáñez et al., 2017), these can be used to estimate optimal silvicultural treatments. Despite the limitations of the available data, and the method considered, we regard the present models as a promising approach to analyse dynamic interaction of disturbances by snow and wind and forest structure, as affected by management interventions or natural changes.

Declaration of Competing Interest
The authors declare that they have no conflict of interest.