Interpreting wind damage risk–how multifunctional forest management impacts standing timber at risk of wind felling

Landscape multifunctionality, a widely accepted challenge for boreal forests, aims to simultaneously provide timber, non-timber ecosystem services, and shelter for biodiversity. However, multifunctionality requires the use of novel forest management regimes optimally combined over the landscape, and an increased share of sets asides. It remains unclear how this combination will shape stand vulnerability to wind disturbances and exposed timber volume. We combined forest growth simulations and multi-objective optimization to create alternative landscape level forest management scenarios. Management choices were restricted to 1) rotation forestry, 2) continuous cover forestry, and 3) all regimes allowed over a harvest intensity gradient from completely set aside landscapes to maximal economic gain. Estimates for the stands’ structural and environmental characteristics were used to predict the stand level wind damage probability. We evaluated averaged wind-exposed standing timber volume and changing forest structure under management scenarios. Intensive rotation forestry reduced tree heights and wind damage risk, but also reduced landscape multifunctionality. Conversely, continuous cover forestry maintained multifunctionality but increased wind damage probability due to taller trees and higher thinning frequency. Overall, continuous cover forestry lowers the total volume of wind exposed timber at any given time compared with rotation forestry. Nevertheless, a selective application of rotation forestry contributes to high economic gains and increases landscape heterogeneity. A combination of management approaches across landscapes provides an efficient way to reduce the amount of wind-exposed timber volume while also increasing habitat for vertebrate and non-vertebrate species and satisfying high timber demands.


Introduction
Managing boreal forest landscapes for multifunctionality promotes the provision of a wide set of ecosystem services and increases forest biodiversity Hölting et al. 2020;Larsson Ekström et al. 2021). Compared with conventional, timber-oriented management, multifunctional management requires the development of novel management methods, and their diversification over the landscapes as well as inclusion of set asides for conservation purposes Triviño et al. 2017;Mönkkönen et al. 2018;Eyvindson et al. 2021). European boreal forests have been managed for centuries, with steadily increasing harvesting rates , negatively affecting forest biodiversity and deadwood volume (Vanha-Majamaa et al. 2007). To balance timber demands with non-timber forest values, a shift in management and planning is required. It remains poorly understood how such diversified management of forest landscapes will affect stand vulnerability to forest disturbances, reflecting economic risks in timber supply and quality in production forests.
Multifunctional management prioritizes a balance between conflicting objectives Eyvindson et al. 2021) and between different stakeholders preferences (Hölting et al. 2019) over the landscape. In Finland, the challenge to simultaneously provide several conflicting benefits is widely accepted, as forests provide many societally important goods and services Mönkkönen et al. 2018;FMAF 2019;Pohjanmies et al. 2021). This includes, e.g., recreation, non-woody ecosystem services, carbon storage, and suitable habitats for forest dwelling species . To balance between the conflicting societal demands and biodiversity, several initiatives for changes in forest management have recently emerged. One alternative is to replace traditional rotation forestry with continuous cover forestry, because of the overall negative effects of the former on biodiversity and ecosystem services (Pukkala 2016;Peura et al. 2018). Another alternative is to increase green tree retention at final harvest (Kuuluvainen et al. 2019;Gustafsson et al. 2020), or to extend rotation length to improve forest carbon stocks and features important for biodiversity ). Yet, a single management type cannot satisfy all of the conflicting objectives simultaneously at the stand level. Instead, landscape-level diversification of forest management is required ). According to Eyvindson et al. (2021), Finnish landscapes should be dominated by continuous cover forestry, combined with a relatively low share of intensive rotation forestry (10-25%) and complemented with set asides to maximize forest multifunctionality. Implementation of diversified management at the landscape scale will affect harvest scheduling, clear-cut sizes, differences in tree heights between neighboring actively managed stands and set asides, and frequency of thinning over the longterm. All of these changes have important implications for wind damage risk (Peltola et al. 1999b;Venäläinen et al. 2004;Heinonen et al. 2009;Zeng et al. 2010;Zubizarreta-Gerendiain et al. 2016).
Wind damage is one of the most important disturbance agents for forest management. Climate change will likely further amplify its impact due to the increasing frequency of strong winds (Nikulin et al. 2011), increasing growing stock (Seidl et al. 2011), and shortening of the frozen soil period resulting in lowered tree ground anchorage during the windiest time of the year Venäläinen et al. 2020). Several forest structural characteristics shape wind damage probability at the tree and stand level. At the single tree level, wind damage probability increases with tree size and increasing wind loads (Hale et al. 2012(Hale et al. , 2015. A mechanistic understanding about the wind load and individual tree resistance to wind allows for the prediction the critical wind speed at which an 'average tree' in a stand uproots or breaks (Peltola et al. 1999b;). The important tree-level factors are the tree's diameter, height, and species (Peltola et al. 1999a;Gardiner et al. 2000), where Norway spruce (Picea abies (L.) H. Karst) is the most wind-vulnerable tree species in Finland (Peltola et al. 1999b;Zeng et al. 2010). Locatelli et al. (2017) found that tree height, stocking density, stem diameter and the size of any upwind gap to be the key factors controlling risk in the model. At the stand level, wind damage probability increases with greater basal area and stand density (Díaz-Yáñez et al. 2019), with proximity to recently thinned stands or clear cuts (Zeng 2006), and a presence of upwind gap (Locatelli et al. 2017). Wind damage is often located near stand edges compared to the stand interior (Zubizarreta-Gerendiain et al. 2012). Therefore, to reduce wind damage risk over the landscape, forest management should aim to minimize the stand height heterogeneity between neighboring stands (Heinonen et al. 2009;Zubizarreta-Gerendiain et al. 2016), omit extensive clear-cuts (Zeng et al. 2007) and increase the structural diversity within the stand Jactel et al. 2017). In addition, forest vulnerability to wind depends on specific site conditions such as soil type, depth, and fertility affecting tree anchorage potential (Day 1950;Gregow 2013).
For stands under even-aged rotation forestry, wind damage risk increases due to the increased forest fragmentation following final cuts and heavy thinning (Zeng et al. 2007(Zeng et al. , 2010Zubizarreta-Gerendiain et al. 2012, 2016, and due to height heterogeneity between neighboring stands (Heinonen et al. 2009). On the other hand, shortening of the rotation length and harvesting younger trees can reduce overall wind damage probability (Kellomäki 2017), but conflicts with biodiversity (Roberge et al. 2016). Continuous cover forestry with the selective harvests of individual trees in an uneven-aged stand can reduce overall wind damage . However, continuous cover forestry will likely require higher thinning frequency per stand than in rotation forestry (minimal length of 15 years between two consecutive thinning, Äijälä et al. 2014), and recent thinnings generally increase the wind damage risk to remaining trees and in neighboring stands (Zeng et al. 2010;Pukkala et al. 2016;Duperat et al. 2020). This effect seems consistent across thinning approaches, but varies in magnitude depending on the stand management history. For example, Pukkala et al. (2016) found that thinning from above (selectively harvesting the tallest trees in the stand) in uneven-aged stands was less risky than thinning from below in even-aged stands. Overall, thinning from below (removing of the smaller trees and maintaining the canopy trees) had the highest increase in wind damage risk from all thinning types . In contrast, in balsam fir stands in Canada, thinning from above increased vulnerability to wind damage and reduced overall stand stability more than thinning from below (Duperat et al. 2020). Therefore, thinning generally appears to increase wind damage risk, but the magnitude of the effect remains unclear and depends on details and stand management type.
Finnish boreal forest stands are currently predominantly managed as even-aged stands with rotation lengths of 40-120 years (Valsta 2017, 60-80 years in Southern Finland, (Äijälä et al. 2014)), with only 6% managed as uneven-aged stands with continuous cover forestry management, and a small share of set asides depending on the forest certification scheme (Peltola et al. 2019). Therefore, the implementation of multifunctional management, dominated by continuous cover forestry, would affect forest structure and the harvest schedule, subsequently affecting the probability of wind damage and exposed timber volume.
To understand how multifunctional management would affect the wind damage risk and volume of wind-exposed timber, we have combined a wind damage probability model ) with a boreal forest growth and management simulator (Rasinmäki et al. 2009). We restricted the selection of management regimes to address a citizen initiative aiming to replace rotation forestry with continuous cover forestry in Finland in state-owned land (VN/1699/2018). The probability of wind damage and quantity of timber at risk was analyzed under three sets of management regimes: only rotation forestry and set asides (hereafter: RF scenario), only continuous cover forestry and set asides (hereafter: CCF scenario), or a combination of all management scenarios (hereafter: ALL scenario). Our research addresses three main questions: 1) How will multifunctional management over a range of harvest intensity gradients affect wind damage risk and the exposure of timber volume? 2) How will the restriction of management regimes affect the probability of wind damage and wind-exposed timber volume? 3) What are the main drivers of the differences in wind damage risk and wind-exposed timber volume in each management scenario?
We hypothesized that (1) increasing harvest intensity will increase wind damage risk due to higher harvest frequency (by thinnings and clear-cuts), and prevalence of Norway spruce as the most economically important species; (2) continuous cover forestry will lower wind damage risk due to the lower heights heterogeneity among neighboring stands, and will expose lower timber volumes to wind damage; (3) rotation forestry will have higher wind damage risk and higher wind-exposed timber volume compared with continuous cover forestry due to a creation of open stands edges. We expected these differences in wind damage risks between management regimes due to differences in frequency and type of thinnings, tree heights, and final harvest. The study aims to improve our understanding on how harvest intensification under the restriction of management regimes will shape wind damage risk and wind-exposed timber volume, as the proxy of potential economic losses.

Methods
Our study followed the workflow in Fig. 1: (i) collect current stand-level characteristics for input; (ii) simulate forest growth under alternative management regimes over 100 years; (iii) optimize landscape-level management to maximize multifunctionality for different levels of harvest intensity (net present income, NPI) under RF, CCF and ALL scenarios; and (iv) calculate the wind damage probability and volume of timber at risk for each stand under alternative scenarios.
Using multi-objective optimization, we produced sets of solutions reflecting the harvest intensity gradient-from no production objectives (completely set aside) to maximal economic gain while maximizing landscape multifunctionality. We evaluated multifunctionality as the combined normalized values for landscape recreational values, carbon storage, and habitats potential for vertebrate and non-vertebrate species . For each scenario, we predicted wind damage probability ) and recorded wind-exposed timber volume between actively managed stands and set asides. The wind damage probability model was developed from the Finnish national forest inventory data , where the dominant management has historically been rotation forestry. We further aimed to identify the dynamics of forest structural characteristics driving vulnerability to wind damage, such as dominant tree species, tree height, time since thinning, and height 1 3 heterogeneity between neighboring stands. We processed the datasets, calculated wind damage probability, and visualized results using R (R Development Core Team 2019) and QGIS (QGIS Development Team 2020).

Study area and input data
Our study area represents a typical Finnish production forest landscape under rotation forestry, with relatively uniform, homogenous forest structure, species composition, and a relatively young age distribution. In total, we used 1470 forest stands aggregated within a single watershed in Central Finland, covering 2242 ha (Fig. 2a). Initial stand conditions in 2016 were obtained as open-source data from the Finnish Forest Centre (available on www. metsa an. fi). The input dataset included stand-level characteristics such as the main species, tree height, basal area, soil type and fertility. The main tree species was pine (50.7% of stands), followed by Norway spruce (33.9%) and other tree species (not specified, Table 1). The dominant soil type was mineral (74.2%) and organic (25.7%). The majority of stands had mineral soil depth of more than 30 cm (99.7%,); 61.5% of stands were located on fertile and 38.5% on poor sites. The mean stand size was 1.5 ± 1.6 ha (mean ± sd). Time since thinning varied between stands, with the majority thinned more than 10 years ago (77%, National Forest Inventory, ).

Forest growth and management simulation
Forest stand development was simulated over 100 years, separated into 20 five-year sequences using the SIMO forest growth simulator (Rasinmäki et al. 2009). SIMO is an open-source tree-level forest simulator, incorporating over 500 models to forecast the development of the forest stand over time. The simulator forecasts tree growth, mortality and ingrowth over time, allowing for an estimate of both dead (including deadwood decay using models from Mäkinen et al. (2006)) and living trees. The simulation processes data at several levels, e.g., predicts the growth of individual trees, and aggregates the results into strata at the stand level. Stand can have one or several strata, where one stratum represents single tree species and size. All stands start simulation as uniform stands (originating from rotation forestry), and develop alternative forest structure over time, given an applied management regime. SIMO allows for the implementation of variations to forest management regimes such as time of final harvest, timing and the frequency of thinnings, and modification of final harvest type such as continues cover forestry. For this Fig. 1 Study workflow from collecting input dataset through forest growth simulation under the variations of rotation forestry (lower case: rf), continuous cover forestry (ccf) and set aside (sa). Management scenarios combine the alternative management regimes (upper case: RF: rotation forestry and set aside, CCF: continuous cover forestry and set aside, ALL: rotation, continuous cover forestry and set aside) with landscape-level optimization over constraints of harvest intensity, creating a gradient from minimal (100% set aside) to maximal net present income (NPI) and multifunctionality. Applied models originate from: forest growth simulator (Rasinmäki et al. 2009), forest management regimes , multi-objective optimization ) and the wind damage probability model  assessment we utilize 58 different management regimes that have been developed and tested in earlier studies . These represent 17 variations of rotation forestry, 40 variations of continuous cover forestry, and one set aside alternative, where no management actions were taken . Rotation forestry regimes differed in the timing of final felling, optional thinning (present/absent), and in a larger number of retained green trees after the final cut . The basic continuous cover forestry regime follows the guidelines from Äijälä et al. (2014). Continuous cover forestry alternatives cover the variations of the site-specific target basal area, and postponing the time of first thinning. For a detailed description and definitions of the rotation and continuous cover forestry management regimes see Eyvindson et al. (2018Eyvindson et al. ( , 2021.

Optimization over harvest intensity gradient
The focus of the multi-objective optimization is to maximize landscape level multifunctionality at a specific value of net present income (NPI) obtained from harvested timber volume. NPI represents the economic value of the forests, discounted to the present value using timber revenues minus silvicultural costs, discounted using a rate of 2%. We opted to use a 2% discount rate to only moderately prioritize near term incomes over distant incomes. This rate is often used in solving a long-term economic problem in forestry (e.g., Brukas et al. 2001;Heinonen et al. 2020). Using higher discount rates may dramatically impact the acquired NPI value, however the trade-off with multifunctionality will not significantly change .
Landscape level multifunctionality was calculated as a sum of four normalized components reflecting important Fig. 2 a Study area represents a forested landscape located in Central Finland comprising 1470 forest stands; b Schematics of the harvest intensity gradient from minimum (Landscape Xn) to maximum harvest intensity (Landscape Xn + m) over the landscape and under com-bination of optimal forest management scenarios. Actively managed stands are either under rotation forestry or continuous cover forestry. Set aside stands complement the landscape under each scenario with different shares across the harvest intensity gradient Table 1 Initial stand characteristics in 2016 for three main tree species: Scots pine (Pinus Silvestris L.), Norway spruce (Picea abies L.) and others (non specifies, mostly deciduous species), as specified in wind disturbance vulnerability model from  The values represent stands counts and their share (%), or mean values over stands ± standard deviation aspects for Finnish nature and people: i) recreational ecosystem services and non-woody services (set W: bilberry, mushrooms, scenic beauty), ii) set X: climate change mitigation, iii) set Y: suitable habitat for vertebrate-, and iv) set Z: non-vertebrate endangered species as surrogates for biodiversity conservation. Individual indices have been previously implemented in SIMO. Specifically, non-woody ecosystem services originate from Miina et al. (2016), for bilberry yield, Tahvanainen et al. (2016), for mushroom yield, and from Pukkala et al. (1995) for scenic beauty. Carbon storage considered the mass of carbon contained within timber, deadwood, and soil (Liski and Westman 1997). Habitat suitability indices for endangered vertebrates originate from Mönkkönen et al. 2014); non-vertebrate species depend on predicted deadwood type (Mäkinen et al. 2006) and habitat availability for endangered fungi and arthropods (Tikkanen et al. 2007). Landscape level multifunctionality follows Eqs. 1, 2a, 2b.
subject to: where D b represents the measured deviation for component set b, B is the set of components, T b is the set of indicators in component set b, f * t , f t * and f t , respectively, represent the ideal, anti-ideal and obtained value for indicator t. Four component sets were used, to reflect multifunctionality of W, X, Y, Z, where W and X were subject to Eq. 2a, and Y and Z were subject to Eq. 2b. For a more detailed description of the optimization, see Eyvindson et al. (2021). The multiobjective optimization represents a goal programming approach, integrating distance metrics to accommodate both equity and efficiency (Eyvindson 2012). Multifunctionality metrics were considered in the optimization at the landscape level rather than at the stand level. All indicators were first calculated at stand level, and then aggregated over the study area to produce a total sum over the landscape level.
We repeated the optimization with three management constraints: exclusively rotation forestry management regimes and set asides (RF scenario), continuous cover forestry and set asides (CCF scenario), or all forestry regimes allowed (ALL scenario). The harvest intensity gradient ranged from no income (0 NPI [€/ha], full set aside) to maximal NPI revenues (minimal set aside) in 21 steps with 5% stepwise increases in NPI. The share of set asides over the landscape increased with decreasing NPI values. The optimization resulted in 63 alternative management plans for RF, CCF and ALL management scenarios (Fig. 2b), where the ALL scenario is dominated by CCF regimes (RF share of 5-25%). Each management scenario represents a unique combination of regimes assigned to individual stands given specific objectives and constraints. Therefore, the management regime for each stand varied under alternative management scenarios, and harvest intensity scenario affected stands assigned as set asides. For a detailed description of the optimization and multifunctionality metric please see Eyvindson et al. (2021).

Probability of wind damage
We used estimates from the wind damage probability model (a binomial generalized linear model with logit-link function ) to predict the probability of wind damage for each stand, time step, and simulation scenario. The model predictors include dominant tree species, dominant tree height, time since thinning, maximal wind speed, temperature sum, presence of stand open edge, soil type, mineral soil depth and site fertility (Table 2). These predictors represent previously identified forest structures and management history affecting wind damage probability, and were available at both high resolutions and at a national scale ). Both wind speed and temperature datasets are available as open access datasets. We aggregated the predicted 10-years return levels of wind speed values (provided at a resolution of 20 m) to the stand level (mean: 10.9 m/s for 10 min., resolution 20 m, (Venäläinen et al. 2017)). The temperature sum represents a mean of the yearly temperature sums of the first 100 days warmer than 5 °C from 1985 to 2018 (mean 1219.7 over 100 degree days, Aalto et al. (2016)). The climatic datasets have a resolution of 10 km and therefore did not vary among stands for our watershed. Both wind speed and temperature datasets were assumed constant over the 100-year simulation time horizon, e.g., without considering the direct effects of climate change.
We linked the simulated stands structural parameters (tree species, tree height, timber volume, time since thinning) with their spatial geometry representation to identify neighboring stands over the landscape. We evaluated stands to have an open edge if its dominant tree height was 5 m higher than the dominant tree height of one or more neighbors , or if the stand was located at the landscape edge, i.e., neighboring a natural or artificial opening of size higher than 16 × 16 m 2 (e.g., road, meadow,

Data analyses
We calculated wind damage risk for each stand and time step along the harvest intensity gradient under the three management scenarios (RF, CCF, ALL). We further averaged the wind damage risk values over the landscape for actively managed and set asides over alternative scenarios. To validate our predicted wind damage risk values based on forest stand simulation (stand resolution, mean size 1.5 ha), we compared simulated values to published wind damage probabilities estimates using Wilcoxon rank-sum test in R software. As the published values (Suvanto et al. 2019, tile N4, freely available at https:// paitu li. csc. fi) were in raster format (resolution 16 m) and at higher resolutions, we first calculated their mean value per each stand to compare them with our stand-level wind damage risk prediction.
To relate wind damage risk to exposed timber volume, we calculated the average maximum stratum timber volume over each scenario, for actively managed stands and set asides. The maximum stratum timber volume represents the stratum of trees having the highest volume, e.g., the tallest trees in the stand. We assumed that higher exposed timber volume will link directly to expected wind damage volume (Zubizarreta-Gerendiain et al. 2018), which is the case especially in large windthrow events (Valta et al. 2019). Our model indicates the vulnerability of the stand to wind damage, but does not provide quantitative information related to the severity of wind damage in the stand. The estimated wind damage risk across management scenarios allowed us to compare management scenarios given prevailing management regime and the harvest intensity gradient.
To shed light on factors driving wind damage risk and wind-exposed timber volume (i.e., predictors dependent on the forest scenario in the wind damage probability model) for the given scenarios, we further investigated how scenarios differ in terms of average share of Norway spruce dominated stands, dominant tree heights, time since thinning, and height difference between neighboring stands along the harvesting intensity gradient.

Results
Our calculated wind damage probability values were consistent with published estimated probabilities for the selected watershed, validating our calculation of simulated datasets for their prediction over harvest intensity gradient (p value > 0.05, Wilcoxon test, Fig. 3, ).
Forest multifunctionality, wind damage risk, and exposed timber volume.
CCF and ALL outperformed RF in terms of landscape level multifunctionality over the entire harvest intensity gradient ( , Fig. 4a). Multifunctionality under RF continuously decreased with increasing harvest intensity. Conversely, ALL and CCF have had hump shape relationships, with the highest multifunctionality at moderate harvest levels (3 k€/ha to 6 k€/ha, Fig. 4a). Multifunctionality and wind damage risk diverged among management scenarios (Fig. 4b). At maximized harvest intensities (> 7.5 k€/ha), the wind damage probability of CCF increased to two times the mean wind damage risk in RF (1.7% in RF to 3.7% in CCF). The ALL scenario, as a combination of CCF and RF regime, was balanced between wind damage probability and multifunctionality, where a 25% increase in multifunctionality only slightly increased wind damage risks even at high harvest intensities (from 2.9 to 3.2%) and therefore limited the related wind damage risk.
Multifunctionality culminated at medium level of standing timber volume for ALL and CCF with typical humped shape (~ 160 m 3 /ha), and at no harvest at RF scenario (Fig. 4c). Landscape multifunctionality under RF increased with increasing timber volume.

Wind damage risk over harvest intensity gradient
Mean wind damage probability differed among management scenarios, and between actively managed stands and set asides, but showed a relatively consistent pattern over the harvest intensity gradient (Fig. 5a). Observed changes in wind damage risk over the harvest intensity gradient could

(c)
Harvest intensity (NPI k€/ha) 0 9 Scenario RF CCF ALL Fig. 4 Comparison of landscape level multifunctionality over a a harvest intensity gradient, b mean probability of wind damage risk, and c mean standing timber volume. The size of each dot is proportional to the NPI -net present income, with larger dot size representing a higher net present income and therefore higher harvest intensity (0 means no income, all stands are set aside). Scenarios: RF -rotation forestry, CCF-continuous cover forestry, ALL -all regimes allowed. Note that Multifunctionality is placed on the y-axis for be attributed to the actively managed part of the landscape as risk remained similar to the set asides across the three management scenarios (Fig. 5a). The wind damage probability was on average two times higher in CCF (4%) and ALL (3.6%) scenarios than under RF (1.7%) ( Table 3). This means that on average, four stands out of 100 per five-year period will face wind damage when continuous cover forestry is either the only or prevailing management regime. The corresponding figure is less than two out of 100 when rotation forestry dominates. Increasing harvest intensity only slightly affects wind damage probability in actively managed stands, with lower risks at low harvest intensities (< 1 k€/ha) for RF and CCF; and at minimized and maximized harvest intensities in ALL.

Timber volume at wind damage risk
Management scenarios differed in the amount of exposed timber volume in actively managed stands, although risk in set asides across management scenarios was relatively consistent (Fig. 5b). In actively managed stands, the exposed timber volume in RF was on average 1.3 times higher than under CCF and ALL scenarios, decreasing from 1.5 to 1.1fold with increasing harvest intensities. In CCF and ALL, average timber volume in actively managed stands remained relatively stable over the harvest intensity gradient, slightly increasing at maximized harvest intensity for ALL.
With increasing harvest intensity more set aside areas are transitioned to production areas, leading to a decrease in average timber volume for the remaining set asides (Fig. 5b). Therefore, high volume stands dominated by economically preferred species (e.g., Norway spruce, Fig. 6a), are selected for management, which is especially prominent in RF (Fig. 6a). This indicates that an increase in harvest intensity reduces the potential financial value of mean volume and dominant tree species of set asides. Timber volumes in set asides therefore decrease with increasing harvest intensity over all three management scenarios.

Predictor variables dynamics
The wind damage probabilities results were tightly linked with the changes in model predictors under varying management scenarios and harvest intensity. This included changes in species composition, tree heights, time since thinning, and tree height heterogeneity between neighboring stands (Fig. 6). We found little variability between scenarios for the proportions of Norway spruce (Fig. 6a), although Norway spruce stands were included as actively managed sooner along the harvest intensity gradient for RF (at 2.5 k€/ha) as compared with CCF and ALL (at 4 k€/ha). RF in actively managed stands differed from CCF and ALL scenarios in terms of average tree height, time since thinning, and mean height differences between neighboring stands (Fig. 6b-d). Actively managed stands under RF were shorter (mean 16 m and 21 m for RF and CCF&ALL, respectively), less frequently thinned (mean thinned every 15 years to < 10 years), and had higher height differences with neighboring stands over the whole intensity gradient than CCF and ALL. Increasing harvest intensity in RF increased mean height differences between neighboring stands and therefore increases landscape level roughness. The crossing of critical threshold value of 5 m height differences between neighbors already occurs at low harvest intensity (1.5 k€/ha). This contrasts with the relatively stable height difference between neighboring stands under dominant continuous cover management over the harvest intensity gradient (CCF and ALL, Fig. 6d).

Discussion
Large scale application of multifunctional forest management requires changes in management objectives and diversification of management approaches Eyvindson et al. 2021). However, we lack a full understanding how these changes in the operationalization of management practices will affect overall vulnerability to forest disturbances. Combining a forest growth management simulator (Rasinmäki et al. 2009) and wind damage risk prediction model , we found that continuous cover forestry is more efficient than rotation ing open edge-neighbors). Thinning in set asides is not allowed; therefore, the mean values are lacking in c). Scenarios: RF -rotation forestry, CCF-continuous cover forestry, ALL -all regimes allowed forestry in reducing the amount of exposed timber volume, while maintaining a high level of a diverse set of ecosystem services and biodiversity.

Effects of harvest intensification
We found that increased harvest intensity increases wind damage risk (our first hypothesis), depends on choice of management. Harvest intensification modifies forest wind damage vulnerability with increased vulnerability under continuous cover forestry scenarios and reduced vulnerability under rotation forestry scenarios (Fig. 4). To fulfill an increasing demand for raw timber material (FMME, FMAF, FME 2014;Hetemäki et al. 2017), variations in rotation and continuous cover forestry allow for shaping wind damage risk and exposed timber volume over time and space. Intensifying continuous cover forestry requires frequent thinnings (Fig. 6c, (Pukkala 2016)), whereas rotation forestry relies on shortening of rotation length and thinning frequency adaptation or exclusion to accumulate timber volume within the stand (Fig. 6b). Current forest management actions in Finland prioritize Norway spruce as the dominant species in the stand (Fig. 6d), while assigning stands dominated by tree species of lower economic importance and of smaller heights and volume are left as set asides (Fig. 6). Old-growth forests and forest structure heterogeneity have been largely reduced over the past decades in Finland due to intensive forest management (Henttonen et al. 2019). To address the loss of biodiversity while maintaining high levels of timber extraction requires increasing the share of protected areas (European Commission 2020) while also improving the network of retention forestry within commercial landscapes (Mielikäinen and Hynynen 2003;Kuuluvainen et al. 2019;Gustafsson et al. 2020). For example, successful conservation of deadwood dependent species requires larger and better connected set aside patches (e.g., > 1 ha for polypores (Ylisirniö et al. 2016)). These set asides should be of sufficient size that their interiors are buffered from adjacent managed stands, because stand microclimate can be affect multiple tree heights toward the stand interior (Dupont et al. 2011;Ylisirniö et al. 2016). Currently, set asides may strongly differ in conservation value between different forest holdings (Sverdrup-Thygeson et al. 2014).

Effects of scenarios on wind damage risk and timber volume
We found that rotation and continuous cover forestry have opposite effects on wind damage probability: rotation forestry reduced, whereas continuous forest cover increased predicted wind damage risk (Fig. 4b). Although the wind damage probability is higher for a given stand under continuous cover forestry (due to the presence of taller trees), than under rotation forestry (Fig. 6b), the results support our hypothesis that continuous cover forestry reduced the volume of exposed timber, and therefore reduced the potential economic loss. This is because fewer trees are found in stands under continuous cover forestry compared with rotation forestry, and therefore the risk of wind damage is ultimately spread across stands (Fig. 4c). This corresponds to the findings from Pukkala et al. (2016), that continuous cover forestry-as uneven-aged forestry-experienced a lower volume of damaged timber compared with rotation forestry following a windthrow event. In terms of wind damage probability, our findings indicate higher wind damage risk for continuous cover forestry than for rotation forestry seemingly contradicting Pukkala et al. (2016). The estimation of probability of wind damage depends on the investigated predictors and extent of observational data. Pukkala et al. (2016) investigated wind damage in stand interiors in experimental long-term study plots following a single wind event. Here, we applied the statistical model of Suvanto et al. (2019), based on occurrences (rather than damage extent) of windthrow events over Finnish national forest inventory plots. Pukkala et al. (2016) and Suvanto et al. (2019) investigated different aspects of wind disturbances, reflecting either their intensity over a short time period or incidence, respectively. Our study combined a forest growth and management simulator (Rasinmäki et al. 2009) with the wind damage risk model of Suvanto et al. (2019) over a century and across management scenarios; however, it omitted the occurrence of wind damage to the stands, which in reality would naturally remove trees from the forest stands and create deadwood over time (Zubizarreta-Gerendiain et al. 2018). The absence of wind disturbance in the simulation model caused continuous tree growth over a 100-year study period (unless harvested), and likely overestimates the exposed timber volume in absolute numbers (Fig. 4). The prediction of wind damage risk based on scenario-shaped forest structures over harvest intensity gradients however sheds light on potential risks in terms of forest economics and landscape multifunctionality.
Inclusion of the wind damage within simulation runs is possible by defining a critical wind speed at which the tree uproots or breaks (e.g., in HWIND model, Peltola et al. 1999b). This allows for the comparison between wind damage volumes over various management regimes (Zeng et al. 2007;Heinonen et al. 2009;Zubizarreta-Gerendiain et al. 2016), or to understand wind damage contribution to restore the habitat for deadwood-dwelling species (Zubizarreta-Gerendiain et al. 2018). However, simulating wind damage within forest simulations is still highly computationally demanding because it must account for stochastic wind events, their intensity, and forest structure under specific management regimes. Compared with simulating the incidence of wind events and forest induced wind damage (e.g., Zubizarreta-Gerendiain et al. 2018), our predictions of wind damage risk based on a statistical model ) allowed inclusion of a wider variety of landscape level management scenarios, over longer time periods, and larger study sites. Further combination of the statistical wind damage risk estimation with a mechanistic understanding about wind damage to the tree and stands level as applied in HWIND should be at the core of future studies investigating wind damage in uneven-aged forests (e.g., continuous cover forestry, set asides) and linking wind damage probability with induced wind damage extent.
We did not confirm our hypothesis that rotation forestry-due to increased forest edges-would consequently increase wind damage risk compared with continuous cover forestry, which has reduced height heterogeneity between neighboring stands (Fig. 6d, Pukkala et al. 2016). The small effect of open edges is likely due to the relatively small influence of edge effects on overall wind damage risk prediction (Suvanto et al. 2019, Table 2), and potentially because the average height difference between adjacent stands can be low in both types of forestry. A threshold of 10 m is often applied in predicting critical wind speed in temperate forests and in southern Sweden (Blennow and Sallnäs 2004;Seidl et al. 2014;Dobor et al. 2020), supported by empirical data following large-scale windthrow event (Bélouard et al. 2012). In boreal forests, stand edge is often delineated by the proximity of recently harvested or heavy thinned stands, creating a wind-vulnerable open gap (Zeng 2006;Zeng et al. 2007;Heinonen et al. 2009). To reduce wind damage risk in rotation forestry, it is crucial to reduce the length of exposed stand edges, e.g., by the aggregation of clear-cuts (Zeng et al. 2007), and minimizing the height difference between neighboring stands (Heinonen et al. 2009). Continuous cover forestry, creating over time uneven-age stand structure with several layers can limit potential wind damage into the stand interior (Hanewinkel et al. 2013), but might increase the probability of wind damage to the remaining dominant trees (Duperat et al. 2020). Another factor likely affecting our wind damage risk estimates was our models' lack of distinction between permanent (near lakes, meadows) and temporary open-edges (clear-cuts), although the former one has inherently higher wind stability ).

Dynamics predictors of wind damage risk
Both forest management scenarios and harvest intensity affected stand structural characteristics over the landscape (Fig. 6). The dominance of Norway spruce across scenarios and harvest intensity reflects its economic importance in the forestry sector. Norway spruce-despite its vulnerability to wind damage, pests, and fires at a mature age and under climate change (Venäläinen et al. 2020)-remains the preferred tree species in Finnish forest management. In addition, Norway spruce regeneration dominates in artificially planted stands under rotation forestry (Äijälä et al. 2014), and also under naturally regenerated stands managed under continuous cover forestry (Pukkala 2016). The increasing share of Norway spruce over the landscape can however further increase disturbance-related damage in future (Ikonen et al. 2020). Therefore, forest management strategies should consider optimizing the species composition of forest regeneration based on local site climatic conditions and soil fertility to minimize connected risks and future uncertainties under climate change (Torssonen et al. 2015).
At the same harvest intensity, rotation forestry reduced average tree height and increased height differences between neighboring forest stands compared with continuous cover forestry (Fig. 6d). Continuous cover forestry requires frequent thinning, whereas thinning can be excluded in low intensity rotation forestry. Exclusion of thinnings in rotation forestry can provide a mitigating solution between biodiversity conservation and timber provision due to higher deadwood volume in unthinned stands (Kuuluvainen et al. 2012). However, this is economically viable only at low harvest intensities (Fig. 6c), and therefore should be complemented by other management regimes.

Conclusion
Rotation and continuous cover forestry differed considerably in wind damage probability across the harvest intensity gradient. Rotation forestry lowers wind damage risk, in part by reducing average tree height, but accumulates high timber volume before the final harvest. This increases economic risks while also compromising non-timber ecosystem services and biodiversity. Continuous cover forestry increases wind damage risk due to taller trees and higher thinning frequencies, but nevertheless reduces the wind-exposed timber volume at any given time because standing timber volume is kept lower. However, a small share of rotation forestry within landscape level management planning contributes to overall high economic revenues and to landscape level structural diversity. A combination of both management approaches within the landscape thus provides an efficient way to reduce the amount of wind-exposed timber volume and increase habitat for vertebrate and non-vertebrate species in boreal forests.