Winter storm risk assessment in forests with high resolution gust speed data

Winter storms pose a major threat to forest management in Central Europe. They affect forests at a large spatial scale and produce large losses in standing and merchantable timber within few hours. The assessment of winter storm vulnerability by statistical modelling serves as an important tool to tackle uncertainities about the damage risk and to inform management decision processes. This study made use of an extensive forest inventory data set from South-West Germany before and after winter storm Lothar in 1999, one of the most severe storm events in Germany over the last decades. Hierarchical logistic models were fitted to relate storm damage probability on individual tree level to features of dendrometry, site, orography, and storm-specific high resolution data of maximum gust speed. We developed two different approaches to implement gust speed as a predictor and compared them to a baseline model with a structured spatial effect function with no gust speed information. Regional and local variability which could not be described by the predictors was modelled by multi-level group effects. Generalisation performance was tested with a spatially and temporally independent data set on storm separation between explicit spatial gust speeds and unknown variability achieved with the parametric multi-level approach led to a higher degree of transparency and utilisability.


Introduction
Windstorms play a dominant role in the disturbance regime of European forests (Seidl et al. 2014). Especially in the season between November and March, extratropical cyclones frequently impact large areas with heavy rainfalls and strong, energetic wind gusts. Exposed to this threat, forest stands were often severely damaged or replaced in a whole region within few hours. The risk of damage confronts forest owners with a large degree of uncertainty in their management decisions. Hence, an assessment of the vulnerability of their forests to winter storm damage would help to better adapt their goals and actions to mitigate future loss .
For vulnerability assessment, it is crucial to understand and quantify the underlying mechanisms which lead to tree damage. However, there are manifold biotic and abiotic factors that are known to have an effect on storm damage (Gardiner 2021). For instance, individual tree features like crown area, root architecture, and height-to-diameter ratio, but also soil rooting depth and exposition of the location do impact storm damage vulnerability. The effectiveness of these factors is hard to grasp because they interact in various complex ways, on different spatial scales, which often result in non-linear behaviour (Messier et al. 2016). Here, modelling could help to break down this complexity by identifying Communicated by Thomas Knoke. the most meaningful factors and describing their effect on tree damage by storms in a plausible way. In their application, well-parameterised models are able to generalise these effects to new settings by robust predictions of storm damage vulnerability for forest stands at different locations in different climatic conditions. Forest storm damage models generally follow a mechanistic or a statistical approach. Mechanistic models aim to describe the physical properties of a wind-tree system. In particular, they describe how effective forces and the reactions of trees result in either windthrow or stem breakage. Mechanistic models allow a close representation of the cause-effect-relationship of tree damage, but require substantial parameterisation effort for specific tree species in a region of interest. For example, a prominent mechanistic model ForestGALES (Gardiner et al. 2000) is parameterised for forests of Great Britain (Hale et al. 2015), Finland (Peltola et al. 1999), or Canada (Byrne and Mitchell 2012). Its application is restricted to a few stands in Germany with comparable species composition and growth conditions. An alternative approach, statistical storm damage modelling, is based on empirical data typically derived from forest inventories and remote sensing. It relates the observed outcome (e.g. a binary variable tree undamaged / damaged) to a set of explanatory variables by the statistical properties of the data. The most prominent algorithms for statistical storm damage modelling were, for example, logistic regression (Schmidt et al. 2010;Kamimura et al. 2016;Suvanto et al. 2019) and machine learning classification techniques (Hart et al. 2019;Albrecht et al. 2019;Suvanto et al. 2019). Statistical models perform well within the range of their training data. Hence, a large data basis that covers wide gradients of all important factors generally leads to more robust effect estimates and to a better generalisation.
For their application, models ideally contain predictors, which are either directly or indirectly affected by management decisions. This allows for linking management options with damage probabilities (i.e. vulnerabilities) and, hence, supports forestry decision making. Common predictors in statistical storm damage models were derived from features of individual trees, stands, site, and orography . Schmidt et al. (2010) showed plausible and meaningful effects for tree species, dendrometric variables, soil water regime and topographic exposition with respect to the main wind direction. In this way, management options such as the choice of site-adapted tree species and target diameter can be attributed to storm damage probabilities. Nevertheless, wind speed as a predictor for storm damage in forest, is less frequently found in statistical modelling. Gusts, in particular, are the direct physical forces effecting tree crowns and leading to their drag, bending, swaying, and, in the end, to stem breakage or uprooting. Due to the mechanistic nature of this interaction, gust speed has commonly been included in mechanistic models (Gardiner et al. 2008;Kamimura et al. 2016). Statistical storm damage models that can be used for robust predictions usually cover large areas and, hence, a continuous spatial representation of gust speed is needed for their parameterisation and application.
For storm damage modelling, gust speed fields must be spatially highly resolved ( < 100 m × 100 m ), accurate, and temporally explicit. However, these properties are not present in the commonly available gust speed data. The wind speed measurement network of weather services is much too sparse and thus does not represent the high variability of gust speed. The rather coarse spatial resolution of reanalysis models (e.g. 0.25 • × 0.25 • for ERA5 (Muñoz-Sabater et al. 2021)) is also unsuitable for statistical storm damage modelling. Besides, it is problematic to use mean wind speed from wind atlases as a predictor because each damage-inducing storm has its specific storm track. Due to these shortcomings, Jung and Schindler (2019) developed a winter storm atlas of the most severe winter storms in Germany from 1981 to 2018 (GeWiSa). GeWiSa includes highly resolved ( 25 m × 25 m ) maximum gust speeds for the 98 most severe storms.
In the absence of an explicit representation of gust speed, Schmidt et al. (2010) used a spatially structured effect (Brezger and Lang 2006) in the predictor set to model local tree damage emergence. This proxy variable most likely describes the spatially autocorrelated gust speed among other, unobserved factors. More recent studies use fine resolution gust speed fields to employ the direct agent in storm damage models. For example, Albrecht et al. (2019) showed moderate to good validation results for models with storm event-specific gust speed data as the only predictor (for single tree damages after storm events Vivian/ Wiebke 1990 and Lothar 1999 in Southwest Germany). However, when modelled together with a set of other dendrometry-and topography-related predictors, the additional gust speed effect did not improve the model considerably. For Finland, Suvanto et al. (2019) used high-resolution spatial data describing the long-term wind regime to model forest stand damage. Although they used rather general 10-year return-rates of gust speeds for parameterisation, the variable showed a significant contribution with a moderate sensitivity.
Statistical storm damage models usually make use of observations either at stand (Suvanto et al. 2019) or at individual tree level (Schmidt et al. 2010;Albrecht et al. 2019). Stand-level approaches average dendrometric variables like tree height, typically derived from forest inventories, at the cost of information loss. Single-tree models, in turn, need to take into account that observations from a forest inventory plot are not independent from each other. Multi-level models account for such structure in data and allow for the quantification of inter-group variation (Gelman and Hill 2007).

3
This variation includes information about storm damage at group level, for instance damage propagation to neighbours, which is usually either neglected or set aside but rarely used to inform predictions with statistical models.
In this study, we present a statistical approach to assess winter storm vulnerability of forests by modelling damage probabilities of single trees. The motivation was to achieve a high degree of generalisation for applications like decision support systems. For this, a large data set of one of Central Europe's most important storm events of the last decades was used. A focus was set on the explicit integration of a fine-grained spatial representation of near-surface gust speed as the direct driver for damage occurence.

Storm damage data
On December 26th, 1999 winter storm Lothar severely affected large parts of East France, Switzerland, and South-West Germany with maximum gust speed exceeding 50 m s −1 (Jung and Schindler 2019; Wernli et al. 2002). Shortly afterwards, Germany's second National Forest Inventory started and was completed in 2002. We used a spatial subset of this inventory data covering the federal state of Baden-Württemberg, the most important forest damage hot spot of Lothar in South-West Germany (Fig. 1B). The nested systematic sampling took place at locations (tracts) arranged on a 2 km × 2 km grid, each with one to four nested sampling sites (tract corners) arranged on a square with 150 m side length. All tract corners within a forest patch were sampled. Twelve tracts were excluded from the analyses because no gust speed information was available. In total, the data set included 63 117 trees (6613 damaged) on 11 000 tract corners organized in 4223 tracts. Measurements followed the standard German National Forest Inventory protocol (Kändler 2009) with a special record in Baden-Württemberg for individual trees damaged. Tree damage definition included wind throw and stem breakage. This binary variable (no damage or damage) per individual tree served as the response variable in our study.
For validation of these models, we processed a second inventory data set. In order to test the ability of the models to generalise, it was chosen to be spatially and temporally independent from storm event Lothar in South-West Germany. On January 18th, 2018 winter storm Friederike hit Great Britain, North France, the Netherlands, Belgium, and Germany. Although its maximum gust speeds were not as extreme as Lothar's, Friederike led to severe forest damage with a regional focus on the federal states of North Rhine-Westphalia, Lower Saxony, and Hesse in Central Germany . Contrary to the storm event Lothar, no inventory of trees damaged by Friederike was undertaken. Therefore, we assessed storm damage for a region in southern Lower Saxony (Fig. 1C) by interpretation of aerial ortho-photography. Aerial imagery from a flight survey conducted a few weeks after the storm event on February 6 and 13, 2018, was used. The flight survey covered in total 315 000 ha containing a forest area of 120 000 ha . Based on digital color-infrared ortho-photos, all forest patches with a coverage of at least 70% damaged crown surface and with a minimum area of 0.2 ha were delineated manually by visual interpretation. A minimum area of 0.2 ha was chosen a) to reduce the delineation error due to shadow casting, and b) to limit the total mapping effort. Reference orthophotography from 2016 was compared in order to exclude areas which had been damaged and/ or cleared before storm event Friederike. To obtain stand damage data for validation, we processed forest sample plot data at enterprise level from the federal forestry administration Niedersächsische Landesforsten and intersected plot locations with damage polygons. This data set encompassed surveys on a total of 10 900 plots of which 374 intersected with delineated storm damage polygons.

Predictors
The set of predictor variables was taken from the predecessor model described in Schmidt et al. (2010). No further variable selection was performed to ensure comparability.
Tree features with an impact on wind-stress mechanics were included in the model. These were tree species and two dendrometric variables: tree height ( ) and diameter at breast height ( ). Tree species were grouped according to Schmidt et al. (2010): Picea abies (L.) H. Karst., Fagus sylvatica L. and Quercus sp. L., Abies alba Mill. and Pseudotsuga menziesii (Mirbel) Franco, Pinus sylvestris L. and Larix decidua Mill., and other deciduous species. Dendrometric variables at the time of the storm event needed to be estimated by single tree growth models. In the case of event Lothar, and values were taken from Schmidt et al. (2010). Based on a complete record of all and measurements from the first National Forest Inventory in 1987, they predicted and for 1999 with a fully initialized single tree growth model. For the values of the dendrometric variables at event Friederike, we processed forest sample plot data at enterprise level from the federal forestry administration Niedersächsische Landesforsten. The repeated surveys took place in the years 1991-2019. At each plot, was measured for trees within two concentric circles with a radius of 6 m for individuals i with 7 cm ≤ i < 30 cm , and 13 m for was measured for one representative individual per plot, survey, and species. To obtain and at the time of storm event Friederike 2018, we first fitted for each tree species a linear mixed model with response variable , variables and stand as predictors, and as group level. These models were used to predict in 2018 for each tree. Then, mean quadratic diameter ( ) in 2018 was calculated per species and plot. Finally, corresponding tree heights, were predicted by height-diameter curves calibrated with the --pairs of representative trees (Lappi 1997;Schmidt 2009). These means of the dendrometric variables were used since damage information was available at plot level only. Gust speed is supposed to be an important driving force for storm damage on trees. Based on empirical data and orographic features, Jung and Schindler (2019) developed the high-resolution gust speed atlas GeWiSA for 98 severe winter storms in Germany in the period 1981-2018. This data was used to obtain the predictor maximum gust speed ( ). The metric is defined as the maximal wind speed lasting for at least 3 s at a height of 10 m during a storm event.
The topographic situation of a site was depicted by a distance-based topographic exposure index (topex) (Scott and Mitchell 2005). It quantifies the aspect and degree of exposure for each site by measuring the angle to the horizon for a given bearing and distance. Exposed sites show negative topex values, plains have values around zero, and sheltered sites are characterised by positive topex values.
We summed up topex values of two bearings spanning a given angle around the main storm wind direction. This ensures an identical effect of topography for both deviations from the main wind direction. Distance and angle choice followed the model selection of Schmidt et al. (2010). In total, we calculated three topex sums for a distance of 1000 m : 1) pointing in the main direction of the storm (topex ) with bearings spanning an angle of 30 • , 2) pointing away from the main storm direction (topex ) with bearings spanning a 30 • angle, and 3) pointing orthogonal to the main storm wind direction (topex ) with 180 • between bearings. The latter replaced two wide-angled topex sums used by Schmidt et al. (2010), which were autocorrelated. The three pairs of topex sums are suited for a comprehensive description of orography with sensitivity to an assumed wind direction. Topex sums describe exposed and sheltered locations. In combination with topex sums, topographic elements like slopes, saddles, and ridges perpendicular to the main wind direction are distinguished. Topex sums accentuate ridges and valley floors in the direction of the storm. All topex sums were calculated on basis of the Copernicus Digital Elevation Model Cop-DEM-GLO90 (Release 2020) (Copernicus 2020).
Information on soil water characteristics was included as a proxy for different rooting zone conditions and, hence, the ability of a tree to anchor in the soil. This variable was derived from the water regime category of the forest site mappings for Baden-Württemberg and Lower Saxony. Forest site mapping assesses pedological, geomorphological, hydromorphological, and geobotanical features of site and soil. In Germany, forest site mapping of the federal states varies in methodology and units, and, hence, the level of detail. In the parameterisation, data from Baden-Württemberg soil water regime was categorized in 28 units grouped in three classes: terrestrial, ground-water influenced, and waterlogged. In order to reduce the number of parameters we pooled these 28 units into storm-specific categories by effect comparison and expert knowledge. In Lower Saxony soil mapping units encompassed a total number of 43 major water categories differentiated into up to three subcategories. These units were manually reassessed and attributed to the pooled storm-specific groups.

Modelling
We followed two different approaches for modelling damage by storm event Lothar using the training data. The first, SSE, was a generalised additive model (Bernoulli distributed, logit link-function) with a spatially structured effect as described in Schmidt et al. (2010). It was fitted as a benchmark model which did not include any explicit gust speed information. For each tree individual i, the binary response after storm event Lothar, where is a binary indicator, was modelled by a logistic regression as follows: Damage probability p i as the only parameter of a Bernoulli random variable was modeled by a set of individual-specific covariates x i and its respective regression parameters . The covariates were: tree , log ( ) , and log ( ) in interaction with , , and the three topex sums, namely , , and . There was no valid information about wind speed available for Schmidt et al. (2010) so that it was neccessary to assume spatial location to be a valid proxy for maximum gust speed, i . For that f i , i , a two-dimensional, isotropic smoothing spline over the spatial coordinates was included in a Generalized Additive Model framework (Wood 2017). Further, no group effects for tracts and tract corners could be estimated because of strong confounding with the two-dimensional spatial trend function. Refitting was necessary because of minor changes in the data basis as well as in the predictor set. First, information on site soilwater regime was updated. Further, we detected a considerably high correlation between two different topex (1) (3) sums spanning a wide angle around the main wind direction. Consequently, we replaced them with one topex sum. In a second, multilevel approach we accounted for the inventory design by which tree individual i is nested in tract j and tract corner k. We therefore included group intercepts for tracts, and for tract corners: Here, we tested two different strategies for replacing the spatial smoothing function of SSE by the predictor . In a parametric model, PGS, we assumed an exponential effect of on the response. In this scenario we added 2 to x . By fixing the exponent to the value 2 here, we followed the rationale of fluid mechanics where the drag force affecting a given area (i.e. tree crown) is proportional to the squared velocity (i.e. gust speed) (Moore et al. 2018).
The second approach for replacing the spatial smoothing function is to include the effect of as a smoothing spline within a Generalized Additive Mixed model framework (SGS). Here, the effect of on tree damage is modelled by a thin-plate regression spline (Wood 2017). This flexible framework allows for the selection of the most plausible functional relationship between and storm damage probability, based on information from the data only, and not by any domain expertise as utilized in the first strategy. In order to elucidate the importance of as a predictor, PGS and SGS were compared against a reference model (NGS) which included all predictors but not .
All model formulation and fitting was conducted within a Bayesian framework. Continuous covariates were centered by their arithmetic mean and scaled to the unit of one empirical standard deviation. We used Normal(0, = 5) as a weakly informative prior distribution (Lemoine 2019) for every population-level parameter in . The estimation of variances 2 u , and 2 v for group level parameters was based on a half Student-t prior with 3 degrees of freedom and a scale of 5. Models were estimated by launching four parallel Markov chains of 2000 iterations from which the first 1000 steps were discarded as warm-up. After model fitting chain convergence was assured by checking the Gelman and Rubin convergence diagnostic R to be less than 1.01 and the effective sample size being greater than 1000 (Gelman et al. 2013).
Predicting with multi-level models for data not belonging to the National Forest Inventory of Baden-Württemberg faces the problem that group membership (in our case tract and tract corner groups) is unknown. One possible option, marginal predictions with group-level effects set to 0, was shown to lead to poorer predictive performance (Skrondal and Rabe-Hesketh 2009;Pavlou et al. 2015). Therefore, we fitted a mixture distribution model with the combination of estimated tract and tract corner effects on a logit scale. We drew 1000 samples from this mixture distribution, added each to the linear predictor combination, and calculated damage probability by inverting the logit scale and taking the mean. This resulted in a marginal prediction informed by group level effects. A similar problem arises when predicting with SSE: the spatial smoothing function is parameterised with coordinates of Baden-Württemberg. When predicting storm damage probability outside this spatial domain for South Lower Saxony, we projected averaged conditions of storm event Lothar onto the model estimates. Therfore, we sampled the spatial effects of 1000 tract coordinates of the parameterisation data set. Each of those were added to the linear predictor combination and transformed to damage probability by the inverse logit function. Finally, the mean of the resulting 1000 probability values was used.

Validation
Sensitivity analyses of predictor effects were performed by using a time series of dendrometric variables estimates for a simulated forest stand at an examplary location at the western margin of the Harz mountains, Germany (Fig. 1C). The simulated forest was represented by a mean stand tree for each species, of which was calculated using a climate-sensitive site index model (Schmidt 2020) and was predicted by a site-sensitive inverted longitudinal height-diameter curve (Schmidt 2009). This development of and over time was calculated for ages 30-130 (30-150 for Quercus sp. + F. sylvatica) in time steps of 10 years. The year of germination of these mean stand trees was set to 1950. Beginning with 2020, climate input data was retrieved from a subset of EURO-CORDEX ensemble containing seven representative projections (ReKliEs-De; Warrach-Sagi et al. (2018)) for scenario RCP 8.5. For each time step and age respectively, the corresponding and values were averaged over these projections. Predictive performance was tested by two approaches. First, all models were compared by a 20-fold cross validation with the parameterisation data set of Lothar. Data were split by keeping all observations of the same tract together. This ensured independence between training and test data. Two measures for predictive performance were used for validation: expected log pointwise predictive density (elpd) (Vehtari et al. 2017) and area under receiver-operator characteristics curve (AUC ). Elpd is a measure for out-of sample predictive accuracy, which uses within-sample fits with higher values for better performance. AUC is a common classification metric which relates the true positive rate to the false positive rate for binary predictions at different probability thresholds. For model comparison, AUC Lothar was averaged over the 20 cross validation runs and AUC Friederike was calculated with the independent validation data for Friederike in Lower Saxony.
All data analyses and plotting were done in the statistical programming environment R (R Core Team 2021). Fitting and cross-validation of models was realized with brms (Bürkner 2017(Bürkner , 2018, which is based on Stan (Carpenter et al. 2017).

Results
Model fits showed significant effects for most predictors, 95% credible intervals did not include a value of 0 (Table 1). Exceptions were the interaction between other deciduous species with as well as two : groundwater-influenced and terrestrial. Storm damage models SSE, PGS, and SGS showed same effect directions for covariates they shared, with smaller amplitudes for SSE than for PGS and SGS, especially for the dendrometric variables. In general, covariate effects of PGS and SGS differed only slightly. Figure 2 A allows for effect comparison of the two different approaches to model . The non-linear effect of SGS showed the strongest monotonic increasing effect in a range from 31 to 46 m s −1 where it increased by 5 units on the logit scale. At values above 46 m s −1 effect variation was very high and its average was decreasing with higher speeds due to lower data density. PGS constrained the effect of to a quadratic relationship. The effect increased monotonically and in total it spanned a large amplitude of more than 10 units on the logit scale. This large amplitude resulted in a steep increase of corresponding damage probabilities (Fig. 2  B). At top between 50 m s −1 and 60 m s −1 , damage probabilities predicted by PGS increased from 0.2 to more than 0.85. This is in contrast to estimates from SGS, whose upper bound of the 0.95 credible interval was less than half of the level of PGS. In the same range, the posterior mean of SGS was with less than 0.15 very low and even decreased with higher values.
Combining the effects of tree and its interactions with tree and in PGS, P. abies showed the highest absolut storm damage probabilities (Fig. 3), followed by the species groups A. alba + P. menziesii, P. sylvestris + L. decidua and Quercus sp. + F. sylvatica. With = 30 m s −1 , all species have similar damage probabilities at = 15 m . With tree = 30 m , P. abies showed a relative damage probability 1.5 times Table 1 Coefficient estimates of different model variants (columns), that were: a spatially structured effect (SSE), a parametric representation including squared gust speed (PGS), a spline function of gust speed (SGS), and excluding the parameter gust speed (NGS) Values are means of the posterior distribution and their 95% credible interval (CI). If CI includes 0 figures are set in italics. Estimates for smoothing splines are given as variance components with 0 for a straight line and the higher its values the higher the degree of wiggliness. Degree of freedom (df) defines the number of estimated parameters     higher than A. alba + P. menziesii, 2 times higher than P. sylvestris + L. decidua, and 3 times higher than Quercus sp. + F. sylvatica. With increasing , P. abies still showed the highest damage probabilities, but the relative difference to other species decreased. For instance, at = 30 m the damage probability of P. abies was 3 times higher than F. sylvatica with = 30 m s −1 , but 2 times higher with = 50 m . Slight differences within tree species groups were caused by species-specific estimations of the dendrometric variables and hence, height-to-diameter ratios. This difference was most pronounced for Quercus sp. + F. sylvatica over all heights, but also P. menziessii + P. abies differed slighly at > 35 m and P. sylvatica + L. decidua at > 30 m Storm damage showed the highest sensitivites for topex sums in and directions (Fig. 4). Negative topex sums indicated storm-exposed sites with the highest damage probabilities. Sheltered sites had positive topex sums and the lowest vulnerabilities. Negative topex sums described topographic ridges stretching in the wind direction and were associated with low storm damage risk. Positive topex sums indicated locations in narrow valley bottoms stretching in the wind direction. Here, high storm damage probabilities were associated, as a result of a tube-effect, with higher expected wind speed. Among all three topex sums, the variant showed the weakest effect amplitude. Pooling of soil water mapping units resulted in four storm-specific categories: terrestrial, shallow terrestrial, groundwater-influenced, and waterlogged soils. Comparing the sensitivity to these categories, periodical waterlogged soils and shallow terrestrial soils showed higher damage probabilities than terrestrial and groundwater-influenced soils (Fig. 5). Low probabilities for terrestrial and groundwater-influenced were a result of effect sizes close to zero (Table 1).
Both group effects on tract and on tract corner level, as well as their combination, showed a bimodal distribution pattern ( Fig. 6A-C): a slightly negative mode with low variation and a distinctly positive mode with large variation. This empirical distribution was fitted by a mixture model of two distributions (see supplementary information SI1): a Normal(−1.4, 0.9) and a Student's t distribution with = 77 (degrees of freedom), = 3.23 (location), and = 3.57 (shape). The mixture weight for the normal distribution was 0.7. The empirical density distribution of combined group level estimates could be approximated by 1000 samples of the mixture distribution (Fig. 6C).
Out-of-sample 20-fold cross validation using elpd showed that model SSE performed best (Table 2). By including as a predictor, PGS and SGS showed better results than model NGS. A similar picture provided the ranking of models by AUC values. In contrast, predictive performance with independent data from storm event Friederike did not differ between models PGS and SGS but showed a higher AUC value for NGS. Again, SSE showed best performance here.

Discussion
The multi-level modelling approach presented here integrates a) site-related factors like orography and soil features, b) species-specific dendrometric characteristics, and c) near-surface gust speed information to estimate a single tree storm damage probability. This probability provides a comprehensive assessment of forest vulnerability to winter storm damage. By replacing the spatially structured effect and PGS included squared as a linear covariable. The effect is displayed A on the scale of the predictor term and B transformed to storm damage probability. Black lines and different shading intensities show median and quantiles of the posterior distribution. All other metric predictors were set to their mean. Categorial predictors were set to tree species Picea abies and soil water category terrestrial. Ticks on the x-axis mark values of the parameterisation data. (Colour figure online) of SSE with in PGS and SGS, we included direct spatial information on the driver of tree damage in the model. By modelling this direct agent for storm damage, we aimed for a high degree of generalisation in the model. However, in our validation results SSE outperformed both models with as a predictor. The spatially structured effect of SSE accounts for autocorrelated heterogeneity on a regional scale (Brezger and Lang 2006). as a direct driver for storm damage most  Fig. 4 Sensitivity of storm damage probability on three different topex-to-distance sum effects from model PGS on storm damage probability over different tree heights. Effects were calculated for the species P. abies with the same and estimates as in Fig. 3.
The remaining topex-to-distance sums were set to 0, was assumed to be 40 m s −1 and was set to terrestrial. (Colour figure online) additional impact. The spatially structured effect, in turn, could interfere with spatially dependent covariates, which is known as spatial confounding (Clayton et al. 1993) and leads to biased effect estimates (Dupont et al. 2022). Indeed, we found less accentuated effect sizes in SSE. In contrast, the spatial predictor and the hierarchical group effects of the approaches PGS and SGS clearly separate between wind-driven damage patterns and unmeasured local effects.
Comparing the different methods to implement , a datadriven smoothing spline and a theory-based assumption of a quadratic relationship, the smoothing spline of SGS showed a diminishing effect at very high . At exposed sites trees adapt to frequently windy conditions, for instance by thigmomorphogenesis (Bonnesoeur et al. 2016) or by a closer height-to-diameter ratio. The latter was modelled explictly by the approaches presented here, but it is likely that the effect was confounded by a higher occurrence probability of trees with close height-to-diameter ratio at sites with top . Further, it is possible that damaged trees are attributed with higher than the actual value that lead to tree damage. The temporal aggregation of to a maximum value over the passage of a storm event does not resolve the total duration, the frequency, or the sequence of gusts affecting the trees. For example, less extreme but repetitive could also lead to swaying and subsequent stem-breakage (Brüchert and Gardiner 2006;Jackson et al. 2019), especially when trees are not adapted to such dynamics (e.g. inner stand individuals ). In contrast, PGS constrains the effect of to an exponential increment, particularly at high above 50 m s −1 . For the predictions it is more plausible to project an increasing effect with higher onto damage projections rather than dampening the effect because of possible tree adaptations. Although the effect in PGS was superimposed by the quadratic relationship and, hence, was not as close to the data as compared to the spline of SGS, both models showed similar validation performance.
Even though is the physical driver for storm damage, it did not improve prediction performance of PGS or SGS when validated with data from Friederike. Albrecht et al. (2019) reported similar findings. They modelled forest damages after two different storm events in Baden-Württemberg, Germany, and found a significant contribution of for only one event Lothar, but not for the earlier winter storm Vivian/ Wiebke in 1990. Conditions such as soil frost or high soil water saturation can vary substantially between winter storm events, which in turn has an effect on the importance of as a damage predictor. In the case of event Friederike, the storm had been preceded by a long period of rainfalls which resulted in high soil water saturation and a lowered root-soil anchorage. Here, tree , , and would probably better explain damage occurrence than nor . Both multi-level models showed high sensitivities of damage probability for covariates related to tree dimension and slenderness in interaction with tree species. The ranking of species according to their storm damage susceptibility is consistent with the findings of the predecessor study by Schmidt et al. (2010) but also with several other studies as reviewed in Gardiner (2021). Especially P. abies and A. alba  of conifers form a larger surface area than those of deciduous species in the winter season. Additionally, the shallow root architecture of P. abies promotes a higher susceptibility to windthrow. Site-related effects of soil water categories and topex sums quantify the abiotic predisposition of a location. Positive coefficients of the soil categories shallow terrestrial and waterlogged underline the importance of limited anchoring due to shallow rooting zones. Although this should also apply to groundwater-influenced soils, the corresponding effect did not influence the damage probability. These soils are less common in low mountain range regions like Baden-Württemberg and they are typically restricted to less exposed valley floors. This might have had lead to confounding with e.g. topex . Additionally, it is possible that the effect differs substantially between tree species because of different root architectures. However, the data basis did not suffice to model tree species-specific effects of . The combination of different topex sums relative to the main storm direction captures a large variety of orographic settings ranging from sheltered or exposed slopes to more complex forms, such as ridges or valley bottoms with different alignments. By step wise model selection, Schmidt et al. (2010) showed that these wind-direction sensitive topex sums better assess the orographic vulnerability to storms than the original topex, which sums up all eight (inter)cardinal directions (Scott and Mitchell 2005). Nevertheless, this implies for the assessment of possible damage vulnerabilities, that assumptions of the main storm directions have to be made. Most winter storms in Central Europe originate in the North Atlantic Oscillation and were driven by the westerlies, so that (south)west-orientated topex values are likely to depict the orographic exposure best.
Multi-level group effects of tracts and tract corners showed high variability compared to the effect sizes of the predictors. This variability integrates all regional and standlevel factors which have not been represented by the other predictors. The bimodal distribution of group level effects included a cluster of positive values. Here, more damage on tract or tract corner level occurred than was explained by the predictor set. Besides missing important measurable covariates, there are other causes which impact storm damage like stochasticity and damage propagation which are hard to quantify. In a storm event wind turbulence and resulting tree motion in combination with variation of individual tree features, such as rooting resistance, crown shape, or stem mechanics, lead to single tree failure and, due to damage propagation, to the loss of whole tree groups or stands (Dupont et al. 2015;Kamimura et al. 2019;Gardiner 2021). While simulations with mechanistic models (Dupont et al. 2015) and with agent based models  successfully emulated patterns of damage propagation, statistical models rarely used information from the stochastic  (2012) also found good performance for non-linear, multilevel taper models when using mean predictions over the distribution of group level effects. For predicting storm damage probabilities with stochastic information of the group level effects, the assumption is that the bimodal distribution found is rather a general winter storm related pattern than a event Lothar specific feature. However, climatic and environmental conditions, as well as the duration and intensity, differ between storm events and impact on the frequency and extent of tree damages. This inter-storm heterogeneity cannot be properly adressed because extensive data sets on damages caused by different events in a region of interest are missing. Nevertheless, the Lothar data set is extraordinary in terms of spatial coverage, number of observations, and covariate gradients covered, so it is most likely that the models depict main aspects of winter storm dynamics in Central Europe.

Conclusions
Winter storm damage probability was modelled in a comprehensive statistical framework based on a large data basis. This was accomplished by combining a deterministic part, including the most relevant features of dendrometry, soil water regime, orography, and a detailed spatial representation of gust speed, with a stochastic component describing unobserved variance at group level. Uncertainties about the origin of spatially correlated effects of its predecessor model SSE are replaced by the clear separation between information from the direct damage agent and stochastic variability in group effects. Due to the lower complexity and parametric formulation, the foundation in physical law, and the plausibility of the effect, PGS is most suitable for assessing storm damage vulnerability in forests. This framework is highly generalisable and allows for the assessment of winter storm vulnerability, for example, in the context of strategic forest planning. Further, the climate-sensitive predictors tree , , and , allow for the projection of vulnerabilites of actual or planned stands to future conditions (Jung and Schindler 2021). Within a forest enterprise, this information can be used to minimise storm damage risk by scheduling time of yield for present stands and optimising the spatial distribution of tree species for future stands.