Simulating wind disturbance impacts on forest landscapes: Tree-level heterogeneity matters

Environmental Modelling & Software 51 (2014) 1e11

Contents lists available at ScienceDirect

Environmental Modelling & Software journal homepage: www.elsevier.com/locate/envsoft

Simulating wind disturbance impacts on forest landscapes: Tree-level heterogeneity matters Rupert Seidl a, *, Werner Rammer a, Kristina Blennow b a b

Institute of Silviculture, Department of Forest- and Soil Sciences, University of Natural Resources and Life Sciences (BOKU), Vienna, Austria Department of Landscape Architecture, Planning and Management, Swedish University of Agricultural Sciences (SLU), Alnarp, Sweden

a r t i c l e i n f o

a b s t r a c t

Article history: Received 26 June 2013 Received in revised form 24 September 2013 Accepted 24 September 2013 Available online 17 October 2013

Wind is the most detrimental disturbance agent in Europe’s forest ecosystems. In recent years, disturbance frequency and severity have been increasing at continental scale, a trend that is expected to continue under future anthropogenic climate change. Disturbance management is thus increasingly important for sustainable stewardship of forests, and requires tools to evaluate the effects of management alternatives and climatic changes on disturbance risk and ecosystem services. We here present a process-based model of wind disturbance impacts on forest ecosystems, integrated into the dynamic landscape simulation model iLand. The model operates at the level of individual trees and simulates wind disturbance events iteratively, i.e., dynamically accounting for changes in forest structure and newly created edges during the course of a storm. Both upwind gap size and local shelter from neighboring trees are considered in this regard, and critical wind speeds for uprooting and stem breakage are distinguished. The simulated disturbance size, pattern, and severity are thus emergent properties of the model. We evaluated the new simulation tool against satellite-derived data on the impact of the storm Gudrun (January 2005) on a 1391 ha forest landscape in south central Sweden. Both the overall damage percentage (observation: 21.7%, simulation: 21.4%) as well as the comparison of spatial damage patterns showed good correspondence between observations and predictions (prediction accuracy: 60.4%) if the full satellite-derived structural and spatial heterogeneity of the landscape was taken into account. Neglecting within-stand heterogeneity in forest conditions, i.e., the state-of-the-art in many stand-level risk models, resulted in a considerable underestimation of simulated damage, supporting the notion that tree-level complexity matters for assessing and modeling large-scale disturbances. A sensitivity analysis further showed that changes in wind speed and soil freezing could have potentially large impacts on disturbed area and patch size. The model presented here is available as open source. It can be used to study the effects of different silvicultural systems and future climates on wind risk, as well as to quantify the impacts of wind disturbance on ecosystem services such as carbon sequestration. It thus contributes to improving our capacity to address changing disturbance regimes in ecosystem management. Ó 2013 Elsevier Ltd. All rights reserved.

Keywords: Forest disturbance Wind model Landscape modeling Windthrow Ecosystem heterogeneity iLand

Software availability The simulation software used here is iLand, the individual-based forest landscape and disturbance model. iLand is a research software tool written in Cþþ, using the open source and cross-platform Qt toolkit. The model runs on Windows, Linux and Mac platforms and includes a graphical user interface for visualizing and interacting with landscape processes in real time (e.g., observing the spread of wind disturbance). Additional JavaScript support provides fine-grained control of model behavior as well as the means for

automation of simulations, allowing great flexibility in applying iLand to diverse research questions. The software is open source and licensed under the GNU General Public License. Model code and executable can be obtained from the website http://iLand.boku. ac.at. The download package includes example files (e.g., species parameters, model drivers, landscape initialization) to run the model out-of-the-box. Furthermore, the iLand website hosts a comprehensive online model documentation based on Wiki technology. 1. Introduction

* Corresponding author. Tel.: þ43 1 47654 4068; fax: þ43 1 47654 4092. E-mail address: [email protected] (R. Seidl). 1364-8152/$ e see front matter Ó 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.envsoft.2013.09.018

Damage from windthrow and wind breakage is currently the most detrimental disturbance factor in Europe’s forest ecosystems

2

R. Seidl et al. / Environmental Modelling & Software 51 (2014) 1e11

(Gardiner et al., 2010). Moreover, wind disturbance regimes have intensified in Europe in recent decades (Schelhaas et al., 2003). This increase in frequency and severity has been attributed partly to changes in the climate system (Seidl et al., 2011a). Considering the climatic changes expected for the future, wind disturbance is expected to increase further in many areas of Europe (Blennow et al., 2010; Schelhaas et al., 2010), inter alia due to facilitated tree growth (taller trees are more susceptible to strong winds) and warmer and wetter winters (i.e., resulting in lowered tree stability to winter storms). Through a devaluation of the wood, the need to harvest prematurely, and negative market effects from large pulses of salvaged timber, wind disturbance can have a strong disrupting effect on timber production and the timber-based economy. For instance, a single wind event in southern Sweden in January 2005 (storm Gudrun) was estimated to have caused an overall economic damage of 2.4 billion Euros in forestry (Skogsstyrelsen, 2006a). For the same wind event, also distinct effects on the forest carbon (C) balance have been reported (Lindroth et al., 2009), with implications for the ability of forests to mitigate anthropogenic climate change. Furthermore, the strong winds of Gudrun have significantly reduced the growth of the remaining forests in southern Sweden in the first three post-storm years (Seidl and Blennow, 2012). These profound effects underline the growing importance of disturbances in the context of sustainable forest management, as an increasing frequency of severe events such as Gudrun have the potential to threaten the continuous supply of many important ecosystem goods and services to society. Consequently, disturbance management plays a central role in the development of strategies adapting forest management to climate change (see e.g., Kolström et al., 2011; Seidl et al., 2011b). Previous analyses have shown that management has considerable influence on the susceptibility of forests to strong winds (Albrecht et al., 2010; Valinger and Fridman, 2011), and that part of the recent intensification of disturbance regimes in Europe can be related to management-driven changes in forest structure and composition (Seidl et al., 2011a). These findings underline the fact, that besides being subject to external drivers such as climate, the disturbance regime is highly sensitive to forest management, and demonstrate the considerable potential for silvicultural risk management strategies. A powerful tool to harness this potential and make forest ecosystems more robust to changing disturbance regimes is modelbased scenario analysis. The modeling of wind disturbance in forest ecosystems has advanced considerably in recent years (see the reviews by Gardiner et al. (2008) and Seidl et al. (2011c)). Processbased wind risk models, estimating the probability of blowdown or breakage for certain stand conditions and wind speeds based on a combination of physical understanding and empirical approximations, have been successfully applied to address questions of wind risk in forest ecosystems (e.g., Gardiner et al., 2000; Blennow and Sallnäs, 2004; Zeng et al., 2004). Furthermore, linking such risk models with simulation models of stand development has enabled the assessment of dynamic interactions in time (e.g., through windinduced changes in stand structure, which in turn influences future susceptibility to disturbance) as well as the assessment of impacts on specific ecosystem services (e.g., Papaik and Canham, 2006; Schelhaas, 2008; Lagergren et al., 2012). Many previous approaches have focused on the stand scale, however, treating forest stands as widely homogenous entities and rarely accounting for within-stand variation in e.g., tree properties and spatial distribution. Such model designs are thus not well suited to address spatially more complex and heterogenous management regimes such as continuous cover forestry or group selection. More broadly, they are at odds with recent developments towards managing for complexity in silviculture (Puettmann et al., 2009). A second limitation of many previous wind modeling

approaches is related to the consideration of wind disturbance as a singular, discrete event rather than a dynamic process over time. Storm events in Europe typically last from a few hours to many days (see e.g., Schiesser et al., 1997), with varying intensities within this time span. Little is known how wind loadings on trees change during a storm event. However, when trees start to fall, previously sheltered trees become exposed, and the fetch (i.e., the stretch of open upwind area) increases. Stand-level models are generally not well suited to incorporate such effects due to their limited spatial extent (but see the recent work by Byrne and Mitchell, 2013). However, also many current landscape modeling approaches do not incorporate spatio-temporal interactions between wind and vegetation in a process-based manner but mostly use prescribed disturbance sizes and shapes and contrast them with the prevailing vegetation to derive disturbance impact (e.g., Scheller et al., 2011). In summary, while many process-based approaches to simulating susceptibility to and occurrence of wind damage exist, models that account for tree-level heterogeneity while considering the dynamic landscape context are still widely missing to date. Here, our objectives were (i) to develop a novel wind modeling approach that simulates wind impact at the level of individual trees, and dynamically accounts for changing stand structure during a wind event, and (ii) to evaluate the new tool for a landscape in Sweden against independent empirical data from the storm Gudrun in January 2005, particularly analyzing the effect of tree-level heterogeneity in simulations of wind impact. More broadly, our aim was to present a process-based yet parsimonious model that simulates wind disturbance as an emerging property of dynamically changing stand conditions, environmental drivers (e.g., wind events, climate change), and the landscape context (e.g., terrain). 2. Methods and materials 2.1. Simulation platform: iLand We used iLand, the individual-based Landscape and disturbance model (Seidl et al., 2012a), as the simulation platform for the development of our wind disturbance model. iLand was chosen since it is spatially explicit and operates on the level of individual trees, allowing for the simulation of heterogenous stand conditions. Furthermore, iLand is designed to simulate dynamic ecosystem processes at the landscape scale, and thus enables us to consider spatial context and disturbance spread explicitly. The following section gives a brief overview over iLand, for more details please refer to Seidl et al. (2012a,b). iLand was developed to dynamically simulate the interactions and feedbacks between climate (i.e., external drivers and potential changes therein), forest vegetation (i.e., tree-level demographic processes such as growth, mortality, and regeneration), and disturbance regimes (i.e., large-scale mortality agents and their interactions in space and time). It is a hierarchical multi-scale model, and explicitly considers bottom-up emergence as well as top-down constraints of processes at different scales (tree, stand, landscape). Trees are simulated as adaptive agents that compete for resources (predominately light, but also water and nutrients are considered), and dynamically adapt to their environment (Seidl et al., 2012a). To aid scalability of individual-tree interactions in landscape-scale applications a patternbased rendering of ecological field theory (Wu et al., 1985) is used in iLand. A measure of competitive strength is derived for every individual and used to partition the stand-level resources available. Production physiology is modeled in a simplified process-based manner using a light use efficiency approach (Landsberg and Waring, 1997), scalar modifiers for the environmental effects on utilizable radiation, and an allocation regime based on allometric equations (Duursma et al., 2007). Tree mortality is based on C starvation (stress-related mortality) as well as age (life-history traits), besides considering disturbance-related tree death. iLand keeps track of dead organic matter in four pools (snags, downed woody debris, litter, soil organic matter), simulating decomposition, heterotrophic respiration, and nitrogen feedbacks to trees following the approach of Kätterer and Andren (2001). Seeds (spatially explicit dispersal in the landscape), favorable environment (phenology approach following Nitschke and Innes, 2008), and light (cf. Seidl et al., 2012a) are required for trees to regenerate in the model. Sapling growth and competition is modeled for mean trees at a 2  2 m resolution, from which trees are recruited into the individual-based structure of iLand after reaching a height threshold of 4 m (Seidl et al., 2012b). The model has been successfully applied over extensive environmental gradients in ecosystems in western North America and central Europe (Seidl et al., 2012a), and was recently used to quantify landscape-level forest C budgets in complex terrain (Seidl et al., 2012b).

R. Seidl et al. / Environmental Modelling & Software 51 (2014) 1e11 2.2. Wind disturbance modeling 2.2.1. Disturbance initiation and occurrence The iLand wind disturbance module presented here is designed as a processoriented model iteratively simulating the impact of a storm event in a forest landscape (see Table 1 for an overview). The occurrence of wind events is therefore simulated based on wind data (mean hourly wind speed, wind direction), which can be derived either from weather stations in the vicinity of a study landscape or from simulations with regional climate models (e.g., Blennow et al., 2010). Temporal scenarios of different wind climates can be accommodated in the simulations by either supplying the model with wind event time series from climate models, or by generating such series by drawing from empirically observed wind distributions. Spatial differences in wind exposure in the landscape, e.g., as an effect of topography, are accommodated by means of a grid of modifiers relative to a reference point (e.g., weather station location), which can be derived either from detailed airflow modeling or by means of statistical approaches. Based on these wind driver data, the actual occurrence and extent of damage are simulated as emergent properties of the model, as described in the following sections. We assume that major wind disturbance events are initiated where canopy rugosity changes abruptly, i.e., where vertical differences between the top heights of neighboring grid cells exceed 10 m in the simulated forest (see also Blennow and Sallnäs, 2004). Cells with top height <10 m are not considered to be directly damaged by wind. A cells’ Moore neighborhood (i.e., its eight neighboring cells) is used to detect changes in canopy rugosity. The wind damage simulation sequence described in detail in the following sections is executed for all the cells in the landscape that are identified as having vertical differences >10 m. These cells are subsequently referred to as “edges” in this contribution. A storm event is simulated iteratively (Table 1), i.e., after each sequence of wind disturbance calculations the forest structure is updated in the model, new edges are detected, and a new iteration of wind damage calculations is initiated (see also Byrne and Mitchell, 2013). This iterative simulation procedure is stopped after the maximum event duration (iteration count) has been reached (see Section 2.2.7 below). 2.2.2. Vertical wind profile To derive the canopy top wind speed from above canopy wind speed a vertical wind profile is calculated, following the approach of ForestGALES (Gardiner et al., 2000). Zero-plane displacement (d0), i.e., the height at which the wind force acts on the forest, is calculated according to Raupach’s drag partitioning model (Raupach, 1992, 1994, Eq. (1)). 1  eðcdl lÞ

0:5

d0 ¼ h$ 1 

! (1)

0:5

ðcdl lÞ

3

with rz the crown radius, hz the height to the live crown, kc an empirical coefficient related to the crown shape, p a constant accounting for porosity (0.5), and n the number of trees per 100 m2 cell. In iLand, rz is available from light competition calculations (Seidl et al., 2012a), while hz is calculated from a fixed, species-specific ratio (also used in the ray tracing algorithm applied in the model). For wind disturbance calculations every 10 m grid cell is represented by the tallest tree occupying it. In order not to overestimate the effect of smaller trees we scale n to the number of idealized trees of equal size as the focal tree, using tree basal area as scaling criterion. Note also that n does not include trees in the regeneration layer (i.e., <4 m height). An application of these stand level methods to sub-stand grid cells was recently successfully demonstrated by Byrne (2011). Surface roughness (z0) is subsequently calculated as Eq. (3), following Raupach (1994): z0 ¼ ðh  d0 Þ$ekgþ0:193

(3)

with k the von Karman constant (0.4) and g the drag coefficient, the latter defined as (Eq. (4))

g¼ 

1 1 ; max l ¼ 0:6 2 Cs þ Ce $2l

(4)

with Cs the surface drag coefficient (0.003) and Ce the element drag coefficient (0.3). The wind speed at the cell top height (Uc) relative to the mean hourly reference wind speed 10 m above d0 (U10) is subsequently calculated according to Eq. (5).

Uc ¼ U10

ln



ln

hd0 z



 0 10 z0

(5)

2.2.3. Individual-tree turning coefficient We build our calculations of tree response to wind on the recent findings of Hale et al. (2012), who showed that e over a range of stand and soil conditions and for two different species (deciduous and evergreen conifers) e the ratio between the maximum turning moment (Mmax) and the square of the hourly wind speed at the top of the canopy (U2c ) was conservative over a wide variety of tree sizes (Eq. (6)). For well acclimated (i.e., not recently thinned) stands they found Mmax ¼ Uc2 $Tc

(6)

with

with h the top height of a given 10 m cell, cdl an empirical parameter (7.5), and l the frontal area index, being defined according to Eq. (2).

2$rz $ðh  hz Þ$kc$p l ¼ 100=n

(2)

Tc ¼ 117:3$dbh2 $h where dbh is the tree diameter at breast height. The turning moment affecting trees is significantly higher at stand edges than within a forest stand. Using data from Gardiner et al. (1997), Byrne (2011) established an equation relating the maximum turning moment to the distance from the stand edge. The factor fedge represents a

Table 1 The sequence of steps to iteratively simulate the impact of wind disturbance events in iLand. The spatial grain of the calculations is 10 m grid cells. See Sections 2.2.1e2.2.7 for details. #

Step

Description

(1)

Start wind disturbance simulation

(2)

Detect edges in the landscape

(3)

Calculate vertical wind profile

(4)

Calculate individual-tree turning coefficient

(5) (6)

Calculate critical wind speeds Simulate wind impact on vegetation

(7)

Maximum event duration reached?

(8)

End wind disturbance simulation

The occurrence of a wind event can be triggered in three different ways in the model, corresponding to different levels of data available for describing the wind regime of a landscape: First, drawing from an empirical wind probability distribution and comparing to a random number (see Scheller et al., 2011); second, using wind speed input data (from either observations or simulations) and evaluating if the simulation year exceeds a predefined percentile threshold for damage (Klawa and Ulbrich, 2003); or, third, deterministically initiating the wind simulations in the model to study particular wind events (this study). Changes in canopy rugosity (i.e., cells with a height difference of >10 m to their neighboring cells) are identified in the landscape. The subsequent steps (3)e(6) are conducted for the thus identified edge cells only. A vertical wind profile at the stand edge is calculated to derive the wind speed at canopy top height. Calculations are conducted for the tallest tree of every edge cell. The individual-tree turning coefficient is calculated following the approach of Hale et al. (2012). Upwind gap size is calculated (relative to tree height) from analyzing grid cells in upwind direction, and sheltering and support from immediate neighbors is taken into account via the iLand light competition index LRI (Seidl et al., 2012a). From the results of steps (3) and (4) critical wind speeds for uprooting and stem breakage are calculated. If the wind speed of the current iteration exceeds the critical wind speed the trees on the focal cell are uprooted or broken. Stand structure and biomass pools are updated accordingly. If the iteration count is below the maximum iteration count (representing the event duration), the simulation is continued by iterating back to step (2), else step (8) is executed. A summary over the iterative wind disturbance simulation is generated (e.g., trees damaged, disturbance size, etc.) and written to the models’ output database. The wind module is terminated and the simulation of forest landscape dynamics is continued from the new, wind-modified vegetation state.

4

R. Seidl et al. / Environmental Modelling & Software 51 (2014) 1e11

(constant) ratio between the maximum turning moments at the stand edge and conditions well inside the forest. Mmax ¼ Uc2 $Tc $fedge

equation, but that increasing competition decreases Tc due to an increasing effect of sheltering from neighbors. In the context of thinning and forest management it is also interesting to note that the absolute sensitivity of Tc to competition increases with tree size (Fig. 1), but that the relative sensitivity (expressed in percent Tc response) decreases with tree size. To apply this equation in the context of iLand a transfer function between CHegyi and LRI (Light Resource Index, the process-oriented competition index used in iLand, see Seidl et al., 2012a) was devised (Eq. (10)).

(7)

with fedge ¼ 5. 2.2.4. Neighborhood effects on tree stability The effect of local upwind vegetation from the edge, also known as fetch, is calculated following data from Gardiner et al. (1997) using the equations of Peltola et al. (1999). Their approximation assumes that the effect of upwind gaps increases with gap size up to 10 tree heights in upwind direction. The thus derived gap factor (fgap) assumes 0.2 for closed stands (and thus cancels out with fedge in closed stands) and reaches fgap ¼ 1 for cells with upwind gaps of 10 tree heights in size. Gaps in iLand are defined as cells with a canopy height of <10 m than the focal cell. The relative gap length rg (gap length/canopy height of the focal cell) is determined by scanning the canopy height in upwind direction using Bresenham’s line drawing algorithm. The maximum turning moment calculation of Eq. (7) is subsequently modified by fgap (Eq. (8)). The additional moment provided by the overhanging displaced mass of the canopy was considered via the factor fCW, set to 1.25 here (see Gardiner et al., 2000). Mmax ¼ Uc2 $Tc $fedge $fgap $fCW

CHegyi

10

2.2.5. Critical wind speeds Critical wind speeds (cws) are calculated separately for uprooting and stem breakage, following the approach of Gardiner et al. (2000). The critical moment for uprooting is given by (Eq. (11))

In addition, the effect of shelter and support from neighboring trees (or lack thereof) needs to be considered (Schelhaas et al., 2007). Hale et al. (2012) proposed the novel idea of using widely applied growth competition indices for this task, and showed a strong relationship of Tc to different competition indices. This is important because Eq. (6) is parameterized for well-acclimated (i.e., not recently thinned) stands, and a thinning would affect the shelter and support from neighbors, but not result in immediate changes in tree size (i.e., the sole explanatory variable of Tc used in Eq. (6)). For dynamic simulation models such as iLand, which are developed to inter alia simulate forest management, it is thus relevant to additionally control for this effect of local sheltering and density. We reanalyzed the data of Hale et al. (2012) to include the effect of local shelter from neighboring trees in the turning moment coefficient equation (Eq. (6)). In this regard it has to be noted that competition indices are by their very nature correlated with tree size. We thus used principal component regression for our reanalysis, i.e., we calculated regression coefficients for the orthogonal transformations of the predictors (i.e., the principal components), and then retransformed these to the original predictors (Morzuch and Ruark, 1991). Also, since the relationship between tree size and competition changes with stand development we included an interaction term in our reanalysis. The resulting equation and coefficients are given in Eq. (9), and its sensitivities are displayed in Fig. 1.

Mcrit ¼ Creg $SWgreen

11

with Creg a species-specific, empirical coefficient derived from tree pulling experiments (e.g., Nicoll et al., 2006) and SWgreen the green stem weight. Given the information is available, also soil-type specific differences in Creg can be incorporated as a modifier to the species-specific base value. The critical wind speed (cws) at canopy top necessary to achieve this moment is derived by equating Eq. (11) with Eq. (8), and transforming for wind speed (Eq. (12)).

Creg $SWgreen Tc $fgap $fedge $fCW

cwsuprooting ¼

!0:5 12

The critical bending moment at which stem breakage occurs was defined by Gardiner et al. (2000) as (Eq. (13)) Mcrit ¼

(9)

with CHegyi the distance-dependent competition index by Hegyi (1974), and with dbh and h both in units of meters. From Fig. 1 it is clear that size dominates the refitted Tc

p 32

$fknot $MOR$dbh3

13

with MOR the species-specific modulus of rupture, and fknot an empirical factor accounting for reduced stability due to the presence of knots (currently set to 1, see

400 200

Tc (kg)

400 200

600

Tc (kg)

400 200

600

95th

600

800

5th

800

800

1000

(c) 1000

(b)

1000 1200

(a)

Tc predicted (kg)

LRI  0:05 0:05 < LRI < 0:5 LRI  0:5

If, for instance, a thinning intervention is simulated with the model, the LRI of the remaining trees increases (i.e., the available light resources increase), which in turn leads to a decreased sheltering effect immediately after the thinning, and an increase in Tc (Eq. (9)). In the years following a thinning, the liberated trees will increase their size and crown volume, and will gradually begin to compete again (i.e., experiencing a decreasing LRI over time), which in the context of wind disturbance will result in an increasing sheltering effect.

(8)

Tc ¼ 4:42 þ 122:1$dbh2 $h  0:141$CHegyi  14:6$dbh2 $h$CHegyi

8 < 3:41 ¼ 3:733  6:467$LRI : 0:5

95th

0

0

0

5th

0

200

400

600

800

Tc observed (kg)

1200

0

2

4

6

dbh²h (m³)

8

10

1

2

3

4

CHegyi

Fig. 1. (a) Observed versus predicted turning coefficients (Tc) for the refitted Hale et al. (2012) model explicitly including the effect of competition (Eq. (9), R2 ¼ 0.950). (b) Sensitivity of Tc to tree size. Solid line is for mean competition index (CHegyi) of the Hale et al. (2012) dataset, dashed lines are for the 5th and 95th percentile of CHegyi, respectively. (c) Sensitivity of Tc to competition as represented by the distance-dependent competition index of Hegyi (CHegyi). Solid line is for the mean tree size in the Hale et al. (2012) dataset, dashed lines give results for the 5th and 95th percentile of tree size, respectively.

R. Seidl et al. / Environmental Modelling & Software 51 (2014) 1e11 also Byrne (2011)). Analog to uprooting the critical wind speed for breakage was derived according to Eq. (14):

cwsbreakage ¼

MOR$dbh3 $fknot $p 32$Tc $fgap $fedge $fCW

!0:5 14

2.2.6. Wind impact on tree vegetation The calculations outlined in the previous sections are conducted for the dominant tree of a 10 m cell, and wind damage occurs if the wind speed exceeds the critical wind speed. Trees are uprooted if cwsuprooting < cwsbreakage or are broken otherwise. If the soil is frozen, as determined by a soil temperature 0  C at 10 cm soil depth, only stem breakage is assumed to occur (see Peltola et al., 2000; Blennow et al., 2010). Soil temperatures are determined for the particular day of year of a wind event from air temperature in conjunction with vegetation and litter characteristics following Paul et al. (2004) in iLand. If the dominant tree of a 10 m cell is thrown, all smaller trees on the cell are likewise assumed to be killed (i.e., assumed to being either hit by the falling tree, or uprooted together with the root plate of the dominant individual). The biomass from uprooted trees is moved directly to the downed woody debris (stem and branch compartments) and litter (foliage compartment) pools, i.e., no snags are created from windthrown trees. For trees broken by wind it is assumed that other individuals <4 m height on the same cell survive, while all taller neighbors are also assumed to be killed by the broken tree. Branches, foliage, and half of the stem compartment biomass are directly moved to the detritus compartments of the iLand soil module, while the remaining half of the stem biomass is assumed to remain as snag. 2.2.7. Dynamic spread of wind damage Damage from a wind event is simulated iteratively in iLand (Table 1). Within a wind event, windthrow and wind breakage create new edges and increase fetch, thus facilitating further wind damage. This process is dynamically simulated in iLand by looping through the simulation sequence described above, i.e., newly created edges are in turn evaluated for exceedance of critical wind speeds and further wind damage. To account for varying wind intensity during a storm event, wind speeds are varied randomly within a certain range with every iteration to mimic the temporal variation in wind speeds. The iterations are stopped when a specified event duration has been exceeded (input parameter). Storm durations can be derived either from climate modeling (e.g., Beniston et al., 2007; Pryor et al., 2012) or drawn from a distribution of historically observed events (see e.g., Schiesser et al., 1997). In summary, by combining wind exposure with detailed information on vegetation structure and composition (i.e., susceptibility), the size, spatial footprint, and severity of a wind event is simulated as an emergent property of the system in iLand.

2.3. Storm Gudrun and the Växjö landscape To test the newly developed wind model and in particular analyze the effects of resolution and spatial heterogeneity in wind disturbance modeling we studied the effect of the storm Gudrun (January 8the9th, 2005) at a 1391 ha forest landscape near the city of Växjö, Kronobergs län, Sweden. Gudrun was the most detrimental storm ever recorded in Sweden, damaging approximately 75 Million m3 of timber (Skogsstyrelsen, 2006a). Peak gust wind speeds reached up to 36 m s1 over land, and the storm endured for approximately 12 h (SMHI, 2005). Kronobergs län in southern central Sweden was the worst hit area, with 16.9% of the standing timber stock damaged by Gudrun (Skogsstyrelsen, 2006a). 2.3.1. Study landscape The 1391 ha Växjö landscape, situated approximately 15 km northwest of the city of Växjö, Sweden (56 530 N, 14 470 E), is characterized by mild, rolling topography with an elevation range of 116e232 m asl. It features boreal conditions and is dominated by Norway spruce (Picea abies (L.) Karst., 54.2%), Scots pine (Pinus sylvestris L., 22.8%), and birch species (Betula pendula Roth. and Betula pubescens Ehrh., 12.4%). The landscape has been managed in a clear-cut system with a portion of the landscape being regenerated via shelterwood. The mean stand age prior to the storm was 46 years, with 21.0% of the landscape older than 60 years. Pre-storm site index data was available from general terrestrial forest inventory. Data on standing volume, species composition, tree height, and age was derived from a satellite-based kNN imputation (25 m horizontal resolution, resampled to 10 m grids) of forest inventory data (Reese et al., 2003). These kNN data were used to identify stands, and to estimate stem numbers, basal area and dbh using relationships by Bondesson (in Karlsson and Westman, 1987). Furthermore, information on within-stand heterogeneity in tree height was extracted from the kNN imputation of forest inventory data (Reese et al., 2003). The sizes of the 297 stand compartments in the landscape ranged from 0.2 to 45.4 ha, with a mean stand size of 4.9 ha. The mean  standard deviation of canopy top height over the landscape was 18.0 m  6.8 m, and the average stem number was 1352  717. Since 96% of the landscape has the same soil type according to the soil map of the Geological Survey of Sweden (http://www.sgu.se) we assumed homogenous soil conditions in this study. Furthermore, based on an analysis of the weather

5

conditions in the three weeks prior to storm Gudrun we assumed soil conditions not to be frozen. According to a satellite-derived classification 21.7% of the landscape (302 ha) were damaged by the strong winds of January 8the9th 2005 (Skogsstyrelsen, 2006b). Most affected were two distinct areas in the center of the landscape which are also the parts of the landscape with the strongest topographic gradients and highest elevation. The median size of damaged patches was 0.38 ha, the mean size was 1.96 ha. The largest contiguous windthrown patch was 116.3 ha in size. Taller (older) forests were disproportionally affected by the storm: 42.2% of forests with a dominant height of >27 m were disturbed, while only 3.4% of forests 11 m experienced wind damage. These detailed satellite-derived damage data (Skogsstyrelsen, 2006b) were used as independent evaluation dataset for the newly developed wind model in this study. 2.3.2. Wind data Gust and ten minute wind speed data for Gudrun were available for the Växjö weather station with hourly temporal resolution (SMHI, 2005). In order to convert from ten minute wind speed data to mean hourly wind speeds we used a factor of 0.95, based on an in-depth analysis conducted by Larsén and Mann (2006) for sites in neighboring Denmark. We temporally defined the Gudrun storm event as the twelve hour period from 1800 hrs on January 8th to 0600 hrs on January 9th 2005 for which the mean hourly wind speed at the Växjö weather station exceeded 10 m s1. The mean hourly wind speed over this twelve hour period was 13 m s1, the maximum mean hourly wind speed reached 16 m s1. To account for changing wind speeds within the twelve hour storm period we randomly varied the mean wind speed of 13 m s1 by 25% in every iteration of the simulation. The main wind direction during Gudrun was west (259  22 ) at the Växjö weather station. Based on data of maximum contiguous windthrown area and storm duration (Nilsson et al., 2007) we estimated the maximum number of iterations for the simulations of the twelve hour storm event to be 108 (i.e., an assumed maximum spread rate of 90 m h1). In order to use this information, recorded for a single weather station location in the vicinity of our study area, at the landscape scale and account for the effect of local topography on wind exposure we here used a simple index of topographic exposure, the topex-to-distance (Hannah et al., 1995; Quine and White, 1998). Hannah et al. (1995) found that observed wind speeds are strongly related to topexto-distance (TTD) values in Scotland, and Ruel et al. (2002) also report good results for applications in Canadian forests. Notwithstanding the higher process resolution of detailed airflow models (e.g., Blennow and Sallnäs, 2004) the TTD was chosen due to its simplicity and ease of parameterization (for a comparison of TTD to more complex airflow models see for instance Suarez et al., 1999). We used wind data from 12 weather stations from southern and central Sweden (SMHI, 2012) in combination with the 30 m resolution digital elevation model of ASTER (2011). For every weather station location the TTD was calculated for distances of 0.5 km, 1 km, 2 km, and 3 km, and related to the 95th percentile of mean ten minute wind speed recorded in 3 h intervals over ten years (1996e2006). TTD for individual distances were only weakly related to wind speed. However, a linear combination of all TTD indices and elevation (using principal component regression to control for correlation between predictors) explained 70.6% of the variation in 95th percentile wind speed. From this relationship, we calculated a topography-dependent windiness modifier by relating the TTD-derived windiness of the reference weather station (i.e., Växjö airport) to the corresponding value at each grid cell of the study landscape. In the simulations this modifier was used to multiplicatively adjust the weather station wind data to the local conditions of the simulated grid cells.

2.4. Study design To evaluate the newly developed model (cf. Bennett et al., 2013) we conducted three simulation experiments. First, to investigate the models behavior over a range of stand conditions, we conducted a local (i.e., “one-factor-at-a-time”) sensitivity analysis of critical wind speeds over a range of stand conditions and tree characteristics (species: Norway spruce). This analysis was particularly designed to elucidate the relative effect of local sheltering in relation to more traditional size- and structure-related susceptibility criteria (see also Schelhaas et al., 2007). Second, we simulated storm Gudrun in the Växjö landscape with the new iLand wind model, and compared simulation results to the satellite-based wind damage observed for the landscape. We analyzed model performance at two different spatial resolutions (10 m and 100 m grid cells) and calculated both non-spatial performance measures (percent of landscape damaged, odds ratio of predicted vs. observed damage, prediction accuracy) as well as spatially explicit measures of correspondence between simulated and observed wind damage (conditional Kappa). In our wind simulations, higher levels of spatial and structure variability increase canopy rugosity and thus the number of starting points for wind damage, modulate the competitive relationship between trees and consequently also the local sheltering against wind (e.g., higher variability in within-stand tree heights will mean less competition and sheltering for dominant trees), and in case of within-stand gaps also elevates the gap factor. In order to elucidate the role of tree-level heterogeneity on the simulation of wind disturbance in forest ecosystems we conducted analyses for three different model initializations: (i) A strict stand-level initialization of the landscape, in which neither structural heterogeneity (i.e., variation in tree size

6

R. Seidl et al. / Environmental Modelling & Software 51 (2014) 1e11

within a stand) nor spatial heterogeneity (i.e., within-stand clumping, clustering, gaps) are accounted for, served as the starting point for our analyses (variant S00). Identical mean trees of every species are positioned at equal spacing across the stands in this variant. Variant S00 thus resembles the representation of forest stands in many established process-based wind risk models (e.g., Gardiner et al., 2008), which conduct calculations for a mean tree per stand. (ii) Variant S10 accounted for the full structural variation in within-stand tree size as recorded in the kNN imputation data, but did not account for within-stand heterogeneity in space. Here, trees within a stand were assigned positions quasi-randomly, using basal area as a weighting criterion to ensure that tree positions were selected realistically (i.e., bigger trees are wider spaced than smaller trees). The same tree positions were used in all iterations of variant S10. (iii) The third initialization variant (S11) accounts for the full structural and spatial heterogeneity within stands. Here, tree positions were determined to resemble the kNN dominant height layer. In other words, the distribution and location of tall and small trees, gaps, and densely populated patches within a stand followed the remotely sensed pattern. S11 thus makes use of the full suite of available data (stand-level inventory, spatial and structural within-stand heterogeneity from satellite-derived kNN) and is the most realistic (and most heterogenous) initialization variant. Note, however, that the average properties of the 297 stands were the same in all variants, and that S10 and S11 differ only in the spatial arrangement of trees (while the simulated individuals are the same). The different levels of heterogeneity studied in the three alternative initialization variants are graphically displayed in Fig. 2. As a third and final analysis step we investigated the sensitivity of the simulation results for the Växjö landscape to different storm characteristics. Using the Gudrun event as a starting point, we varied important drivers such as wind speed, wind direction, and event duration by 20% to assess the effect different storms might have had on the landscape. We furthermore also conducted a run in which we assumed the soil conditions to be frozen, to evaluate the effect of soil frost on damage from wind disturbance. These analyses are designed to give further insight into the sensitivities of the model, but are also indicative for the disturbance response that can be expected under a changing climate and wind regime. For all analyses, ten replicates were simulated and mean values are presented except where stated otherwise. Creg was set to 111 and 134 Nm kg1 for Norway spruce and Scots pine, respectively, and the MOR values assumed for these species were 30.6 and 39.1 MPa (Gardiner et al., 2000). The analysis of all simulation results was conducted using the R project for statistical computing (R Development Core Team, 2011), in particular using the libraries sp, spdep, and SDMTools (Bivand et al., 2008; VanDerWal et al., 2011).

3. Results 3.1. Sensitivity analysis of critical wind speeds In general, the cws (i.e., the wind speed required to uproot or break a tree) for stem breakage was higher than for uprooting in the

model, with the exception of tall trees with very small diameters (Fig. 3). The model was sensitive to increasing tree height as well as to increasing height to diameter ratio, with the cws decreasing with both factors. The effect of local sheltering, here incorporated via the iLand competition index LRI, was significant for both uprooting and breakage, but overall of a lesser magnitude than the effects of tree size and shape. Generally, locally sheltered trees had higher cws in the model compared to openly growing individuals. However, more influential than the effect of local sheltering and support from neighboring trees was the effect of upwind gaps (i.e., broader-scale shelter). 3.2. Simulating the impact of storm Gudrun Neglecting within-stand heterogeneity (S00) resulted in a considerable underestimation of wind damage in the simulations. Prediction accuracy was only moderate, and the spatial damage patterns emerging from the simulation were only weakly correlated with the satellite-based observations in variant S00 (Table 2). Accounting for structural variation within stands (S10) increased the simulated wind damage and improved predictions both in terms of quantity as well as spatial patterning. However, performance improved even more strongly if the full structural and spatial within-stand variation was considered (S11). Both overall damage level as well as spatial patterning were strongly improved by also accounting for within-stand spatial heterogeneity. An exception to this finding were severely damaged areas, for which the correspondence between simulation and observation did not improve further after accounting not only for structural (S10) but also for spatial heterogeneity (S11) in the simulations (Table 2). The spatial pattern of wind damage simulated under variant S11 matches the observations well (Fig. 4). The model was able to reproduce the two main damage clusters in the Växjö landscape, but overestimated damage along the western (i.e., windward) border of the landscape. Furthermore, the size distribution of disturbed patches emerging from the iterative simulations corresponded well with observations (Fig. 5a). With regard to impact iLand showed a tendency of overestimating wind damage for small trees and underestimated damage for tall trees (Fig. 5b).

Fig. 2. Maps of the model initialization variants studied for the 1391 ha Växjö study landscape (grain: 10 m grid cells). Dominant tree height (m) prior to the storm Gudrun is indicated, derived from a combination of forest inventory and kNN imputation. S00: neither structural nor spatial within-stand variation are accounted for e trees were initialized as average tree per stand with regular spacing; S10: trees were initialized with full within-stand structural variation (based on kNN data) but random spatial variation (i.e., with tree positions assigned quasi-randomly); S11: full structural and spatial variation, both size variation and tree positions were derived from kNN data.

R. Seidl et al. / Environmental Modelling & Software 51 (2014) 1e11

7

(b) 25

25

(a) h=25m

dbh=30cm

br 20

20

ur

15

br ur 10

cws (m s−1)

15

ur

br ur

5

5

10

cws (m s−1)

br

10

20

30

40

50

20

25

dbh (cm)

35

40

(d) 25

25

(c)

dbh=30cm, h=25m

br ur

15

br ur

5

5

br ur

10

br ur

cws (m s−1)

15

20

20

dbh=30cm

10

cws (m s−1)

30 h (m)

50

60

70

80

90

100

110

120

h dbh ratio (dim.)

50

100

150

200

250

upwind gap size (m)

Fig. 3. Sensitivity of critical mean hourly wind speed (cws) ten meters above zero plane displacement height to (a) diameter at breast height, (b) tree height, (c) height to diameter ratio, and (d) upwind gap size. Solid lines are for trees at the edge of dense stands (LRI ¼ 0.05), dashed lines for individuals without shelter and support from neighbors (LRI ¼ 1.0). ur ¼ cws for uprooting, br ¼ cws for stem breakage. Displayed sensitivities in panels (a) through (c) assume an upwind gap size of 150 m. Species parameters are for Norway spruce.

3.3. Sensitivity to storm characteristics An analysis of the sensitivity to storm characteristics (initialization variant S11) showed that wind speed was by far the most important driving variable (Table 3). A 20% increase in wind speed would have almost doubled the average size of wind-damaged patches at the Växjö landscape. Frozen soils, on the other hand, would have reduced the damaged area considerably. In contrast, the sensitivity to wind direction and event duration was small. Further analysis showed that most of the damage occurs within the first third of the simulation (Fig. 6), which explains the low sensitivity to the duration of the storm event. 4. Discussion 4.1. Limitations and potential for further development We here have presented a new and innovative approach to simulate wind disturbance in forest ecosystems at individual tree resolution, with the location, area and cohort damaged by wind being derived dynamically by means of process-based simulations. Model behavior was found well in line with general expectations

and previous sensitivity analyses with process-based wind risk models (e.g., Gardiner et al., 2000; Byrne, 2011). Simulations for the storm Gudrun at the Växjö study landscape showed that simulated wind damage corresponds well with satellite-based observations (Table 2), and documents the utility of our approach. However, several limitations remain in the current model which could be addressed in future improvements of the approach. An acclimation effect of trees to wind, for instance, is currently not included in the model. Yet previous empirical work has shown that trees frequently exposed to strong winds adapt their allocation regime and improve their stability, e.g., by increased root growth (Cucchi et al., 2004; Danjon et al., 2005; Nicoll et al., 2008; Reubens et al., 2009). Not accounting for these aspects of acclimation is a likely reason for the overestimated wind damage along the western edge of the Växjö landscape (Fig. 3), where the study landscape borders agricultural land and a lake, respectively. A related issue is that we currently focus the detailed wind damage calculations on stand edges. This simplification sharply decreases the number of simulated cells for which these calculations have to be performed and drastically decreases computation time of the model (which is of importance e.g., in scenario analyses requiring long time horizons and a large number of replicates). This

8

R. Seidl et al. / Environmental Modelling & Software 51 (2014) 1e11

Table 2 Comparison of simulation results to the observed wind damage for the storm Gudrun at the Växjö landscape. S00, S10, and S11 are simulation experiments with identical wind forcing but different levels of within-stand heterogeneity e see Fig. 1. Values are averages over 10 replicated simulations. Asterisks indicate statistical significance of Kappa values at a ¼ 0.05.

Percent of landscape damaged (%)a Deviation to observed damage (%) per hectare (mean  sd)a Odds ratio (dimensionless) predicted/observeda Accuracy of predicted damage (%)b Accuracy of predicted severe (>50%) damage (%)b Conditional Kappa for damaged areab Conditional Kappa for severely (>50%) damaged areab a b

S00 No structural variation No spatial variation

S10 Full structural variation Random spatial variation

S11 Full structural variation Full spatial variation

14.0 6.4  34.1 0.587 36.5 17.2 0.076* 0.094*

16.7 4.2  34.2 0.721 44.1 21.7 0.119* 0.125*

21.1 0.5  34.3 0.963 60.4 19.4 0.151* 0.096*

10 m horizontal resolution. 100 m horizontal resolution.

focus on edges is supported by the success of previous risk models applying a similar approach (Blennow and Sallnäs, 2004; Blennow et al., 2010). In contrast to these previous approaches, however, we here define edges as changes in canopy rugosity at a horizontal resolution of 10 m cells, and thus also conduct wind loading calculations for, e.g., individual seed trees or for tall trees (>10 m taller than the surrounding canopy) in vertically structured stands. In general, since fetch is dynamically calculated (Eq. (8)), the approach presented here would be applicable also to calculate wind loading for closed stand conditions. While detailed airflow models exist to calculate the effect of terrain and local obstacles on wind loading (e.g., Mortensen et al., 1998), we have chosen a simplified statistical approach to derive the differences in landscape-scale wind exposure for the current case study. The utility of the TTD for this purpose was recently documented for instance by the analysis of Albrecht et al. (2010). In general, considering the sensitivity of the model to wind speed (Table 3) the documented ability to identify hotspots of wind damage for Gudrun with satisfactorily fidelity (Fig. 4) underlines the utility of our statistical approach. However, it has to be noted that statistical relationships as the one used here to derive landscape-scale variation in windiness are specific for a certain area

and their generality and transferability are likely low. While the wind modeling approach presented in Section 2.3 could also be driven by results of detailed airflow modeling, previous studies found no clear advantage of such more complex approaches compared to statistical methods relying primarily on topographic exposure (Suarez et al., 1999). In addition, the empirical approach chosen here has the advantage of high computational efficiency. The resulting good overall runtime performance of the model (wind simulations for the Vaxjö landscape take only a few seconds) supports the applicability of the approach for large landscapes and in replicated scenario analyses. With regard to replicates, a moderate number of replicated simulations yielded robust results in this study (cf. Fig. 5), where only the within-event wind direction and speed were allowed to vary stochastically around a deterministically set mean value (see Table 1, step 1). 4.2. Implications and applications We here have shown that tree-level heterogeneity is important for modeling the impact of wind on forest ecosystems. Our results with a novel process-based wind model show that neglecting structural and spatial heterogeneity in forest stands will likely lead

Fig. 4. Maps of observed and simulated wind damage percentage (% of 10 m grid cells damaged per hectare) at the Växjö landscape. Simulations relate to the experiment accounting for the full structural and spatial variation (S11). An overlay of 20 m elevation isolines provides additional information on the topography of the study landscape. Simulation results are averages over 10 replicates.

R. Seidl et al. / Environmental Modelling & Software 51 (2014) 1e11

<0.5 0.5−1

1−2

2−5

5−10

0.4 0.3 0.2 0.1

proportion of height class damaged

(b)

0.0 0.1 0.2 0.3 0.4 0.5 0.6

proportion of damaged patches

(a)

9

>10

10

15

patch size (ha)

20

25

dominant height (m)

Fig. 5. (a) Patch size distribution of areas damaged by wind. (b) Damage proportion over dominant height classes. Observations are indicated in dark gray, simulations in light gray. Simulation results are given for the S11 initialization variant (see Fig. 2), and represent averages over 10 replicates, with error bars indicating the extreme values over all replicates.

Wind speed Event duration Wind direction Soil frozen

Change in driver (%)

Change in disturbed area (%)

Change in average patch size (%)

20 þ20 20 þ20 20 þ20

41.2 þ53.3 2.4 þ0.8 2.0 1.1 17.1

39.8 þ80.4 4.2 þ2.5 3.2 þ0.6 19.2

0.20 0.10 0.00

Table 3 Sensitivity of the Växjö landscape to storm characteristics. Key drivers were changed individually relative to the values observed for the storm Gudrun (described in detail in Section 2.3). Simulations were conducted with the S11 initialization variant (see Fig. 2), and results are averages over 10 replicates.

model is not only able to estimate effects on indicators such as ecosystem timber and C stocks (e.g., Seidl et al., 2012b), but can also track the legacies of such disturbances dynamically over time. In this regard it has to be noted that disturbance impacts are an emergent property of the simulated vegetation dynamics in iLand. While many current forest landscape simulators impose wind damage of a certain size and severity (e.g., Scheller et al., 2011), simulating these processes from the dynamic interactions between vegetation and wind is likely to improve the robustness of impact assessments in the context of changing environmental conditions. With regard to robustness it is important to note that the current approach e as is the case for most other process-based wind risk models e is highly sensitive to the wind data driving the model. Based on current process understanding, even small changes in the extreme wind climate can have big impacts on forest ecosystems. This underlines the importance of improved data on peak wind speeds, both for the observation period as well as for future projections. While future changes in wind are still highly uncertain, some studies suggest an increased frequency and intensity of strong winds (e.g., Nikulin et al., 2011; Pryor et al., 2012), which would have the potential to further intensify the forest disturbance regime. Another important factor in the context of climate change and wind disturbance is soil frost: Most storm events in Europe occur in winter, and frozen soils have been found to considerably increase the wind stability of trees (Silins et al., 2000; Gardiner et al., 2010). If a warming climate leads to decreasing periods of

proportion of damage

to underestimated wind damage. This finding is particularly noteworthy considering that our study landscape consisted of forest under an even-aged management system (but see Suzuki et al. (2012) on the variation of spatial structure within stands of the same age and development stage). However, our result of higher simulated wind damage in more heterogenous stand conditions does not allow inferences on the effect of different management regimes on wind stability, as all initialization variants studied here represent generic, homogenized renderings of the same management system. The current approach is, however, not limited in this regard, and could for instance also be applied to study the effects of changing the management regime from even-aged to uneven-aged forestry in future applications (see e.g., Seidl et al., 2008). This presents a major improvement over many previous wind risk models, which implicitly assume even-aged and homogenous stands (see the discussion by Gardiner et al. (2008)). Furthermore, we have integrated the tree-level wind risk model into a dynamic landscape modeling framework. This integration allows progress beyond the assessment of the wind risk to forests at a single point in time, and enables the simulation of dynamic spatio-temporal feedbacks via changes in stand structure and composition (i.e., a wind disturbance event can change the structure and composition of a forest landscape, which in turn will lead to a modified wind risk in subsequent years). Moreover, dynamic interactions between disturbance agents (e.g., wind and bark beetles) could be included in this framework in the future (see Seidl et al., 2011c). An additional asset of implementing the presented wind risk model into the iLand simulation framework is the ability to quantify the impacts of disturbances on ecosystem services. The

0

20

40

60

80

100

iterations Fig. 6. Temporal evolution of simulated wind damage during the Gudrun event. The figure shows the relative contribution of every iteration to total simulated damage for each of the ten replicates (gray) and their average (bold black).

10

R. Seidl et al. / Environmental Modelling & Software 51 (2014) 1e11

soil frost and changed soil moisture regimes in winter, an increased susceptibility to wind has to be expected, a trend that has already been empirically documented for past decades (Usbeck et al., 2010). 5. Conclusions Disturbance regimes are intensifying in Europe, and managing disturbances is thus an increasingly important issue for ecosystem management. In order to support management decision making in this regard models are needed that are robust (founded in processunderstanding) yet flexible (with regard to the incorporation of different management strategies). We here have presented a process-based simulation model of wind disturbance, accounting for tree-level heterogeneity by means of an individual-based forest landscape model. The presented approach accounts for local sheltering by neighbors, and dynamically simulates the progression of wind damage during a wind event, with disturbance size and severity an emerging property of the simulation. Analyses with this novel model showed that neglecting spatial and structural withinstand heterogeneity has a considerable impact on simulated wind damage, and thus supports the notion that complexity matters for understanding and modeling disturbances (Turner et al., 2013). Our model-based analyses furthermore highlight a considerable sensitivity of disturbance impacts to changes in extreme wind climate as well as soil frost. The newly presented model can be used to study the effect of different silvicultural systems on future wind risk as well as to quantify the impacts of disturbance on ecosystem services, and thus contributes to building the capacity to address changing disturbance regimes in ecosystem management. Acknowledgments We thank S. Hale, Forest Research, UK, for providing data on tree size, competition, and turning moment coefficients for reanalysis, the Department of Forest Resource Management, Swedish University of Agricultural Sciences, for providing kNN data on forest cover, M. Andersson, Swedish University of Agricultural Sciences, for preparing the tree cover data, and the Swedish Forest Agency for providing satellite data on wind damage after the storm “Gudrun”. We furthermore thank B. Gardiner, INRA, France for helpful comments on the intricacies of mechanistic wind modeling. Two anonymous reviewers are acknowledged for providing detailed and insightful comments that helped to improve the manuscript further. This work was partly funded by two European Community’s Seventh Framework Program Marie Curie Fellowship to RS (grant agreements 237085 and 334104). KB acknowledges funding from the research program Models for Adaptive Forest Management (MOTIVE) within the European Community’s Seventh Framework Program (project no. 226544). References Albrecht, A., Hanewinkel, M., Bauhus, J., Kohnle, U., 2010. How does silviculture affect storm damage in forests of south-western Germany? Results from empirical modeling based on long-term observations. Eur. J. For. Res. 131, 229e 247. ASTER, 2011. Global Digital Elevation Map. Advanced Spaceborn Thermal Emission and Reflection Radiometer. http://asterweb.jpl.nasa.gov/gdem.asp (accessed 10.12.12.). Beniston, M., Stephenson, D.B., Christensen, O.B., Ferro, C.A.T., Frei, C., Goyette, S., Halsnaes, K., Holt, T., Jylhä, K., Koffi, B., Palutikof, J., Schöll, R., Semmler, T., Woth, K., 2007. Future extreme events in European climate: an exploration of regional climate model projections. Clim. Change 81, 71e95. Bennett, N.D., Croke, B.F.W., Guariso, G., Guillaume, J.H.A., Hamilton, S.H., Jakeman, A.J., Marsili-Libelli, S., Newham, L.T.H., Norton, J.P., Perrin, C., Pierce, S.A., Robson, B., Seppelt, R., Voinov, A.A., Fath, B.D., Andreassian, V., 2013. Characterizing performance of environmental models. Environ. Model. Softw. 40, 1e20.

Bivand, R.S., Pebesma, E.J., Gomez-Rubio, V., 2008. Applied Spatial Data Analysis with R. Springer, New York, p. 374. Blennow, K., Andersson, M., Sallnäs, O., Olofsson, E., 2010. Climate change and the probability of wind damage in two Swedish forests. For. Ecol. Manage. 259, 818e830. Blennow, K., Sallnäs, O., 2004. WINDAda system of models for assessing the probability of wind damage to forest stands within a landscape. Ecol. Model. 175, 87e99. Byrne, K.E., 2011. Mechanistic Modeling of Windthrow in Spatially Complex Mixed Species Stands in British Columbia (PhD thesis). University of British Columbia, 185 pp. Byrne, K.E., Mitchell, S.J., 2013. Testing of WindFIRM/ForestGALES_BC: a hybridmechanistic model for predicting windthrow in partially harvested stands. Forestry 86, 185e199. Cucchi, V., Meredieu, C., Stokes, A., Berthier, S., Bert, D., Najar, M., Denis, A., Lastennet, R., 2004. Root anchorage of inner and edge trees in stands of Maritime pine (Pinus pinaster Ait.) growing in different podzolic soil conditions. Trees 18, 460e466. Danjon, F., Fourcaud, T., Bert, D., 2005. Root architecture and wind-firmness of mature Pinus pinaster. New Phytol. 168, 387e400. Duursma, R.A., Marshall, J.D., Robinson, A.P., Pangle, R.E., 2007. Description and test of a simple process-based model of forest growth for mixed-species stands. Ecol. Model. 203, 297e311. Gardiner, B., Blennow, K., Carnus, J.M., Fleischer, P., Ingemarson, F., Landmann, G., Lindner, M., Marzano, M., Nicoll, B., Orazio, C., Peyron, J.L., Reviron, M.P., Schelhaas, M.J., Schuck, A., Spielmann, M., Usbeck, T., 2010. Destructive Storms in European Forests: Past and Forthcoming Impacts. Final Report to European CommissioneDG Environment. EFIATLANTIC, European Forest Institute, 138 pp. Gardiner, B., Byrne, K., Hale, S., Kaminura, K., Mitchell, S.J., Peltola, H., Ruel, J.C., 2008. A review of mechanistic modelling of wind damage risk to forests. Forestry 81, 447e463. Gardiner, B., Peltola, H., Kellomäki, S., 2000. Comparison of two models for predicting the critical wind speeds required to damage coniferous trees. Ecol. Model. 129, 1e23. Gardiner, B.A., Stacey, G.R., Belcher, R.E., Wood, C.J., 1997. Field and wind tunnel assessments of the implications of respacing and thinning for tree stability. Forestry 70, 233e252. Hale, S.E., Gardiner, B.A., Wellpott, A., Nicoll, B.C., Achim, A., 2012. Wind loading of trees: influence of tree size and competition. Eur. J. For. Res. 131, 203e217. Hannah, P., Palutikof, J.P., Quine, C.P., 1995. Predicting windspeeds for forest areas in complex terrain. In: Coutts, M.P., Grace, J. (Eds.), Wind and Trees. University of Edinburgh, pp. 113e130. Hegyi, F., 1974. A simulation model for managing jack pine stands. In: Fries, J. (Ed.), Growth Models for Tree and Stand Simulation. Royal College of Forestry, Stockholm, Sweden, pp. 74e90. Karlsson, C., Westman, S.-E., 1987. Skogsuppskattning skogsinventering. Liber Förlag, Stockholm, ISBN 91-40-11916-5. Kätterer, T., Andren, O., 2001. The ICBM family of analytically solved models of soil carbon, nitrogen and microbial biomass dynamicsddescriptions and application examples. Ecol. Model. 136, 191e207. Klawa, M., Ulbrich, U., 2003. A model for the estimation of storm losses and the identification of severe winter storms in Germany. Nat. Hazards Earth Sys. Sci. 3, 725e732. Kolström, M., Lindner, M., Vilén, T., Maroschek, M., Seidl, R., Lexer, M.J., Netherer, S., Kremer, A., Delzon, S., Barbati, A., Marchetti, M., Corona, P., 2011. Reviewing the science and implementation of climate change adaptation measures in European forestry. Forests 2, 961e982. Lagergren, F., Jönsson, A.M., Blennow, K., Smith, B., 2012. Implementing storm damage in a dynamic vegetation model for regional applications in Sweden. Ecol. Model. 247, 71e82. Landsberg, J.J., Waring, R.H., 1997. A generalised model of forest productivity using simplified concepts of radiation-use efficiency carbon balance and partitioning. For. Ecol. Manage. 95, 209e228. Larsén, X.G., Mann, J., 2006. The effects of disjunct sampling and averaging time on maximum mean wind speeds. J. Wind Eng. Ind. Aerodyn. 94, 581e602. Lindroth, A., Lagergren, F., Grelle, A., Klemedtsson, L., Langvall, O., Weslien, P., Tuulik, J., 2009. Storms can cause Europe-wide reduction in forest carbon sink. Glob. Change Biol. 15, 346e355. Mortensen, N.G., Landberg, L., Troen, I., Petersen, E.L., 1998. Wind Atlas Analysis and Application Program (WASP). Risø National Laboratory, Roskilde, Denmark. Morzuch, B.J., Ruark, G.A., 1991. Principle components regression to mitigate the effects of multicollinearity. For. Sci. 37, 191e199. Nicoll, B.C., Gardiner, B.A., Rayner, B., Peace, A.J., 2006. Anchorage of coniferous trees in relation to species, soil type and rooting depth. Can. J. For. Res. 36, 1871e 1883. Nicoll, B.C., Gardiner, B.A., Peace, A.J., 2008. Improvements in anchorage provided by the acclimation of forest trees to wind stress. Forestry 81, 389e398. Nikulin, G., Kjellström, E., Hansson, U., Strandberg, G., Ullerstig, A., 2011. Evaluation and future projections of temperature, precipitation and wind extremes over Europe in an ensemble of regional climate simulations. Tellus 63A, 41e55. Nilsson, C., Goyette, S., Bärring, L., 2007. Relating forest damage data to the wind field from high-resolution RCM simulations: case study of Anatol striking Sweden in December 1999. Glob. Planet. Change 57, 161e176. Nitschke, C.R., Innes, J.L., 2008. A tree and climate assessment tool for modelling ecosystem response to climate change. Ecol. Model. 210, 263e277.

R. Seidl et al. / Environmental Modelling & Software 51 (2014) 1e11 Papaik, M.J., Canham, C.D., 2006. Species resistance and community response to wind disturbance regimes in northern temperate forests. J. Ecol. 94, 1011e1026. Paul, K.I., Polglase, P.J., Smethurst, P.J., O’Connell, A.M., Carlyle, C.J., Khanna, P.K., 2004. Soil temperature under forests: a simple model for predicting soil temperature under a range of forest types. Agric. For. Meteorol. 121, 167e182. Peltola, H., Kellomäki, S., Hassinen, A., Granander, M., 2000. Mechanical stability of Scots pine, Norway spruce and birch: an analysis of tree-pulling experiments in Finland. For. Ecol. Manage. 135, 143e153. Peltola, H., Kellomäki, S., Väisänen, H., Ikonen, V.P., 1999. A mechanistic model for assessing the risk of wind and snow damage to single trees and stands of Scots pine, Norway spruce, and birch. Can. J. For. Res. 29, 647e661. Pryor, S.C., Barthelmie, R.J., Clausen, N.E., Drews, M., MacKellar, N., Kjellström, E., 2012. Analyses of possible changes in intense and extreme wind speeds over northern Europe under climate change scenarios. Clim. Dyn. 38, 189e208. Puettmann, K.J., Coates, K.D., Messier, C., 2009. A Critique of Silviculture. Managing for Complexity. Island Press, Washington, DC, p. 206. Quine, C.P., White, I.M.S., 1998. The potential of distance-limited topex in the prediction of site windiness. Forestry 71, 325e332. R Development Core Team, 2011. R: a Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria, ISBN 3900051-07-0. http://www.R-project.org/ (accessed 10.12.12.). Raupach, M.R., 1992. Drag and drag partition on rough surfaces. Bound. Layer Meteorol. 60, 375e395. Raupach, M.R., 1994. Simplified expressions for vegetation roughness length and zero-plane displacement as functions of canopy height and area index. Bound. Layer Meteorol. 71, 211e216. Reese, H., Nilsson, M., Pahlén, T.G., Hagner, O., Joyce, S., Tingelöf, U., Egberth, M., Olsson, H., 2003. Countrywide estimates of forest variables using satellite data and field data from the National Forest Inventory. Ambio 32, 542e548. Reubens, B., Pannemans, B., Danjon, F., de Proft, M., de Baets, S., de Baerdemaeker, J., Poesen, J., Muys, B., 2009. The effect of mechanical stimulation on root and shoot development of young containerised Quercus robur and Robinia pseudoacacia trees. Trees 23, 1213e1228. Ruel, J.C., Mitchell, S.J., Dornier, M., 2002. A GIS based approach to map wind exposure for windthrow hazard rating. North. J. Appl. For. 19, 183e187. Schelhaas, M.J., Hengeveld, G., Moriondo, M., Reinds, G.J., Kundzewicz, Z.W., ter Maat, H., Bindi, M., 2010. Assessing risk and adaptation options to fires and windstorms in European forestry. Mitig. Adapt. Strateg. Glob. Change 15, 681e 701. Schelhaas, M.J., 2008. The wind stability of different silvicultural systems for Douglas-fir in The Netherlands: a model-based approach. Forestry 81, 399e414. Schelhaas, M.J., Kramer, K., Peltola, H., van der Werf, D.C., Wijdeven, S.M.J., 2007. Introducing tree interactions in wind damage simulation. Ecol. Model. 207, 197e209. Schelhaas, M.J., Nabuurs, G.-J., Schuck, A., 2003. Natural disturbances in the European forests in the 19th and 20th centuries. Glob. Change Biol. 9, 1620e1633. Scheller, R.M., Hua, D., Bolstad, P.V., Birdsey, R.A., Mladenoff, D.J., 2011. The effects of forest harvest intensity in combination with wind disturbance on carbon dynamics in Lake States Mesic Forests. Ecol. Model. 222, 144e153. Schiesser, H.H., Pfister, C., Bader, J., 1997. Winter storms in Switzerland north of the Alps 1864/1865e1993/1994. Theor. Appl. Climatol. 58, 1e19. Seidl, R., Fernandes, P.M., Fonseca, T.F., Gillet, F., Jönsson, A.M., Merganicova, K., Netherer, S., Arpaci, A., Bontemps, J.D., Bugmann, H., Gonzalez-Olabarria, J.R.,

11

Lasch, P., Meredieu, C., Moreira, F., Schelhaas, M.J., Mohren, F., 2011c. Modelling natural disturbances in forest ecosystems: a review. Ecol. Model. 222, 903e924. Seidl, R., Schelhaas, M.J., Lexer, M.J., 2011a. Unraveling the drivers of intensifying forest disturbance regimes in Europe. Glob. Change Biol. 17, 2842e2852. Seidl, R., Blennow, B., 2012. Pervasive growth reduction in Norway spruce forests following wind disturbance. PLoS ONE 7 (3), e33301. Seidl, R., Rammer, W., Jäger, D., Lexer, M.J., 2008. Impact of bark beetle (Ips typographus L.) disturbance on timber production and carbon sequestration in different management strategies under climate change. For. Ecol. Manage. 256, 209e220. Seidl, R., Rammer, W., Lexer, M.J., 2011b. Adaptation options to reduce climate change vulnerability of sustainable forest management in the Austrian Alps. Can. J. For. Res. 41, 694e706. Seidl, R., Rammer, W., Scheller, R.M., Spies, T.A., 2012a. An individual-based process model to simulate landscape-scale forest ecosystem dynamics. Ecol. Model. 231, 87e100. Seidl, R., Spies, T.A., Rammer, W., Steel, E.A., Pabst, R.J., Olsen, K., 2012b. Multi-scale drivers of spatial variation in old-growth forest carbon density disentangled with Lidar and an individual-based landscape model. Ecosystems 15, 1321e1335. Silins, U., Lieffers, V.J., Bach, L., 2000. The effect of temperature on mechanical properties of standing lodgepole pine trees. Trees 14, 424e428. Skogsstyrelsen, 2006a. Stormen 2005 e en skoglig analys. Meddelande 1, 2006. Jönköping, Sweden. ISSN 1100e0295. Skogsstyrelsen, 2006b. Kartering av stormfälld skog som behöver återbeskogas. Rapport, 14 pp. SMHI, 2012. Observationsdata för 52 synoptiska stationer, 1961e2010. Sveriges meteorologiska och hydrologiska institute. http://www.smhi.se/klimatdata/ meteorologi/vind (accessed 10.12.12.). SMHI, 2005. Januaristormen 2005. Sveriges Meteorologiska och Hydrologiska Institut. Faktablad 25, November 2005. Suarez, J.C., Gardiner, B.A., Quine, C.P., 1999. A comparison of three methods for predicting wind speeds in complex forested terrain. Meteorol. Appl. 6, 329e 342. Suzuki, S.N., Kachi, N., Suzuki, J.I., 2012. Variability of local spatial structure in a wave-regenerated Abies forest. Ecol. Res. 27, 893e901. Turner, M.G., Donato, D.C., Romme, W.H., 2013. Consequences of spatial heterogeneity for ecosystem services in changing forest landscapes: priorities for future research. Landsc. Ecol. 28, 1081e1097. Usbeck, T., Wohlgemuth, T., Dobbertin, M., Pfister, C., Bürgi, A., Rebetez, M., 2010. Increasing storm damage to forests in Switzerland from 1858 to 2007. Agric. For. Meteorol. 150, 47e55. Valinger, E., Fridman, J., 2011. Factors affecting the probability of windthrow at stand level as a result of Gudrun winter storm in southern Sweden. For. Ecol. Manage. 262, 398e403. VanDerWal, J., Falconi, L., Januchowski, S., Shoo, L., Storlie, C., 2011. SDMTools: Species Distribution Modelling Tools: Tools for Processing Data Associated with Species Distribution Modelling Exercises. R Package Version 1.1-5 http://CRAN. R-project.org/package¼SDMTools (accessed 10.12.12.). Wu, H., Sharp, P.J., Walker, J., Penridge, L.K., 1985. Ecological field theory: a spatial analysis of resource interference among plants. Ecol. Model. 29, 215e243. Zeng, H., Peltola, H., Talkkari, A., Venäläinen, A., Strandman, H., Kellomäki, S., Wang, K., 2004. Influence of clear-cutting on the risk of wind damage at forest edges. For. Ecol. Manage. 203, 77e88.