Effects of environmental factors and forest management on landscape-scale forest storm damage in Turkey

Abstract

Key message

Windstorms have recently caused noteworthy destruction across the Northern forests of Turkey. The intensive forest management practices applied for more than 60+ years may unknowingly have resulted in wind-sensitive forests in the region. After a storm, the forest service salvages the losses, but no further precaution is taken against future storms. To our knowledge, there has not been any research looking into the cumulative effects of environmental factors on storm damage in Turkish forests. Maxent, which is an ecological niche model, might help decision-makers in developing forest management strategies against storms given its ease of use, known successful performance, and flexible variable evaluation approach. This study revealed that management preferences were mainly responsible for forest storm damage in Kastamonu Province, Turkey.

Context

Excessive wind cause serious damages to individual trees and forest stands. When unintentionally coupled with the forest management preferences, catastrophic levels of damage might be unavoidable.

Aims

The main objective was to assess the environmental factors contributing to the impact of a strong windstorm that occurred between March 14 and 15, 2013 and resulted in 1.5-million m³ timber losses in the Kastamonu Regional Directorate of Forestry.

Methods

Maximum entropy modeling (Maxent) and geographic information systems (GIS) were used to evaluate the factors contributing to the forest damage.

Results

Stand type, diameter class, and elevation were the most important variables affecting the level of wind damage. The pure and mixed coniferous stands were the hardest hit when compared with the deciduous stands. The damage increased as the density of forest roads grew.

Conclusion

It was concluded that windstorms pose serious threats to Turkish forests. Storm damage risks must therefore be integrated into forest management. In order to better understand the environmental factors contributing to the destructive effects of windstorms in forests, it would be best to focus on the telltale signs pointing the wrong-doing in forest management preferences at larger environmental scale rather than looking for reasons behind the occurrences of scattered small-scale damage.

1 Introduction

The most important principle in forestry is now accepted since the beginning of 1990s as forest sustainability (Wiersum 1995). Forestry aims to provide optimal economical gain while ensuring the continuous provision of ecological, social, and cultural services (Bruenig 2016). There are biotic, abiotic, and anthropogenic factors threatening this sustainability. Windstorms are regarded as one of the leading factors in forest damage (Reyer et al. 2017; Gregow et al. 2017). Mountain forests dominated by coniferous tree species were found susceptible to extreme wind speeds (Jung et al. 2016). Maximum gust speed was found to be the most critical factor for forest destruction (Schindler et al. 2016). In total, 110 million (mil) m³ of timber were lost to storms in Sweden in the twentieth century (Nilsson et al. 2004). However, studies showed that the damage could be reduced by 50% by decisive forest management practices in these conifer forests (Albrecht et al. 2015). Partial harvesting creating small openings of < 1 ha in British Columbia, Canada, did not reduce the intensity of storm damage in western hemlock (Tsuga heterophylla Raf. Sarg.) and red-cedar (Thuja plicata Donn. Ex D. Don in Lamb) forests, but it increased the susceptibility of other coniferous and deciduous tree species (Coates et al. 2018). In another study, an uneven-aged forest structure was found more resistant than an even-aged structure against a major winter storm in Switzerland (Hanewinkel et al. 2014). The common denominator in the majority of such studies is access to long-term experimental forest growth and yield monitoring data (Albrecht et al. 2012a; Albrecht et al. 2012b; Brazdil et al. 2018). Schütz et al. (2006) examined the interactions between maximum wind speeds captured with Doppler radar and topography, tree attributes, past silvicultural treatments and neighborhood characteristics. In another study, Donis et al. (2018) painstakingly sampled 169 1-km² study plots to reveal the windstorm damage triggering factors at the stand level. If accompanied by snow, winter storms tend to break coniferous trees (Altunel 2019), as opposed to uprooting them, particularly Norway spruce (Picea abies (L.) Karst.) (Peltola et al. 1999). However, heterogeneity in stand structure helps them better withstand storms (Diaz-Yanez et al. 2017). Introducing deciduous trees to predominantly monoculture coniferous forests would raise the forests’ resistance to storm damage (Jactel et al. 2017; Valinger and Fridman 2011). However, the level of such mixing needs to be understood; otherwise, the expected outcome in a future storm might not be achieved (Priebe et al. 2018). As trees get older, they also usually become more susceptible to such disturbances. Therefore, adaptive management strategies such as conversion of pure and even-aged stands into mixed and uneven-aged stands, shortening the rotation cycles, and tending frequencies should be developed to increase the resilience of stands against the destructive effects of windstorms (Jonsson et al. 2015). The stand characteristics, preliminary thinnings, soil type, and topography were found significant in all types of storm damage that occurred in Finnish and Estonian forests (Ilisson et al. 2005; Suvanto et al. 2016). However, in the absence of dependable long-term datasets, technological innovations, time, and man-power, niche models might offer insight into the factors behind wind damage within the context of ecological modeling (Muller-Landau et al. 2006), due to their repeatability, cheapness, and easy applicability. The interest in such a model, maximum entropy modeling (Maxent) (Phillips et al. 2006; Phillips et al. 2017), has lately increased due to an abundance of its applicable fields ranging from ecology (Mi et al. 2017), to the probability of forest fires started by lightning (Chen et al. 2015), earth science (Rahmati et al. 2016), and forest inventory (Mateo et al. 2018). However, to our knowledge, there has not been any research that has determined the model’s suitability in determining environmental factors contributing to storm damage in forests.

Windstorms cause forest damage and their range in severity and extent affect the strategic forest management objectives. Current forest management preferences and climate change might aggravate the situation even further (Lindner et al. 2010); therefore, alternative methods should be devised to improve the resilience of forests. Although some unexpected smaller-scale winds have always created occasional damage, larger windstorms hit Kastamonu forests in 1993, 1998, and 2001–2002 causing 35,000-m³, 40,000-m³, and 90,000-m³ timber losses, respectively. The most destructive and comprehensively recorded windstorm, causing 1.5-million m³ timber destruction, occurred within Kastamonu Regional Directorate of Forestry (KRDF) between March 14 and 15, 2013.

We assumed that other than the wind and its derivatives, i. e., duration, direction, gust speed, etc., there must be other environmental and management-related factors affecting the level of damage in the forests. Therefore, the main objective of this study was to investigate which factors influence the results of severe windstorms on forests. Moreover, the specific aim was to compose a windstorm susceptibility map for Kastamonu Province so that forest management prescriptions could be better developed.

2 Material and methods

2.1 Study area

This study was conducted on a 12,861-km² domain, in the Western Black Sea Region of Turkey, within the KRDF (Fig. 1).

Fig. 1

Location of the study area

The KRDF includes two neighboring provinces (Kastamonu and Sinop) and is ranked as the most productive region in Turkish forestry with a standing volume (SV) of 195 million m³ and an annual increment (AI) of 5.7 million m³ (URL-1 2018). It should be noted that the province of Sinop was excluded from the analyses, because some of the variables used for Kastamonu were not available for Sinop. Constituting the larger part of the KRDF, the Kastamonu Province has 67% of its area covered with forests. It has a SV of 142 million m³ and an AI of 4.3 million m³. Besides, Kastamonu produced 1.4 million m³ (68%) of the regional directorate’s total annual timber cut of 2.05 million m³ in 2017. The elevation ranges between the Black Sea coast and the summit of Ilgaz Mountain chain at 2578 m. The majority of the province’s vegetation is covered with the forests ranging from red pine (Pinus brutia) in the south-facing slopes of the low to medium elevations to pure and mixed stands of Corsican pine (Pinus nigra), Scots pine (Pinus sylvestris), and fir (Abies nordmanniana subsp. equi-trojani), as well as a wide range of deciduous, coniferous, and herbaceous plant species in the medium to high elevations (Gunal 2013). Two mountain chains, Kure and Ilgaz Mountains, running parallel to the coast in the north and in the south, respectively, create numerous watersheds varying in size and with varying climates. The northern chain cuts the warmer sea weather from reaching the inner parts of the province, thus creating harsher climatic conditions in the interior. The northern slopes of the Ilgaz Mountain chain receive a higher amount of precipitation, creating a more humid environment and fluctuating daily temperatures. When such temperature differences are present, these can create strong thermal winds (De Wekker et al. 2018). The average annual temperature and precipitation within the region are 10 °C and 450 mm, respectively. Prevailing winds come from the south-westerly directions approximately 40% of the year. The average wind speed is 2.2 m s⁻¹ with occasional gusts up to 30 m s⁻¹; these figures were the culmination of an almost 55 years of data collection from 16 meteorological stations from across the province (Ucar and Balo 2010).

2.2 Methodology

Maxent is a machine learning algorithm that estimates the probability distribution (for wind damage in this case) with the maximum entropy and closest to a uniform distribution (Phillips et al. 2006). It evaluates each given location, already defined in a rasterized environment, with respect to the questioned phenomenon while considering the variables to be analyzed. It, then, seeks the same string of variable combinations within the questioned domain. The maximum entropy probability distribution corresponds to the model’s probability surface varying between 0 and 1 (Erturk 2017). It has been favored by researchers simply because it provided a rather successful model prediction, irrespective of the sampling size (Hernandez et al. 2006). Although this is the first study utilizing Maxent for assessing forest wind damage, other studies, using neural network approach (Hanewinkel et al. 2004) and machine learning techniques (Hart et al. 2019), have been used to assess the factors controlling wind damage in forests. The model was built over a raster cell structure, so that it was expected to estimate the likelihood of the phenomenon in each cell as a function of the environmental variables (Royle et al. 2012), which were made available within the specified domain at the same sampling resolution. Maxent was favored for being less sensitive to the number of input variables and insensitive to the environmentally observed errors resulting from the variation in spatial location (Wisz et al. 2008). However, studies have shown the model accuracy and balance to be higher when the sampling size increased (Chen and Lei 2012). A sampling size of 1 hectare (ha) was used in this study.

2.2.1 Modeling

At the compartment level, the damaged area(s) were known within 2 months following the storm because all sub-departments of KRDF, e.g., the forest enterprises and their field directorates, completed a damage inventories. The damage was compartment-wise classified as “discrete” (sporadic and individual), “group” (sporadic but in groups of two, three, or more patches), and “cluster” (widespread in blocks). In addition, all field personnel were instructed to make a sketch of the compartment displaying what type of classification was assigned and where within the compartments they occurred. Location was determined by global positioning devices (GPS) creating a spatial distribution of the damage. These digitally stored data were used to make a completely new polygon fashioned layer of the study area. This layer can be considered as the fundamental building block of this study and was classified into two distinct classes, “damaged” and “undamaged” (Suvanto et al. 2016). This approach eliminated the need to set up a threshold defining what part of the forest was damaged and what was not.

This new layer constituted the basis from which the sampling scenarios were constructed to evaluate the sensitivity of any given location in the study area and to assess the environmental factors contributing to damage. First, we converted the damage from irregularly formed polygons to a specified sampling size of 1 hectare (100 × 100 m) with each designated by a point with a centroid (x, y). Thus, a point cloud was revealed as the sampling dataset averaging the terrain, management, and other related features in a compact, better manageable, 1-hectare scale. This procedure gave n = 2203 point locations from all the damage-assigned polygons (called single-point, presence-only type modeling) (Fig. 2(a)). However, the damage sampling was only from around the central location of each damage-assigned polygon varying in size tremendously according to the nature of the damage at that location. This was done because the model, Maxent, requires a set of such points as sampling input. Due to the random and sporadic nature of the damage, a spatially biased model was created with a bias toward the worst affected regions of the province. A fundamental limitation of presence-only type data is the sample selection bias (where some areas in the landscape are sampled more intensively than others) and this creates a much stronger effect in presence-only models than that in presence-absence models (Elith et al. 2011).

Fig. 2

The distribution of the damaged compartments (single vs. multi-point sampling)

Later, to offset the effect of sampling bias, and to avoid representing oddly shaped damage-assigned polygons with a single spatial location, we tried a second scenario by intersecting the entire study area with a 1-hectare grid. This was done by placing a 100-m × 100-m mesh over the study domain in ArcGIS 10.6.

In both sampling scenarios, the grid(s) simply displayed the damage locations. Beyond the positional accuracy, no further characteristic was assigned to them because the study was trying to associate the locations of these grid(s) with environmental variables.

The damage-assigned polygons depending upon their sizes included a number of hectare representing points. These still represented a dimension of 1 hectare but we believed they would better represent the damaged compartments because the damage area would be more systematically sampled. The points falling inside the damage-assigned polygons increased to n = 27,379 (called multi-point, presence-absence type modeling) along with the rest of the points falling within the province’s remaining forest areas (no damage-assigned points, n = 791,815) (Fig. 2(b)), which yielded enough points to run 10 iterations to strengthen the model performance during analyses.

The procedure transformed the polygon-shaped damage distribution and the remaining undamaged areas to a series of 1-ha points representing multi-point sampling locations and it helped to incorporate more sampling variability into the model and turned the random nature of the damage distribution into a systematic distribution (Kramer-Schadt et al. 2013; Yackulic et al. 2013). It is important for better model performance to include absence data (no damage) into the model (Radosavljevic and Anderson 2014). Trying to represent a damage-assigned polygon with a single 1-ha cell as tried in the single-point scenario caused less sampling variation in large damage areas and more representation in small damage areas. Thus, the multi-point scenario was devised to balance the uncertainties.

In order to capture those details, the sampling size might further be reduced to even smaller cell sizes; however, the practicality of the model might become questionable given the size of the region in this study.

Finally, storm damage–assigned points were input as sampling (presence) and the remaining points as absence into the model. In either setting, both the sample points and the environmental variables were prepared using UTM (Universal Transverse Mercator) projection and WGS-84 (World Geodetic System, 1984) datum at 1-ha spatial grid resolution. Maxent version 3.4.1 was used with the logistic option because the finalized damage sampling was in binary form (damaged/undamaged) and there were numerous categorical variables to deal with (Liu et al. 2018). ArcGIS 10.6 was used in the data preparation stages of all the variables and during the intermediate analyses.

2.2.2 Environmental variable acquisition

Province-wide fifteen distinct environmental variables were generated utilizing data from four different sources (Table 1). To coincide with the above sampling, the same “1-ha” grid structuring was again used to rasterize all the environmental variables. Thus, both in single and multi-point sampling, all the environmental variables aligned with the 1-ha sampling grid(s).

Table 1 The investigated environmental variables

Shuttle Radar Topography Mission (SRTM) 3-arc second (90 m) elevation data (DEM) were acquired from Consortium for Spatial Information data portal (CGIAR-CSI 2019). The study area was encompassed within six DEM tiles. They were downloaded and mosaicked to form the DEM of the study area. Using the finalized DEM, first “elevation,” second “slope,” third “aspect,” and finally “ruggedness” variables were calculated. The elevation ranged from 0 to 2578 m and percent slope ranged from 0% (flat) to 406%, while aspect varied between 0 and 359°. Terrain ruggedness was calculated as a derivative of the slope. Rasterized slope data are made up with a number of adjoining cells, each assigned with a slope value. Each cell within the data can be compared with its surrounding 8 or 24 cells. The variability among these cells is known as the ruggedness value for the particular cell. If all slope values within these cells matched, the cell in the middle got a value of 1. The more the variability among the investigated cells, the higher the terrain ruggedness (Riley et al. 1999). A nine-cell raster structuring was used in this study, and the data was formed as categorical (scores of 1 to 9).

Additionally, the regional forest road network was obtained from the vectorized forest road database of KRDF, converted to road density and rasterized as continuous data. It varied between 53 and 311 m/ha. Wind speed and precipitation were acquired from the state meteorology service in the form of daily field readings and were rasterized utilizing the inverse distance weighted technique (IDW) (Friedland et al. 2017), which was the distance-related influence of an input point across the interpolated surface and treated as continuous data. The wind speed varied between 1.8 and 21.3 m s⁻¹ on March 14, 2013. On the other hand, precipitation varied between 1 and 22.2 mm/m² on March 13, 2013. The previous day’s record of precipitation was used to see whether such a variable might have been a factor in damage probability. The variables for major soil types (MST) and soil depth (SD) were extracted from the vectorized provincial soil classification database. In contrast to the abovementioned continuous data, they were first classified, then rasterized, accordingly, to better differentiate the soil-related characteristics in the modeling. MST was classified in three categories as zonal [1], intra-zonal [2], and azonal [3]. SD was divided into 4 categories as rather shallow (0–20 cm [1]), shallow (20–50 cm [2]), less deep (50–90 cm [3]), and deep (> 90 cm [4]).

The remaining management-related variables, i. e., stand type, mixture form, forest type, diameter class, crown closure, and stratification, were extracted from the vectorized forest management database of KRDF, which included the records from the last inventory in 2012 (Table 2).

Table 2 Forest management specific variables

The majority of the environmental variables listed in this study were included because the literature supported our selection (Felton et al. 2016; Anyomi and Ruel 2015; White et al. 2015). A few of them, such as terrain ruggedness, road density, and stratification have not been tested previously in such a study, so we wanted to examine if they would make a difference. The results were further elaborated through jackknife analysis for comparing the environmental variable share in the model’s success and determining the variables exhibiting high correlation to one another. The process is that each variable is independently dropped from the process while a series of models is simultaneously executed (Phillips et al. 2006) to display the contribution of each variable on the overall success of the model (Phillips et al. 2017).

The model training and control datasets were adjusted as 75% and 25% of “sampling points in either modelling scheme” (Erturk 2017), respectively. In total, eight abiotic and seven anthropogenic environmental variables (6 continuous + 9 categorical-type data) were considered. Depending upon the model variable percent contribution and the result of jackknife statistics, the most marginal variables that were not registering any percent contribution to the model were omitted, and iterations were carried out until all the remaining variables contributed to the model. Variable response curves were additionally considered to fully understand the contributing variable’s contribution to the model. The variables “terrain ruggedness,” “soil type,” “mixture form,” and “stratification” were dropped from the analyses for not registering any contribution after 10 iterations in the first run. All the remaining variables were kept.

Before further analyzing the results, we checked the correlations between the variables, utilizing Pearson correlation coefficients (Weiss et al. 2016). A correlation matrix was generated through ArcGIS 10.6-Spatial Analyst-Band Collection Statistics. Crown closure to diameter class and forest type to diameter class pairs showed moderate correlation with r = 0.619 and r = 0.589, respectively. Although all three variables were the direct results of forest management preferences, forest type was therefore dropped from the analysis (Table 3). This further reduced the already weak contribution of the variable, slope, so it was also dropped from the analysis. The results were calculated with the remaining variables.

Table 3 Pearson correlation matrix

3 Results

The first model produced an AUC value of 0.824 with the single-point sampling. It was acceptably successful (Hosmer and Lemeshow 2000) for the investigated phenomenon despite the fact that the likelihood of an unintentionally induced operator bias was present. Whether the model sampling scheme represented by the points had been automatically generated or manually created inside the damaged compartments, there was a lack of within-compartment variability in the single-point sampling. The area of the damaged polygons ranged from 0.25 to 60 ha, which translated to a single hectare representing point in the smallest compartment and 61 such points in the biggest compartment. The number of damage-assigned polygons, which was bigger than 1 ha, was 2098, which meant that just 5% of the total number of such polygons (n = 2203) were suitably representable with a single 1-ha sampling. Thus, a considerable amount of within-compartment variability would have been lost. These reasons prompted the use of the 1-ha grid multi-point sampling (second method) instead. Although the AUC value of the resulting second model dropped to 0.739, we decided to adopt it as having sufficient success. Although there was a loss of 7% in discriminatory power, the stronger model contributions of a few leading variables convinced us to continue with the multi-point sampling approach.

Assessing the contribution of all the remaining environmental variables to the model was the second stage of model performance evaluation. The three most important and the three least important model-contributing variables were as follows: stand type (36.5%), diameter class (29.9%), and elevation (20.7%); and crown closure (0.3%), soil depth (0.5%), and aspect (0.3%), respectively (Table 4).

Table 4 Environmental variable contribution list

The corresponding percentages were the heuristic estimates of the relative environmental variable contributions to the model of forest storm damage probability. While determining the percentage, in each iteration of the training dataset, the additions in normalized gain is appended over the contribution of the corresponding variable, or taken off from it if the change is negative (Phillips et al. 2017).

When the model contributions of all the remaining variables were checked, it was apparent that they were distinct and could well be categorized as high, intermediate, and low level effects on the model. Stand types [1], [2], [3], and [5], which are found almost everywhere within the coniferous forest, and diameter classes [2], [3], [4], and [5] related to management preferences, were the hardest hit during the storm. Although lower in rank, but still highly effective in model contribution, elevation was another important environmental variable associated with the damage. Aside from a couple of instances in coastal locations, the bulk of the damage occurred within the forests located between 400 and 1100 m above sea level. Because these elevations are the prime forest habitat within KRDF, and the likelihood of a strong wind is always present, the need to adopt more forest storm damage mitigating strategies is clear.

Although lower in rank, but as expected, wind speed, above 4 m s⁻¹, had an impact on the level of damage. It might be the reason why “elevation” was one of leading model-contributing variables. Data provided by a better-situated meteorological station network and utilizing an airflow model like WAsP (Mortensen et al. 2003) could have increased the importance of wind speed.

Another lower ranking but important variable was the road density. There has not been an adequate study directly linking road density to forest destruction during the windstorms. However, the creation of linear edges, increasing forest fragmentation and gust speeds by channeling inside the forest areas, has been reported as problematic for forests (Dupont et al. 2018; Schwartz et al. 2017). The road pattern is related to logging and transportation resulting from the need for timber extraction. This has, in turn, culminated in many linear gaps and resulting edges inside the forests. Our results also showed that as the road density increased, the storm damage also increased, which was consistent with previous research (Torun 2018).

The jackknife diagram showed the average values of the environmental variables emerging after iterating the storm damage probability model 10 times (Fig. 3). The variables, such as stand type, diameter class, elevation, and road density, were the most influential, absence of which would cause discriminatory power loss in the model (blue bars represented the variable’s contribution in the overall model performance (red bar) vs. turquoise bars which represent the total model gain when the particular variable is excluded from the model). For example, if crown closure is excluded, the model’s gain did not change. Although making some contribution, the exclusion of the variables, soil depth, and, to some degree, precipitation would not make a large difference in the model’s performance. The effect of wind speed, by itself, was just below the most influential variables in the overall model performance. Aspect was the only variable, which did not contribute to the model.

Fig. 3

The jackknife diagram results of all the contributing environmental variables

Finally, to understand how the variables, which were the main building blocks behind the model’s logic, were influential on the resulting possibility, response curves were generated. They showed the logistic prediction shift as each contributing variable was modified while maintaining the remaining variables around their average sampling values. Thus, the curves manifested the effect of altering just one variable only, while the model evaluated the advantage of a number of other variables changing together (Phillips et al. 2017). The red or blue bars and lines were the average logistic prediction of each variable at the end of 10 iterations and the standard deviations occurring after each iteration, respectively (Fig. 4).

Fig. 4

Response curves of some model-contributing environmental variable(s)

The model output, which could be named as “the storm damage probability map” showed the forest storm damage probability of any 1-hectare spatial resolution within Kastamonu Province (Fig. 5). The color graduation from green to red shows the increasing wind damage probability during southerly windstorms and ranged between 0 and 1. We then reclassed the data into four easily discernable classes and named them as low probability (0 to 0.25), medium probability (0.2501 to 0.5), high probability (0.5001 to 0.75), and extreme probability (0.7501 to 1). When the area calculations were made over the map, 35.9% of Kastamonu forests were found to have high to extreme storm damage probability according to the data used in this study (Table 5).

Fig. 5

Forest storm damage probability map of Kastamonu Province, Turkey

Table 5 Windstorm probability acreages and percentages

4 Discussion and conclusion

Windstorms are the most important abiotic factor causing forest damage in many parts of the world. Turkey at this point, while preserving and steadily extending its forest area, has lately started to be hard-hit by such phenomenon. Despite the considerable damage experienced throughout the country, neither a comprehensive course of action nor a study looking into the underlying factors of the storm damage on forests has so far been conducted. A popular, but never previously used to asses storm damage, niche model Maxent was evaluated to assess the importance of fifteen different abiotic and anthropogenic variables, e.g., elevation, slope, aspect, terrain ruggedness, wind speed, precipitation, road density, stand type, mixture form, forest type, diameter class, crown closure, stratification, soil classification, and soil depth. Maxent has been the method of choice for habitat suitability in ecology and relies on well-defined single-point random sampling. However, it was impossible to describe a sporadic phenomenon like forest storm damage, applying such a sampling scheme. In single-point sampling, since the model only relied upon a number of random hectares representing points coming from the damaged areas, the setting created the near ideal data input for the model and the resulting AUC value scored high (0.824). However, the random nature of the damage and the resulting sampling transformation limited the model’s ability to work with the presence-only type modeling approach because there was no way of representing the undamaged areas with any type of methodologic sampling scenario.

In multi-point sampling on the other hand, since it generated a greater number of points from the damaged areas from the same number of compartments, more heterogeneity was incorporated into the model, and the resulting AUC value scored somewhat lower (0.739). Despite losing some discriminatory power, we opted to continue with the multi-point sampling. The smaller sampling size of a hectare (ha) compared with the almost industry-standard square kilometer (km²) in Maxent modeling (Liu et al. 2018; Remya et al. 2015; Kumar et al. 2014) meant inevitably that the model would lose some discriminatory power.

Five of the model-tested nine environmental variables were abiotic, which suggested that no matter what the foresters do, some storm damage was inevitable. However, two leading and critically important management-induced anthropogenic variables were among the damage defining variables. Stand type and diameter class were the two most influential environmental variables contributing to forest storm damage and are directly affected by forest management with stands categorized based on diameter classes and crown closure percentage in Turkey. The interventions when periodically applied to stands aims to gradually widen the distance among the neighboring trees, promoting the trees to grow and the crowns to extend at ease toward canopy closure. Although sustainably managed, the abovementioned interventions might sometimes unintentionally cause a tendency to create even-aged diameter distributions resulting in less diversity and vertical stand stratification, which can have negative effects on the stability of stands in the future (Drossler et al. 2015; Bebi et al. 2017). There is no way of knowing to what extent and when the last or preceding intervention(s) were applied. However, those interventions probably contributed to the storm damage to KRDF in 2013.

On the other hand, forest roads, which have long been used as the only means of taking all the management-related services to the forests by the Turkish Forest Service, have dramatically increased in distance and density in the region. They are an integral part of forest management (Lugo and Gucinski 2000); however, the more they have been studied, the more their indirect effects have been reported (Zeng et al. 2010; Ruck et al. 2012; Schindler et al. 2012). For example, the more forest roads that are built, the more linear edges are created in the forests, thus creating a lot of vulnerable forest edges (Zeng et al. 2010).

The KRDF forests are mostly natural although there are some examples of artificially afforested coastal plantations of “Pinus pinaster Ait.” and “Pinus brutia Ten.” The climatic conditions are ideal for tree growth so that even formerly cultivated but later abandoned areas are swiftly transformed by the adjacent forests. State forests and the stands are predominantly covered with Corsican pine, Scots pine, and fir. The heavily damaged stand types of [1], [2], and [3] (coniferous, coniferous + deciduous: largely coniferous, and coniferous mix) were all covered with these tree species with the exception of marginally damaged stands with broad-leaved species, beech (Fagus orientalis L.), chestnut (Castanea sativa Mill.), and ash (Fraxinus excelsior L.) in stand type [2]. The damage was less in stand type [5] (deciduous + coniferous mix: largely deciduous), considerably lower in [4] (pure deciduous) and almost none existent in [6] (deciduous mix). When the homogeneity is reduced in stands, their ability to withstand hazards is raised (Seidl et al. 2014). Since Turkish forestry’s main agenda has primarily been economic gain, the predominantly coniferous (Corsican, red, Scots pines; firs, spruce, cedar) species, although indigenous to the regions, are favored for their quick turnover and high yield. This approach, most of the time, leads to pure even-aged stand structures (Kara 2018) and higher susceptibility to windstorms. Thus, it would be wise to think that uneven-aged stand formation either pure or mixed compositions would be better able to withstand storm damage (Pukkala et al. 2016; Valinger and Fridman 2011). Developing forests with low risk from windstorms but with little commercial value is possible, but will not be of interest to forest managers with income targets. Therefore, a balance between the risk and the projected benefits might be possible by introducing some heterogeneity into the stands. Although disputed by Foster and Boose (1995) as equivocal, the model results showed that precipitation before the storm might possibly have saturated the already winter damp soil, further weakening the shallow root systems of predominantly coniferous trees (Schaetzl et al. 1989; Usbeck et al. 2010; Coutts et al. 1999).

As expected, the wind speed was among the variables found to have an impact. Although reported and validated as the leading abiotic factor in global forest destruction (Gromke and Ruck 2018), it only appeared among the intermediately effective variables in this study, because the distribution of the meteorological station network in and around Kastamonu was not ideal when the storm struck in 2013. Thus, the importance of both the wind speed and precipitation, and their real impact on the storm damage, may have been somewhat underestimated.

Three of the highly to intermediately important variables were anthropogenic and the other two were abiotic. It is very unlikely that we can alleviate the detrimental effects of the investigated abiotic variables in this study, given the present course of global climate change. However, there is a lot potential for changes in forest management and lessons to be learned from past mistakes. The storm threat must be integrated into forest management, and certain precautions, e.g., avoidance of uneven-aged monocultures, limiting forest road expansion projects, not excessively constraining the diameter classes just because industry needs wood, better monitoring of the climate, and past forest management interventions, must be made. The results have shown that the climate, by itself, is not to blame and improved forest management strategies must be adopted to ensure the sustainability of these precious resources in the long run.