Developing predictive models of wind damage in Austrian forests

Abstract

• Context

Among natural disturbances, wind storms cause the greatest damage to forests in Austria.

• Aim

The aim of this study is to quantify the effects of site, stand and meteorological attributes on the wind disturbance regime at the operational scale of forest stands.

• Methods

We used binomial generalized linear mixed models (GLMMs) to quantify the probability of damage events and linear mixed models (LMMs) to explain the damage intensity at the forest stand level in four management units with a total forest area of approximately 28,800 ha.

• Results

Timber stock volume, stand age, elevation, previous disturbances, wind gust speed and frozen state of soil contributed in explaining probability of wind damage. While the model of disturbance probability correctly classified 90 % of all cases in the data set (specificity 95 %, sensitivity 26 %), the model for damage intensity explained only low percentages of the variation in the observed damage data (full model R ² = 0.38, fixed effects-only model R ² = 0.09; cross-validation in the four forest management units yielded similar R ² values).

• Conclusion

The developed models indicated that decreasing the proportion of Norway spruce (Picea abies [L.] Karst), limiting stand age and reducing the timber stock in course of tending treatments in stands exposed to wind disturbance can mitigate the risk and the expected damage intensity. High gust speeds and salvage cuts after earlier damage increase the probability of further wind disturbance events.

1 Introduction

In Europe, during the period between 1950 and 2000, an average of 18.7 mill. m³ of timber were damaged by wind annually (Schelhaas et al. 2003). This makes storms leading to uprooting and stem breakage of trees the most detrimental natural threat to European forests (i.e. approximately 66 % of total damage from storms, fire, bark beetles and snow). In Austria, post-windthrow salvage logging fluctuated between approximately 1 and 10 mill. m³ of timber per year in the period 1990–2012, which corresponds to shares of 4 to 50 % of the annual cut (Prem and Beer 2012; Anonymous 2013). Peak years of windthrow salvage due to large-scale stand-replacing events mainly during the winter season were 1990 (7 mill. m³), 2007 (9 mill. m³) and 2008 (10 mill. m³). Beside these severe storm events, a high proportion of the timber salvage was due to small-scale or low intensity disturbance events. However, these less intense disturbance damage events accumulate to substantial losses in timber value and cause additional costs for harvesting and further follow-up costs regarding planting, tending and other silvicultural measures. Also, after large-scale stand-replacing events, management plans become obsolete and need to be updated. Beside the adverse economic consequences in timber production, windthrow can negatively affect other forest ecosystem services like protection against rockfall and avalanches (Brang et al. 2006), drinking water preservation (Weis et al. 2006) or in situ carbon sequestration (Thürig et al. 2005).

The wind disturbance regime is driven by the interplay of forest characteristics and weather (Dale et al. 2000). Species composition, stand height, stand edges, canopy roughness and tree attributes such as crown length and slenderness correlate with wind damage (e.g. Valinger and Fridman 1999; Mitchell et al. 2001; Olofsson and Blennow 2005; Sellier and Fourcaud 2009). Beyond a wind speed of 45 m/s, stand-replacing damage is almost certain, regardless of stand condition (Gardiner et al. 2010). With decreasing wind speed, the effects of tree and stand characteristics on damage intensity become more apparent (Xi and Peet 2011).

It is expected that frequency of storm events in Central Europe may increase in a warmer climate (Lindner and Rummukainen 2013); however, uncertainty is high (Stocker et al. 2013), and storm tracks may shift northwards resulting in new areas exposed to wind damage. An intensifying wind disturbance regime may exert a positive feedback on bark beetle disturbances through the provision of abundant breeding habitat (e.g. Marini et al. 2013). Additionally, a warmer climate may benefit bark beetles, which may then complete two or even three life cycles per year (Jönsson and Bärring 2011). Other disturbance agents like snow breakage and regular harvests can also modify the structure of forest stands and increase their susceptibility to wind disturbance. Because of the magnitude of potential forest loss and economic relevance, interest has been growing to identify stand and site attributes which explain the variation in damage and to develop quantitative models to assess the vulnerability of forests to wind damage as a prerequisite for targeted risk management. The literature of storm damage in Central European forests is extensive, and many studies scrutinized various driving factors of large-scale storm damage (e.g. Dobbertin 2002; Schütz et al. 2006; Schindler et al. 2009; Schmidt et al. 2010). However, intermediate and small-scale endemic wind disturbances are less widely researched, although their cumulative effect can be significant on forest ecosystem services (e.g. Nagel and Diaci 2006; Klopcic et al. 2009).

Major data sources for such studies include salvage records kept by forest enterprises (e.g. Hanewinkel et al. 2008; Klopcic et al. 2009), regional to national scale damage statistics either based on large-scale forest inventories (e.g. Jalkanen and Mattila 2000) or semi-quantitative salvage reporting schemes on the administrative district or province level (e.g. Thom et al. 2013). Recently, the use of damage estimates derived from remote sensing information (e.g. aerial photographs and satellite images) has attracted much attention (e.g. Lanquaye-Opoku and Mitchell 2005; Usbeck et al. 2012). However, each of these approaches has some limitations. Standard management records provide local operational context and report damage also of a few cubic metres of timber only; however, detailed information on stand and site variables is generally missing, and as a consequence, variables that can explain the damage are scarce. Spatial coverage is usually limited, as book-keeping rules vary greatly among forest enterprises. On the other hand, a great advantage of large-scale forest inventories is large spatial coverage. Operational context of individual inventory plots is missing, however, similarly to regional damage monitoring systems. These rely on some kind of qualitative damage assessment in a highly aggregated form. Considering remote sensing, it has been used mainly to assess damage extent after large-scale events so far (e.g. Schindler et al. 2012). Beside these observational approaches, there has been experimental work on tree pulling (e.g. Nicoll et al. 2006) and mechanistic modelling to determine critical wind speed for either uprooting or stem breakage and then calculating the probability of the occurrence of such wind speed by assessment of the local wind climate attributes (e.g. Peltola et al. 1999; Gardiner et al. 2000). While the latter approaches provide a clearly defined link to weather phenomena and are thus potentially applicable for climate change impact assessments, most empirical studies contain only local relationships without general transferability to other regions or conditions.

In this study, our objective was to develop quantitative statistical models to estimate (1) the probability of wind damage events and (2) the intensity of the damage in Eastern Alpine mountain forests. We employed a large empirical database considering forest and site characteristics, weather data, forest management and other disturbance agents. The focus was at the stand scale due to its importance for operational forest management.

In particular we hypothesized that

• (a) by utilizing data with huge spatial coverage, established empirical relationships are robust over a wide range of conditions;

• (b) the use of weather-related predictor variables improves model performance and reliability.

2 Material and methods

2.1 Study area

Data from four management units (FMU) of the Austrian Federal Forests (AFF) were available for the current analysis. AFF is responsible for the management of 588,000 ha of forest in Austria (i.e. 15 % of total forest area of the country). The four FMUs were Traun-Innviertel, Steyrtal, Waldviertel-Voralpen and Steiermark (Fig. 1); for the current study, 28,870 ha of forest were considered.

Fig. 1

Location of the four forest management units under study in Austria

Elevation ranges from submontane (400–800 m a.s.l.) to high montane (approximately 1,600 m a.s.l.) vegetation belts (Table 1) in all the FMUs except in Traun-Innviertel, which is exclusively located at elevations below 700 m. Norway spruce (Picea abies [L.] Karst) is the most abundant tree species in these AFF management units, a tree species which is considered to be one of the most vulnerable to wind damage in the region (Hanewinkel et al. 2013). Other conifers are present in the area, but with a much smaller percentage (Abies alba Mill., Larix decidua Mill., Pinus sylvestris [L.]). The main broadleaved species is Fagus sylvatica [L.]; other broadleaves have just minor proportions of up to 4 % of basal area. All four FMUs have a uniform age class distribution up until the usual rotation length of the main species (100–140 years). Age classes above that show a strongly decreasing trend, with the oldest stands being approximately 220 years old.

Table 1 Characteristics of four Austrian Federal Forests management units which provided data for the analysis

Forest area within the FMUs is structured into several administrative levels. Districts usually have an area of a few thousand hectares. A district contains several dozen compartments, an administrational unit that—in a mountainous environment—confines an area with similar site properties within major topographic borders (27.1 ha on average in the study area). The smallest operational unit is the sub-compartment (i.e. the forest stand), which usually has an area of several hectares (6.6 ha on average in the study area). Sub-compartments have quite homogeneous site and stand conditions and are the basic silvicultural planning and treatment unit. In Austria, 10-year management plans include operational silvicultural prescriptions at stand level and yield regulation at the district level (i.e. determination of allowable annual cut).

All four FMUs suffered great losses of timber due to storms in 2007 and 2008, with proportion of stands being damaged as high as 28 % (FMU Traun-Innviertel in 2008). FMU Steyrtal also experiences substantial storm damage in 2003. See Fig. 2 for a detailed comparison of wind salvage and total timber removals in the FMUs in the period 2003–2008.

Fig. 2

Average harvested timber volumes over all stands in the four forest management units of the study area (black column, total timber removals including salvage; striped column, wind salvage)

2.2 Database

To build the database for the analysis, stand and site data of management plans of 15 districts within the four FMUs were combined with related harvest records and a gridded weather data set covering the FMUs under study (Table 1).

2.2.1 Forest management plans

The management plans covered the decade from 1999 to 2008. For all stands, the plans described attributes including yield class of the most abundant tree species in a stand, volume and age by species and provide a qualitative description of the mixture type. However, the plans did not include details on compositional and structural features. Since management plans describe the initial state of stands at the beginning of the respective 10-year planning period only, the annual development of stand attributes (in our case timber stock volume at tree species level) over time was projected by means of yield tables (Marschall 1975) and removals reported in the harvest records. Stands younger than 20 years of age were not included in the analysis database, as such young stands were not considered to be vulnerable to wind damage.

2.2.2 Harvest records

Harvest records contained all timber removals specifying the year of harvest (without exact date), total extracted volume and the reason for the removals distinguishing regular harvests and salvage due to various damaging agents. However, no indication of the spatial distribution of the fellings inside the sub-compartments was included in these records. Regular harvests and salvage due to wind (no differentiation between uprooting and stem breakage in the records), snow and bark beetles were each cumulated for periods of up to 4 years prior to any year in the 10-year planning period to account for damage history in the stands. This reduced the length of the 10-year time series of damage data to 6 years. The reported harvested volume was multiplied by 1.2 when related to standing stock to account for standard practices with regard to treatment of harvest residues (Pretzsch 2010).

2.2.3 Weather data and related information

The assessment of the relationship between damage and weather data including storms necessitates local time series of weather variables at the level of the investigated forest stands at a temporal resolution that allows the identification of the driving weather stimuli. For the current analysis, air temperature and wind speed data were provided for all four FMUs. Daily time series of air temperature (minimum, mean, maximum), on a mesh with a width of 100 m over the FMUs, were generated. The data were interpolated from the network of weather stations of the Austrian weather service (ZAMG; http://www.zamg.ac.at/cms/de/klima/messnetze/wetterstationen). The entire network consists of more than 200 automated stations located throughout Austria. A second-order polynomial fit based on the four seasons that tracks vertical temperature gradients was applied to capture the behaviour of the air temperature field within the complex Alpine topography. The interpolation routine distinguished between three regions covering the study FMUs and depended on longitude, latitude and elevation. As no gap filling was applied to the observations, the interpolation relied on the original measurements. The number of weather stations used varied between 15 and 20 depending on year and FMU. Based on Paul et al. (2004), daily air temperature was used to calculate whether the soil was frozen any day of the year. If the uppermost 10-cm layer of the soil was calculated to be frozen, then it was assumed to have a stabilizing effect on the trees against windthrow. Calculating the soil temperature was based on mean annual and summer air temperature and the minimum and maximum air temperature of the current day. Leaf area index, understory vegetation and litter mass of the soil were also taken into account, for which average values of the study area were used defined by expert knowledge. For more details on the soil temperature model, see Paul et al. (2004).

Wind is perhaps the most difficult weather parameter to be generated on a grid, especially in a complex orography as the European Alps. One obvious reason is that wind measurements are representative for only a very small area within which the measurements are taken. The highly discontinuous propagation of the wind field in space makes it almost impossible to homogenize observed time series data by comparing them to other series, farther away. Another inherent problem is that wind measurements carried out at one station are nonhomogeneous in time, as any change in the roughness length of the surrounding topography which may be caused by a growing tree has a differing impact on the measurements. So, wind observations are fraught with problems, and hence, it is difficult to interpolate measurements in space and time.

In the current study, we used integrated nowcasting through comprehensive analysis (INCA) to provide wind data for the study FMUs (Haiden et al. 2011). INCA uses digital elevation data of 1 × 1 km grid size. In the case of wind, the nowcast starts with a three-dimensional analysis based on a first guess obtained from a numerical weather prediction (NWP) model output that is enhanced by the consideration of further observations at weather stations. The wind fields are calculated by transforming 10-m wind observations to the NWP model-level wind using an elevation-dependent factor and by applying an inverse distance squared interpolation routine on the observed corrections. Additionally, an iterative relaxation algorithm is enforced to warrant mass-consistent fields. Wind vectors at grid points near to stations are kept at the observed values during the relaxation procedure. Thus, the INCA data set has been designed to match the observed values. In this study, maximum daily wind speeds were computed from 24-hourly wind speeds (10-min wind speeds at full hours) covering the period 2003–2008. Ten-min wind speeds were transformed to 2-s gust speeds (VMAX) using multiplication by a gust factor of 1.65 (Cvitan 2003). Such short-term gusts are commonly considered as major determinant of wind damage in forests (Mayer 1987). To visualize the spatial heterogeneity of the wind speed data, the number of days with a 2-s gust speed above 30 m s⁻¹ is shown for the study FMUs in Fig. 3.

Fig. 3

Number of days with daily maximum 2-s gust speed [m s⁻¹] exceeding 30 m s⁻¹ during the period 2003–2008 for the four forest management units

Weather-related data were linked to forest compartment centroids using GIS software (ESRI 2012). All stands within a forest compartment were assumed to have the same weather attributes. All stand-level attributes available for the analysis are shown in Table 2.

Table 2 Available stand-level attributes for modelling wind damage in the study area

2.3 General modelling approach

The modelling process was structured in two main steps. First, modelling the probability of a wind damage event, then as the second step, the damage intensity given that a damage event had occurred in the stand. A mixed model framework was used because of the hierarchical structure of the data set. Spatial units (sub-compartments) were nested within bigger units (compartments, districts, FMUs) and carried a 10-year time series data set; hence, spatial and temporal autocorrelation of model residuals had the potential to cause bias in model estimates (Pinheiro and Bates 2000). With the use of a mixed model framework, random effects were added to the model equations in the model fitting process. These random effects induced a simple correlation structure for observations in the same data cluster (in our case sub-compartment, compartment and district), and therefore, the effects of confounding variables related to these clusters could be estimated (i.e. the random effects) and separated from the actually studied predictor variables (i.e. the fixed effects) during the analysis.

As an exploratory step to attain insight on the interrelationships of variables in the data set, Pearson correlation coefficients among the continuous stand-level variables were calculated. Candidate predictor variables were used in an automated procedure with all possible combinations to fit model equations to the data set for both the probability and the intensity models. Afterwards, Akaike information criterion (AIC) values (Akaike 1973) of the resulting models were compared. This is a method that rewards the goodness-of-fit of a model and penalizes for the increasing number of predictors at the same time. The penalty for the inclusion of one additional parameter was 2 AIC units following Arnold (2010). If further posterior analysis (see below) did not reveal inadequacies, the model with the lowest AIC value was chosen (Burnham and Anderson 2002). Partial effect plots (also called marginal effects or least square means) were used to assess the effect of the predictor variables in the models. Partial effects measure the change in the expected value of the response variable as a result of a change in a certain predictor variable while keeping all the other covariates fixed at the median values of the respective variables in the database (see Table 2).

The robustness of the models was tested by 10-fold cross-validation, in which the data set was randomly broken into ten partitions, and models were fitted to data consisting of all but one partition that served as the test group. This procedure was repeated ten times with a different test group each time; then, goodness-of-fit was evaluated (Mosteller and Tukey 1968). Cross-validation tests were also implemented with the four FMUs being used as partitions.

Versions of the models that included fixed effects only were also used for posterior tests. We did this by multiplying the design matrices of the models by the fixed effects calculated in the model fitting process. This served as an assessment of how well the fixed effect variables could explain wind damage in a general context without using local implicit information related to the four FMUs.

In the model fitting process, predictor variables describing accumulated salvage and regular harvests in the previous 4 years (W_t, B_t, S_t and R_t) were all natural log transformed. Considering these variables, we manually fitted the models with different time periods (t = 1–4 years). Extending the number of the preceding years improved the fit of the models. The 4-year period proved to be the best compromise between increasing model fit and decreasing the number of years available for model development at the same time. In addition, when using 4 years to accumulate the disturbance history of each stand, 6 years remained for model fitting, which met the requirement regarding the minimum level of a random effect, in our case years within the random effect “sub-compartment”. The variance of a random effect cannot be estimated correctly in case of fewer than six levels (Crawley 2002).

In the analysis, different groupings of site type (SITE) were also tested. This aimed at combining sites with similar soil moisture and nutrient supply (as defined in the site classification system of the AFF; Weinfurter 2004). However, none of these groupings improved model fit.

The software package R was used for the statistical modelling (R Core Team 2013). The automated fitting of models for model selection was done with the dredge function from the R package MuMIn (Barton 2013). The plotLMER.fnc function of the R package languageR was used for the partial effect plots (Baayen 2011).

2.4 Modelling the probability of disturbance events

Wind damage probabilities were modelled by logistic regression. This was done in a framework of a generalized linear mixed model (GLMM). Salvage values below 1 % of standing stock were considered registration errors and the respective data record as a nonevent (see Klopcic et al. 2009; Overbeck and Schmidt 2012). To translate the linear predictor of the model to probabilities, Eq. (1) was used

$$ {\uppi}_i=\left({ \exp}^{\alpha +\beta \times {X}_i+{\gamma}_i}\right)/\left(1+{ \exp}^{\alpha +\beta \times {X}_i+{\gamma}_i}\right) $$

(1)

where π _i is the expected annual probability of the occurrence of a damage event in the ith row of the design matrix of the model, α is the intercept, β is the vector of fixed effect parameters, X _i is a row from the design matrix of the model and γ _i is the sum of random intercepts that account for the spatial and temporal cluster effects in the observed damage related to forest stand, compartment and district level. GLMMs were fitted with the lmer function of the package lme4 (Bates et al. 2012).

Linearity of the relationship between damage events and the explanatory variables was assessed by plotting the partial residuals of the model (see Zuur et al. 2008) and fitting smoothed curves using the loess function of the basic R package. Data transformation was used to account for eventual nonlinearities in the data. Classification table (i.e. confusion matrix) and derived parameters such as sensitivity, specificity and the area under receiver operator characteristic curve (AUC) were used for assessing goodness-of-fit of the models. Sensitivity is the ratio of the true positive and the sum of true positive and false negative predictions. Specificity is the ratio of the true negative and the sum of true negative and false positive predictions. Sensitivity and specificity values were plotted against a range of cut-off points which serve as means to differentiate between events and nonevents.

The AUC shows the probability that a randomly selected observed positive event has a higher predicted probability value than a randomly selected observed negative event (Fawcett 2006). It can range from 0.5 (no predictive ability) to 1 (perfect discrimination) and is independent from cut-off values. The somers2 function of the R package Hmisc was used for the calculation of AUC values (Harrell 2012).

2.5 Modelling the intensity of damage events

For modelling damage intensities [m³ ha⁻¹ year⁻¹] by linear regression, data records classified as “damage” (threshold of 1 % salvage rate; see previous section) were used as input in a linear mixed model (LMM) framework. The response variable was natural log transformed to improve normality and homogeneity.

The general equation of the model was

$$ {\pi}_i=\alpha +\beta \times {X}_i+{\gamma}_i $$

(2)

where π _i is the expected annual intensity [m³ ha⁻¹ year⁻¹] of a damage event in the ith row of the design matrix of the model, α is the intercept, β is the vector of fixed effect parameters, X _i is a row from the design matrix of the model and γ _i is the sum of random intercepts that account for the spatial and temporal cluster effects in the observed damage related to forest stand, compartment and district. LMMs were fitted with the lmer function of the package lme4 (Bates et al. 2012).

Goodness-of-fit was evaluated by the coefficient of determination (R ²). It indicates how much variation in the data is explained by the fitted model. We calculated root mean squared error (RMSE) to see how close predictions were to observed values. Normality and homogeneity were tested by plotting residuals against predicted values. Histograms were used to assess normality of residuals. Residuals plotted against explanatory variables were inspected for linear relationships. Test indices were also calculated for models with fixed effects only.

3 Results

3.1 The disturbance regime

For the 6,015 stands that were available for the analysis during the 6-year period, the ratio of events to nonevents (i.e. years with a wind damage versus years without a damage) was approximately 1:16. With regard to the wind damage intensities, 43 % of the disturbance events caused damage smaller than 10 m³ ha⁻¹, 46 % between 10 and 50 m³ ha⁻¹ and 11 % more than 50 m³ ha^-1. The mean intensity of the wind damage events was 25.9 m³ ha⁻¹. Table 3 shows the proportion of wind disturbance events which were preceded by other disturbance or regular harvest in the previous 4 years in the same stand.

Table 3 Proportion of wind disturbance events in the period 2003–2008 with at least one salvage cut (wind, bark beetles or snow) or regular harvest occurring in the same stand in the four preceding years

3.2 Probability of disturbance events

The final model included predictor variables ELEV (elevation), AGE (mean stand age), VOL (timber stock volume before harvests of actual year), W4, B4, S4, VMAX and SF (soil state on the day with the highest 2-s gust speed; 1, frozen; 0, not frozen). All other variables were found uninformative. YC (yield class of the main species of the stand) was excluded from the analysis to avoid multicollinearity (see Zuur et al. 2008), as it strongly correlated with several other explanatory variables (e.g. Pearson correlation was 0.70 with VOL, −0.56 with ELEV). Adding interaction terms to the models did not improve the goodness-of-fit.

The signs of parameter estimates of the predictor variables were all positive, except for SF (Table 4). Figure 4 presents the effects of individual predictors on the estimated probability of a wind damage event. The AUC value of the model was 0.84. When used with the fixed effects only, the AUC was 0.71, which is considered “acceptable” in statistics literature (Hosmer et al. 2013). Partial residual plots showed that the assumption of linearity was valid for all the candidate predictor variables (not shown here). In the 10-fold cross-validation of the selected model, AUC values varied between 0.68 and 0.76, thus indicating good stability. When the four FMUs were used separately as test data, AUC values were 0.75, 0.72, 0.70 and 0.80 for Traun-Innviertel, Steyrtal, Waldviertel-Voralpen and Steiermark, respectively. AUC values in the current study were in line with other studies of storm damage in Central Europe in which this index was used to evaluate the classification into damaged and undamaged stands (AUC = [0.78–0.79] in Schindler et al. 2009; AUC = 0.76 in Klaus et al. 2011; AUC = [0.73–0.74] in Schindler et al. 2012).

Table 4 Estimated model coefficients of the generalized linear mixed model for the probability of occurrence of wind damage and the linear mixed model for damage intensity [m³ ha⁻¹ year⁻¹]

Fig. 4

Partial effects in the wind damage probability model. The Y-axis denotes the annual probability of a wind damage event. The solid lines represent the partial effects of the predictor variables (X-axis), the starting and end points denote the minimum and maximum values in the database used for model development. For a description of variables, see Table 2

The cut-off point analysis (Fig. 5) for the fixed effects-only model showed a proportion of 65 % correctly classified cases for both the disturbance events and nonevents at a cut-off value of 0.04 (i.e. both sensitivity and specificity had a value of 0.65 at this cut-off). Nevertheless, it may be more practical to take a look at these values at a cut-off value, where the numbers of predicted events and nonevents are closest to the observed ones (i.e. disturbance frequencies in the observations and predictions are similar). This cut-off value was 0.13, and the related sensitivity and specificity were 0.26 and 0.95, respectively (the related share of overall correct classifications was 90 %). This means that nonevents were predicted with good reliability, and in case of damage events, the model identified correctly one fourth of the observed events.

Fig. 5

Cut-off values and the related sensitivity, specificity and classification rate values for the probability model of wind disturbance events (predictions made without random effects)

To assess the importance of weather-related predictor variables in the models, we also fitted the final model in a version where these were excluded. When omitting VMAX and SF from the probability model, AUC values decreased slightly, to 0.83 in case of fixed and random effect model and to 0.70 in case of the fixed effect-only model. Cross-validation indicated only a small decrease in robustness (not shown here).

3.3 Intensity of disturbance events

When modelling wind disturbance intensity, many of the predictor variables were found uninformative in the model selection process, except for PA (proportion of Norway spruce), VOL and W4 (Table 4). VOL showed higher partial effect than the other predictor variables (Fig. 6). The linear relationship between the predictor variables and the response variable was confirmed by the residual plots. Residuals were normally distributed. Normality and homogeneity improved substantially by ln-transforming the response variable. Considering goodness-of-fit, there was a big difference between R ² values of the full model and the one refitted with fixed effects only (Table 5) with R ² values of 0.38 and 0.09, respectively. Related RMSE values were 38.49 m³ ha⁻¹ year⁻¹ and 43.77 m³ ha⁻¹ year⁻¹, respectively.

Fig. 6

Partial effects in the wind damage intensity model. The Y-axis denotes the intensity of a wind damage event. The solid lines show the partial effects of the predictor variables (X-axis), the starting and end points denote the minimum and maximum values in the database used for model development. For a description of variables, see Table 2

Table 5 Performance statistics of the generalized linear mixed model for the probability of occurrence of wind damage and the linear mixed model for the damage intensity [m³ ha⁻¹ year⁻¹]

In the 10-fold cross-validation, the fixed effect-only model had R ² values between 0.02 and 0.20. In case the four FMUs were used separately as test data, R ² values were 0.08, 0.09, 0.08 and 0.05 for Traun-Innviertel, Steyrtal, Waldviertel-Voralpen and Steiermark, respectively.

4 Discussion

4.1 Disturbance drivers and model quality

Classification performance tests showed that the wind damage probability model performed moderately well also in case of a “fixed effects-only” version. This indicates that the model could be used for prediction purposes outside of the study area (see Bolker et al. 2009). Model results clearly indicated a higher probability of wind disturbance for old, highly stocked forests at higher elevations in years that were preceded by high intensity disturbances (see Fig. 4). The inclusion of AGE may represent the positive effect of increasing tree height on damage probability, as the correlation between AGE and VOL was fairly small. ELEV may capture the combined effects of increasing wind speed with increasing elevation and increased exposure to winds at higher elevations due to reduced shielding effects of surrounding topography.

The strong effect of preceding disturbance events may be due to newly created stand edges and increased canopy roughness (Rochelle et al. 1999; Schütz et al. 2006), which is in line with findings of other studies (e.g. Klopcic et al. 2009; Thom et al. 2013). W₄, B₄ and S₄ in our analysis had very similar effect sizes in the case of the probability model (see Table 4), which means that changes in stand structure induced by salvage of wind, bark beetle and snow damage are apparently similar and that they may be interchangeable (they also have the same measurement unit) in the model. Therefore, they could also be used in an aggregated form when estimating wind damage probability. Intensity of individual salvage cuts was in most cases below 20 % of standing volume and thus at levels similar to regular thinnings and regeneration cuts. The fact that regular harvests (R₄) were not found important in explaining wind damage indicates a certain difference between the effects of regular harvests and salvage cuts. Whether this is the result of careful forest management or some other factor, could not be tested in the current study.

The four FMUs in the current analysis represent a broad range of sites and stands in Austrian commercial forests. Salvage practices may vary among regions, which affects the relationship of stand and site conditions and management activities on one hand, and damage from wind on the other. The robustness of the developed models as revealed by the cross-validation among FMUs indicated a similar management strategy in all four units of the AFF under study.

Overall, standard predictors describing forest site yielded only small effects on damage probability (only ELEV was found to be informative). This may be explained by a rather low representativeness of site descriptors for entire stand polygons due to substantial small-scale variation of site and soil conditions in mountain forests. Moreover, the spatial variability of site attributes was apparently much larger than that of wind damage, reducing their explanatory power. In the light of these arguments, the inclusion of SF (indicating frozen soil) in the model highlights the importance of the interplay of soil and weather conditions in determining susceptibility to wind damage. However, SF also highlights the challenge of including attributes in model development which are subject to substantial uncertainty, regarding both space and time dimensions. With SF, this is particularly true to the required match of VMAX and unfrozen soil status. Related uncertainty would have even increased further in the case of including a stand-specific soil water balance calculation which requires spatially accurate estimates of water holding capacity of the soil.

The small effect size of VMAX in the probability model and its exclusion from the intensity model confirmed the hypothesis that incorporating wind speed as a driving factor in our models was a challenge. Analysing whether the relatively coarse spatial resolution of the wind speed data (1 × 1 km) or the inherent variability in storm damage events was the major reason for the weak effect of VMAX in the probability and intensity models was beyond the scope of the current study. An alternative to using wind speed data is describing the local topographic exposure of forest stands to wind. This method is widely employed in storm damage research and uses indices like topex, topex-to-distance or other indices that include information on the aspect and slope of stands (see e.g. Quine and White 1998).

There are other empirical studies of storm damage that link wind field characteristics to wind damage (e.g. Schütz et al. 2006; Schindler et al. 2009, 2012). However, they do this for large-scale singular storm events, where the exact timing of a storm is known, and linking the gust speed estimate to the damage event is less problematic. In our case, excluding the weather-related variables (VMAX and SF) from the probability models caused only a small decrease in explanatory power. Nevertheless, the sign of both VMAX and SF and the effect size of SF in the probability model clearly indicated conceptual consistency and the relevance of including explicit weather-related attributes in wind disturbance models.

In case of the intensity model, the big difference of R ² between the full model (fixed and random effects) and the fixed effects-only model indicated that the variables used in the study explained only a small proportion of the variance in wind salvage and that implicit local effects could not be revealed with the available database. This may be due to features of the local wind climate or stand conditions which are not accounted for by the available data. Local adaptation of trees to higher wind speeds can decrease wind disturbance susceptibility (Nicoll et al. 2008). With regard to forest stewardship-related attributes, Thom et al. (2013) found that attributes such as ownership or road density contributed significantly to explaining the variation in damage data at the level of administrative districts.

4.2 Limitations of database and study design

The current study used data from forest management plans. It is important to note that originally, the underlying data had not been compiled for scientific purposes, and therefore, just the typical standard set of site and stand characteristics was available for the analysis. Recording errors of timber removals can lead to inconsistent data sets. For instance, the assignment of removals to one of several causes is such a crucial issue. The important advantage of the database is its spatial coverage and that the same data collection procedures have been used throughout all the studied FMUs.

Linking a disturbance event to a highly stochastic predictor such as gust speed may impose several problems. For the analysis, exact dates of the disturbance events were not available. In case of windthrows, it is possible that a winter storm damages trees in November or December, but the event only gets registered and the timber salvaged during the next year. This can lead to cases in the model matrix when high gust speeds (>30 m s⁻¹) are linked to no damage (compare Nilsson et al. 2007) in the stand, or high amounts of damaged timber to very low gust speeds (e.g. <5 m s⁻¹).

The topography used by INCA is given on a grid with a resolution of 1 km. As such, the terrain model may fail to approximate highly complex structured topography sufficiently well. However, since there is a lack of observational data, the INCA data set provides the nearest approximation of the wind conditions in the study regions over a longer period. Further wind downscaling techniques combined with in situ measurements could improve the analysis of the actual wind conditions within the forest areas. Despite these limitations, we intended to include VMAX in the analysis to test if it contributes to explaining damage from wind disturbances. The relevance of developing predictive models sensitive to changes in wind-related drivers becomes evident when analysing climate change impacts on forests and related ecosystem services (e.g. Lindner et al. 2010). The effect size of VMAX in the probability model was relatively small, and it was not even included in the intensity model. However, with a higher spatial resolution of accurate wind speeds and an exact recording of the dates of wind salvage, empirical models for wind damage prediction using wind speed as a predictor could very likely be improved significantly.

Another issue related to the INCA wind speed data was the availability of data from 2003 onward whereas the disturbance damage data was available for the period 1999–2008. However, no straightforward approach was available to fill this gap in the wind speed data. In addition, 4 years from the 10-year data record were used to build predictor variables characterizing harvest and damage history of stands. Nevertheless, these predictor variables proved to be important in explaining damage events, similar to the findings of Thom et al. (2013), who used salvaged damage from the preceding 2 years as explanatory variable in a recent study of the wind and bark beetle disturbance regime at landscape scale.

An important aspect of the disturbance regime, the spatial extent of disturbance events inside a stand, or at higher aggregation level inside a compartment, could not be considered. Substantial efforts and additional data (mainly deducible from remote sensing sources) would have been needed to make the database spatially explicit. In our models, effects of damaged neighbour stands were not included as fixed effects for similar reasons. However, the latter practice would also have been contrary to our aim to explore how well predictions could be made from stand-level information only.

5 Conclusions

The developed models indicate that decreasing the proportion of spruce, limiting stand age and reducing timber stock in course of tending operations in stands vulnerable to wind disturbance are options to mitigate the risk and the expected damage intensity. Furthermore, results clearly indicated that previous disturbances from wind, bark beetles and snow increase the risk of wind damage. This, therefore, underscores the importance of stand stability and resilience in forest management.

It was found that standard data from management plans have the potential to be used to identify stands at high risk of wind damage. However, expectations that salvage volumes due to wind disturbance damage at the stand level can be accurately projected from standard data appear too optimistic.

The big difference in performance between models with fixed and random effects and models without random effects highlighted the importance of unknown attributes that have yet to be discovered in order to improve wind disturbance modelling. The limited improvement in model performance by including wind speed data points was a weakness on one hand, but indicated that there is potential to further improve empirical wind disturbance models by improved matching of damage events to weather data.