Simulating damage for wind storms in the land surface model ORCHIDEE-CAN (revision 4262)

Earth System Models (ESMs) are currently the most advanced tools with which to study the interactions between humans, ecosystem productivity and the climate. The inclusion of storm damage in ESMs has long been hampered by their big-leaf approach which ignores the canopy structure information that is required for process-based wind throw modelling. Recently the big-leaf assumptions in the large scale land surface model ORCHIDEE-CAN were replaced by a three dimensional description of the canopy structure. This opened the way to the integration of the processes from the small-scale wind damage 5 risk model ForestGALES into ORCHIDEE-CAN. The integration of ForestGALES into ORCHIDEE-CAN required, however, developing numerically cheap solutions to deal with: (1) landscape heterogeneity, i.e., recent forest edge to distinguish the inner and outer area around the gaps for different wind gust parametrisation; (2) downscaling spatially and temporally aggregated wind fields to obtain more realistic wind speeds that would represents gusts; and (3) downscaling storm damage within the 2500 km pixels of ORCHIDEE-CAN. This new version of ORCHIDEE-CAN was parametrised over Sweden. Subsequently, 10 the performance of the model was tested against data for historical storms in Southern Sweden between 1951 and 2010, and South-western France in 2009. In years without big storms, here defined as a storm damaging less than 15×10m of wood in Sweden, the model error is 1.62×10m which is about 100 % of the observed damage. For years with big storms, such as Gudrun in 2005, the model error increased to 5.05×10m which is between 10 % and 50 % of the observed damage. When the same model parameters were used over France, the model reproduced a decrease in leaf area index and an increase 15 in albedo, in accordance with SPOT-VGT and MODIS records following the passing of Cyclone Klaus in 2009. The current version of ORCHIDEE-CAN (revision 4262) is therefore expected to have the capability to capture the dynamics of forest


Introduction
During the last 15 years, Western Europe has been severely affected by storms with the six most damaging European storms ever recorded hitting France (Lothar and Martin, December 1999; Klaus in January 2009), Sweden (Gudrun, January 2005), and Germany (Lothar, December 1999;Kyrill, January 2007) as well as neighbouring countries (the Netherlands, Belgium, Switzerland, Czech Republic, Slovakia, and the Baltic States) (Gardiner et al., 2010).The short-term impact of wind on forests includes billions of Euro in damage to wood stocks, loss of valuable protected old forest stands, increased fire occurrence (Miller et al., 2007), pest vulnerability (Komonen et al., 2011), and a temporary decrease in the productivity of the remaining forest stands (Merrens and Peart, 1992;Everham et al., 1996;Seidl and Blennow, 2012).Wind throw is reported to be the cause of 57 % of forest disturbances in Europe and is thus more significant than stand-replacing disturbances through pest attack (26 %) and through fire (16 %) (Schelhaas et al., 2003;Seidl et al., 2014).Furthermore, a global literature review indicates that wind disturbance triggers fire activity in warm and dry climates but induces pathogens and/or insect outbreak in warm and wet climates (Seidl et al., 2017).
Similar to fire (Randerson et al., 2006), the direct and indirect effects of wind disturbances may contribute to the top of the atmosphere radiative forcing (O'Halloran et al., 2012).The direct effects, such as a reduction in leaf area index (Juárez et al., 2008), transpiration (Negrón-Juárez et al., 2010), an increase in the surface albedo (Planque et al., 2017), and decrease in roughness (Zhu, 2008), have been shown to impact the regional climate, i.e. following the storm Klaus in 2009, cloud cover frequency was observed to decrease over Les Landes in south-western France (Teuling et al., 2017).The indirect effects include a reduction of the gross productivity by damage to the rooting system, increased tree mortality due to facilitating insect or pathogen outbreaks (Sturrock et al., 2011), or a change in the forest ecosystem by shifting canopy structure (Lin et al., 2017).
Storm-induced disturbances are likely to provide feedbacks on climate through direct effects such as increasing greenhouse gas emissions (Lindroth et al., 2009), increasing surface albedo (Planque et al., 2017), and decreasing local cloud frequency (Teuling et al., 2017), as well as indirect effects such as increased natural disturbances, a reduced logging rate in subsequent years, increased weathering, and increased C, N, and cation leaching (Futter et al., 2011;Köhler et al., 2011).Increased weathering and leaching could even extend the effects of wind throw from the land to the oceans since terrestrial processes have been found to play an important role in the lateral C fluxes to the oceans through inland waters (Battin et al., 2009;Regnier et al., 2013).
ESMs can be seen as a mathematical representation of the major biophysical processes of the natural world (Sellers et al., 1986;Henderson-Sellers et al., 1996;Sellers, 1997;Bonan, 2008) and are currently the most advanced tools to study the interactions among humans, their use of vegetated ecosystems, and the climate (Jackson et al., 2005;Swann et al., 2012;Naudts et al., 2016).Although Earth system modelling groups dedicate considerable resources to studying the effects of fires (Lasslop et al., 2014), forest management (Naudts et al., 2016), land cover changes (Swann et al., 2012;Devaraju et al., 2015), and shifting cultivation (Wilkenskjeld et al., 2014;Yue et al., 2018), storm-induced disturbances and their climate feedback are not yet explicitly dealt with in ESMs.The objective of this study is to develop the model capability for the ESM IPSL-CM through its land component ORCHIDEE-CAN to simulate the effects of wind storms on the land surface by building on a good understanding of ecosystem-level processes (Hale et al., 2012(Hale et al., , 2015)).Until the direct and indirect climate effects of wind storms have been quantified, implementing storm damage in ESMs is justified by its precursory effect on other natural disturbances such as fires and insects (Seidl et al., 2014).Abrupt mortality from drought, wind storms, fires, pests, and their interactions will need to be accounted for if ESMs are to be used to quantify the effects of future climate on forest dynamics and forest resilience.
Several classifications mostly based on maximum wind speeds within a cyclone have been proposed to define different types of storms ranging from depressions to hurricanes.In general, storm damage strongly depends on the frequency and intensity of storms.For example, when gusts within a cyclone exceed 20 m s −1 , uprooting and stem breakage is to be expected, resulting in severe damage.When wind speeds remain below 17.1 m s −1 the damage is expected to be less devastating but nevertheless substantial (Scatena et al., 2004).The implemented approach, however, did not require a classification of storms as it can simulate the transition from no storm damage at low wind speeds to storm damage at high wind speeds (Sect.2.2.4).

Model description and parameterization
ORCHIDEE-CAN revision 2566 (see Sect. 2.1) was further developed by implementing the modifications and additions listed below (see Sect. 2.2 to 2.3), resulting in revision 4262.From revision 4262 onwards, ORCHIDEE-CAN has the capability to simulate tree mortality from wind storms.The notation used to describe the model is listed in full in Table 1.

ORCHIDEE-CAN (revision 2566)
ORCHIDEE (Krinner et al., 2005;Bellassen et al., 2010) is the land surface model of the IPSL (Institute Pierre Simon Laplace) Earth system model.Hence, by conception, it can be coupled to a global circulation model.In a coupled setup, the atmospheric conditions affect the land surface and the land surface, in turn, affects the atmospheric conditions.
Geosci.Model Dev., 11, 771-791, 2018 www.geosci-model-dev.net/11/771/2018/However, when a study focuses on changes in the land surface rather than on the interaction with climate, it also can be run off-line as a stand-alone land surface model.The standalone configuration receives atmospheric conditions such as temperature, humidity, and wind, to mention a few, from the so-called "forcing files".Unlike the coupled set-up, which needs to run on the global scale, the stand-alone configuration can cover any area ranging from the global domain to a single grid point.Although ORCHIDEE does not enforce a spatial or temporal resolution, the model does use a spatial grid and equidistant time steps.The spatial resolution is an implicit user setting that is determined by the coarsest resolution of the forcing data and the boundary conditions, i.e. the vegetation distribution, climatological forcing data, and the soil map.If higher-resolution drivers are available the model can then be run on that scale.If site-level drivers are available then simulations on the site scale are feasible.ORCHIDEE can run on any temporal resolution; however, this apparent flexibility is rather restricted as the processes are formalized at given time steps: half-hourly (i.e.photosynthesis and energy budget), daily (i.e.net primary production), and annual (i.e.vegetation dynamics).Hence, meaningful simulations have a temporal resolution of 15 min to 1 h for the energy balance, water balance, and photosynthesis calculations.
ORCHIDEE builds on the concept of meta-classes to describe vegetation distribution.By default it distinguishes 13 such meta-classes (one for bare soil, eight for forests, two for grasslands, and two for crop lands).Each meta-class can be subdivided into an unlimited number of plant functional types (PFTs).When simulations make use of species-specific parameters and age classes (as is the case in this study), sev- eral PFTs belonging to a single meta-class will be defined.Biogeochemical and biophysical variables are calculated for each PFT.ORCHIDEE-CAN (revision 2566) (Naudts et al., 2015;Ryder et al., 2016;Chen et al., 2016;McGrath et al., 2016) is one of the branches for the ORCHIDEE model development, which was selected to simulate large-scale wind-throw and storm damage because, contrary to most land surface models, ORCHIDEE-CAN simulates dynamic canopy structures, a feature essential to simulate the likelihood of wind throw and the subsequent damage.Changes in canopy structure resulting from wind throw are then accounted for in the calculations of the carbon, water, and energy exchange between the land surface and the lower atmosphere.
In ORCHIDEE-CAN, tree height and crown diameter are linked to tree diameter through allometric relationships.Individual tree canopies are simulated as spherical elements with their horizontal location following a Poisson distribution across the stand (Naudts et al., 2015).A forest is represented by a user-defined number of diameter classes.Each diameter class represents trees with a different mean diameter and height and therefore informs the user about the social position of trees within the canopy.The difference in social position within a stand is the basis of intra-stand competition, which accounts for the fact that trees with a dominant position in the canopy are more likely to intercept light than suppressed trees and therefore contribute more to the stand-level photosynthesis and biomass growth (Deleuze et al., 2004).The allocation scheme is based on the pipe model theory (Shinozaki et al., 1964) and its implementation by Sitch et al. (2003).The scheme allocates carbon to different biomass pools (leaves, fine roots, and sapwood) while respecting the differences in longevity and hydraulic conductivity between the pools (Naudts et al., 2015).
At the start of a simulation, each PFT contains a userdefined number of diameter classes.This number is held constant, whereas the boundaries of the classes are adjusted throughout a simulation to accommodate temporal evolution in the stand structure.By using flexible class boundaries with a fixed number of diameter classes, different forest structure can be simulated.An even-aged forest, for example, is simulated with a small diameter range between the smallest and largest trees.All trees thus belong to the same stratum.An uneven-aged forest is simulated by applying a large range between the diameter classes.Different diameter classes will therefore represent different strata.Each diameter class contains a single modelled tree.The modelled tree is replicated to give realistic stand densities.Following this, tree growth, canopy dimensions, and stand density are updated.Throughout a simulation, individual tree mortality causes stand density to decrease.In ORCHIDEE-CAN individual tree mortality is caused by self-thinning and forest management.In the absence of these processes, a constant rate of so-called environmental background mortality is applied.The inclusion of the so-called environmental background mortality im-plicitly accounts for mortality through fires, pests, and wind throw.Following the development of the wind-throw and storm damage module in revision 4262, mortality from wind throw is now explicitly accounted for and thus no longer included in the so-called environmental background mortality.
Furthermore, age classes are used after land cover change and forest management events to explicitly simulate the regrowth of the forest.Following a land cover change, biomass and soil carbon pools (but not soil water columns) are either merged or split to represent the various outcomes of a land cover change.This dynamic approach to stand and landscape structure is exploited in other parts of the model, i.e. precipitation interception, transpiration, energy budget calculations, the radiation scheme, absorbed light for photosynthesis, and, since revision 4262, tree mortality through wind storms.

ORCHIDEE-CAN (revision 4262)
In ORCHIDEE-CAN (revision 4262) the biomass of the different pools, leaf area index, crown volume, crown density, stem diameter, stem height, and stand density are simulated as the accumulated growth and are passed to the wind-throw module.The wind-throw module calculates the critical wind speed based on the principles applied in ForestGALES (Gardiner et al., 2000) and storm damage based on the approach developed and tested by Anyomi et al. (2017).Figure 1 summarizes the major components of storm damage calculations in ORCHIDEE-CAN.A more detailed description of the different components in the figure is presented in following sections.

Critical wind speeds (ForestGALES)
The presence or absence of storm damage in a forest stand can be modelled with the concept of critical wind speed.If the wind speed exceeds the critical wind speed of a forest, the forces applied by the wind speed are sufficient to overturn the whole tree or break its stem.The exact value of the critical wind speed depends on the canopy structure (Vollsinger et al., 2005;Hale et al., 2012), the tree species, the soil properties, and the root profiles (Nicoll et al., 2006).In this study, the physics formalized in ForestGALES (Gardiner et al., 2010;Hale et al., 2015), a hybrid mechanistic forest wind damage risk model, were added into ORCHIDEE-CAN (revision 4262) to simulate the critical wind speeds of all forest stands.The model simulates the critical wind speed for two types of damage: tree overturning and stem breakage.The critical wind speed for overturning is calculated as where CWS is the critical wind speed (m s −1 ), and the subscript ov denotes the critical wind speed for overturning.κ is von Karman constant and D is the inter-tree spacing (m).
C reg is a regression coefficient that was derived from tree- pulling experiments (N m kg −1 ; Nicoll et al., 2006).SW is the green weight of the bole of the tree (kg).SW is calculated by multiplying the model-simulated above-ground biomass with green density for different tree species (see Table 2).f CW is the enhanced momentum caused by the overhanging displaced mass of the canopy.In ORCHIDEE-CAN f CW was set to 1.136 (unitless), as suggested by Nicoll et al. (2006) from analysing extensive tree-pulling data.In other words, the applied turning moment is a constant fraction of the total turning moment.f edge is the edge factor (unitless) and a detailed description of the factor is given in Sect.2.2.2.G is the gust factor (unitless) and its calculation is also described in section Sect.2.2.2; h is the average tree height (m), d is the zero-plane displacement (m), z 0 is the roughness length (m), and all are simulated by ORCHIDEE-CAN (Naudts et al., 2015).Note that d and z 0 depend on the wind speed at canopy height u h , and hence iterations are required to solve Eq. (1) (see below).
The critical wind speed for stem breakage is calculated as where the subscript bk denotes stem breakage, MOR is the modulus of rupture (Pa) of green wood and was parameterized for different tree species (see Table 2), and dbh (m) is the tree diameter at breast height as simulated by ORCHIDEE-CAN; f knot is a factor to reduce wood strength due to the presence of knots (unitless).Similar to Eq. ( 1), d and z 0 in Eq. ( 2) also depended on the wind speed measured at canopy height u h , and hence iterations are required to solve Eq. ( 2) (see below).
The relationships between the aerodynamic parameters (d and z 0 ) and the vegetation structure follow the analytical relationships proposed by Raupach (1994): where c d1 is a regression parameter (c d1 = 7.5), C w is crown width (m), and C d is crown depth (m).Individual trees in the model are simulated as spherical elements, and the canopy width and canopy depth are thus identical; γ is a parameter that depends on the canopy characteristics (γ = and h is the atmospheric stability correction function ( h = ln(2) − 1 + 1 2 ).Tree crowns, branches, and stems are considered as porous and flexible materials that will streamline and thus change their shape with changing wind speed (u h ).Streamlining was parameterized through the parameters C and n, which were reported for wind tunnel experiments with different tree species (Rudnicki et al., 2004;Vollsinger et al., 2005) (see Table 2).The maximum value of h is set at u h = 10 m s −1 and the minimum value of C D is set at u h = 25 m s −1 (limits are based on wind speed range in Mayhead, 1973).The species-specific streamlining effect for wind speeds outside of this range was calculated by holding u h constant to its lower or higher threshold.In other words, C D was implemented as a kind of step function.
Critical wind speeds are calculated as the solutions to a non-linear set of equations for overturning, i.e.Eqs. ( 1), (3), and (4), and another set of equations for stem breakage, i.e.Eqs. ( 2), (3), and ( 4).An initial wind speed, u h =25 m s −1 , is applied to Eqs. ( 3) and ( 4) to obtain an approximation for CWS ov by applying Eq. ( 1).Similarly, an initial wind speed  (unitless) 6.0 6.0 6.0 6.0 6.0 of u h =25 m s −1 , is applied to Eqs. ( 3) and ( 4) to obtain an approximation for CWS bk by Eq. ( 2).Subsequently, u h is set to the value of CWS ov (or CWS bk ) to estimate the aerodynamic parameters (d and z 0 ) for the next iteration.The iteration process is stopped if the difference in CWS between two iterations falls below 0.01 m s −1 or the number of iterations exceeds a threshold of 20.
Whereas ORCHIDEE-CAN is designed to simulate both even-and uneven-aged stands, ForestGALES is currently limited to simulating the critical wind speeds for even-aged forests.Although this difference in design is thought to have few consequences, it is considered essential in the calculation of the ratio between tree height and tree spacing (the so-called inter-tree spacing, D).In an even-aged stand both tree spacing and tree height are homogeneous and thus well defined at the stand level.This is no longer the case for tree height in uneven-aged stands.This issue is accounted for by calculating a critical wind speed for each diameter class separately.Although this approach addresses the possible heterogeneity in tree height, it requires a value for inter-tree spacing, by definition a stand characteristic, to be calculated for each diameter class.To calculate the inter-tree spacing for each diameter class, firstly the total woody biomass at the stand level is calculated.Subsequently, this total woody biomass is divided by the biomass of the modelled trees in each diameter class.The outcome is considered the virtual inter-tree spacing of the diameter class D and was used in Eqs. ( 1) and (2) to calculate the critical wind speeds for each diameter class.
By default three diameter classes are used to describe the heterogeneity within a forest stand.ORCHIDEE-CAN then calculates the critical wind speed for breakage and overturning based on the vegetation structure parameters for each diameter class.When using three diameter classes, as is the case in this study, a total of six critical wind speeds are thus calculated for each forest in each grid box.Subsequently, the lowest critical wind speed is used to determine the dam-age type for each diameter class.The number of damaged trees in each diameter class is then calculated by multiplying the damage rate (D β ; see Sect.2.2.4) with the tree numbers within each diameter class.The total number of trees damaged by a storm was calculated as the sum of the damaged trees in each diameter class.

Gustiness and edge effect (ForestGALES)
ORCHIDEE-CAN is driven by half-hourly wind fields.For storm damage, such a time step is already too large because the half-hourly wind field hides the extreme wind gusts that occur within a half-hourly time step.Storm damage is determined more by the extreme gusts than by the average wind speed.This scaling issue is dealt with by explicitly simulating the gustiness through the so-called gust factor.The gust factor G was parameterized as a function of inter-tree spacing to tree height ratio and edge distance to the tree height ratio.These dependencies and their parameter values are based on wind tunnel experiments (Gardiner et al., 2000;Hale et al., 2012Hale et al., , 2015)).
where D h is constrained between 0.075 and 0.45, and x is parameterized as a function of leaf area index (LAI), i.e. x = 28 h LAI for the inner area near recent gaps and x = 9 h for the outer area forest away from such gaps (see Sect. 2.2.3).The length scale parameter x was derived from large eddy simulation (Pan et al., 2014).G adj is a linear scaling factor introduced to make the predicted gusts (based on wind tunnel experiments) better match the observed field measurements.
The edge effect was considered at the landscape level.Forests within a modelled grid were separated into two regions, inner area and the outer area, which can be simulated by the model (see Fig. S1 in the Supplement).At the inner area, the effect of vegetation structure on wind speed is accounted for through the edge factor f edge in Eqs. ( 1) and ( 2).The calculation of f edge follows the approach proposed by Gardiner et al. (2000).
At the outer area, the edge effect is negligible such that f edge is set to 1.0.

Vegetation structure
Vegetation structure is simulated at the landscape and the stand level (see Fig. 1).At the landscape level the simulations distinguish between forests with a newly formed forest edge and forest with established edges.Edges result from natural or anthropogenic stand-replacing disturbances.First, the surface area of stand-replacing disturbances is cumulated over the last 5 years (A 5 ), a time horizon corresponding to the time required for forests near newly formed edges to adapt to the increased gustiness (Gardiner and Stacey, 1996).By prescribing the average gap size (A gap ) to 2 ha and assuming gaps are square shaped and the gustiness is affected over a distance of 9 times the canopy height (9 h) (Gardiner and Stacey, 1996), the forest area that experiences an increased gustiness due to proximity to recent gaps (A inner ) is calculated as where the factor of 1 4 accounts for the fact that only the downwind edge perpendicular to the wind will experience an increased gustiness.The second term is the inner area generated by a single gap and the third term is the total number of gaps in the grid cell.The forest area that has no edges in its proximity A outer is calculated as the residual: where A grid is the area of the grid box being modelled.

Storm damage
With wind speeds approaching the critical wind speed, damage such as defoliation and branch damage become more likely.Once the wind speed exceeds the critical wind speed, overturning and stem breakage are possible but their likelihood increases with further increasing wind speeds.A sigmoid damage function is used to simulate the rate of storm damage to a forest.A similar approach has been applied and tested by Anyomi et al. (2017) for estimating storm damage as a function of the daily maximum wind speed.This relationship is formalized as where D β is the damage rate (unitless) and thus the share of trees that will be killed, D max is an observed maximum damage rate, which was set to 1.0, R f is a relaxation parameter to adjust the damage rate given by a certain wind speed below the model-calculated critical wind speed, and a value 6.0 was applied for all tree species.U max is the maximum daily wind speed from the forcing or model calculation.Subsequently, the lowest out of the six calculated critical wind speeds (see Sect. 2.2.1) is used to determine the damage type for each diameter class.The number of damaged trees in each diameter class is then calculated by multiplying the damage rate (D β ) with the tree numbers within each diameter class.The total number of trees damaged by a storm was calculated as the sum of the damaged trees in each diameter class.

Salvage logging
Damaged wood due to storms is often left on-site in unmanaged forests; however, salvage logging is often carried out for a managed forest in order to recover some of the economic losses and avoid large-scale insect outbreaks triggered by wind disturbance.When dealing with the effects of wind damage on the biomass pools of forests, the anthropogenic response to storm damage needs to be accounted for.ORCHIDEE-CAN distinguishes managed and unmanaged forest.In unmanaged forests all carbon contained in trees killed by wind storms ends up in the litter pools.For managed forests, a prescribed harvest ratio determines the wood that is salvage logged and the wood left on-site.In species that are prone to bark beetle attacks following wind throw, the volume left on-site is very small.In Sweden, a maximum volume of 5 m 3 ha −1 of newly damaged logs is allowed.However, following large-scale storm damage this threshold has been temporarily lowered to 3 m 3 ha −1 in order to reduce the risk of spruce bark beetle outbreaks.Following Gudrun, it was observed that no more than 1.8 to 1.1 m 3 ha −1 was left on-site.Given that the current implementation of storm damage was designed to deal with large wind storms, with a fair risk for subsequent bark beetle outbreaks (Kärvemo, 2015), applying a very high efficiency for salvage logging, i.e. 99 %, appears justified (Schroeder et al., 2006).

Soil characteristics
Although ORCHIDEE-CAN distinguishes three soil texture classes, i.e. sand, clay, and loam, the current approach to simulating soil water hydrology following de Rosnay ( 2003) Furthermore, ForestGALES distinguishes shallow-, medium-, and deep-rooting species.This classification was applied in ORCHIDEE-CAN through the parameter (humcte) describing the vertical root profile.In ORCHIDEE-CAN, the rooting density is assumed to following a function that decreases exponentially from the top to the bottom of the soil layers and is considered independent from site conditions or stand age.If 90 % of total root mass was found above a depth of 2 m, the species was considered shallow rooted.The effect of rooting depth on critical wind speeds was accounted for by using a different regressing coefficient C reg for shallow-and deep-rooting species (see Table 2).Under the current parameter settings of the rooting profile, ORCHIDEE-CAN considers all tree species to be shallow rooted.Note that in ORCHIDEE-CAN the rooting profile is also a critical parameter with which to calculate the drought stress of trees.
When the soil is frozen, ORCHIDEE-CAN only allows wind storm damage from stem breakage.In other words, there is no overturning of trees when the soil is frozen.The soil temperature at 0.8 m below the surface was used as the threshold to decide whether the soil was frozen or not.

Downscaling wind fields
In this study, simulations are forced by 6 h CRU-NCEP climate reanalysis (Viovy, 2011;Maignan et al., 2011).The internal weather generator of the ORCHIDEE-CAN model interpolates this reanalysis to obtain the half-hourly (30 min) mean wind speed used in the calculation of storm damage.Interpolation of the 6 h reanalysis is expected to dampen the wind speed and therefore wind damage calculated by ORCHIDEE-CAN would be underestimated.To overcome this issue a tuning parameter called mean wind ratio (MWR) was introduced.
The MWR converts the mean wind speed from 6 h CRU-NCEP reanalysis into a mean wind speed at the 30 min time step.ORCHIDEE-CAN then uses these daily maximum estimated values to calculate damage rates that may occur in a forest stand near or away from forest edges.Note that in this study the values for the mean wind ratio are specific for the CRU-NCEP half-degree forcing at a 6 h time step.The mean wind ratio will thus need to be re-parameterized if the wind driver is replaced by any other forcing.
The MWR was estimated from 38 European eddycovariance sites covered by forests for which the meteorological data were freely available.For the period 1996 to 2007, 208 site-year combinations were retained for which over 60 % of the half-hourly measurements were available.The remaining data gaps were filled with the ERA-Interim reanalysis data (Dee et al., 2011;Vuichard and Papale, 2015).The 208 site-years were extracted from the CRU-NCEP reanalysis.Subsequently, the ratio between the maximum 30 min mean wind speed within a 6 h period and the 6 h mean wind speed obtained for the same location and time frame from the CRU-NCEP reanalysis was calculated.
where the subscripts Fluxnet and CRU-NCEP denote the data source, U Fluxnet is the maximum value of the 30 min mean wind speeds (m s −1 ) within a 6 h time frame, and U CRU-NCEP is the 6 h mean wind speed from the CRU-NCEP dataset.Finally, maximum wind speeds in the time series of each forest site were stratified according to the wind force catalogue in Beaufort wind scale (BWS), an empirical measure that relates wind speed to observed conditions on land (National Meteorological Library and Archive, 2010).Data were then analysed to account for the relationship between the mean wind speed and the ratio between the observed and temporal average mean wind speed.
The quality of the fitted relationship was calculated as the root mean square error (RMSE).RMSE is a statistic that is widely applied to quantify the difference between values predicted by a model and the values actually observed: where Y denotes the observed values, F denotes the modelled values, and the subscription i is the sample index.RMSE is used throughout this study to quantify the goodness of fit between observed and predicted values.

Critical wind speeds for five tree species
The calculation of critical wind speeds is parameterized for five common tree species in Europe: Norway spruce (Picea sp.), Scots pine (Pinus sylvestris), beech (Fagus sylvatica), birch (Betula sp.) and oak (Quercus ilex and Quercus suber).
In ORCHIDEE-CAN these five tree species are simulated as separate PFTs over Europe.Parameters for the other PFTs within and outside of Europe were based on ForestGALES, which has been parameterized for 21 tree species including the five species simulated by ORCHIDEE-CAN.Table 2 lists the parameters used in ORCHIDEE-CAN to calculate the critical wind speeds.Spruce, pine, and birch make up almost the entire forest cover of southern Sweden.Pine is the single most important species in Les Landes.These regions were only used to test the model.The anticipated simulation domain of this new development is Europe.The five species for which the model was tested make up 67 % of the European forest cover.In terms of taxonomic families the representativeness increases to > 90 % (Koeble and Seufert, 2001).
3 Observational data, model tuning, and evaluation

Storm damage observations
A 60-year-long record of storm damage statistics over Sweden was extracted from the country-level European Forest Institute storm damage database from 1951 to 2010 (Nilsson et al., 2004;Schlyter et al., 2006;Bengtsson and Nilsson, 2007;Gardiner et al., 2010).We refined this dataset with regional information.Following Gudrun in January 2005, the Swedish Forest Agency mapped the spatial distribution of the storm damage for damage classes ranging from 5 to 50 m 3 per ha in steps of 5 m 3 for the first damage class and 10 m 3 for all subsequent classes (see Fig. 7).The storm made landfall near the south of Norway, went through the north of Götaland, and resulted in extensive forest damage in the central area of Götaland.The area around the cites of Ljungby and Växjö was reported as having the greatest damage of about 30 m 3 ha −1 .The spatial extent of storm damage was retrieved from ocular inspection of aircraft images processed by the Swedish Forest Agency.In January 2009, Klaus made landfall in south-western France near Les Landes forest and damaged 43.1 × 10 6 m 3 of wood in France.The dynamics of surface albedo and LAI between 2001 and 2010 were extracted from MODIS and SPOT-VGT satellite images, receptively (Planque et al., 2017).These remote sensing time series were overlaid by ORCHIDEE-CAN simulations and compared at two locations between 2001 and 2010, thus resulting in a total of 20 data points (shown as pink arrows in Fig. S2).

Model tuning for storm damage
Although all parts of the wind-throw module come with their own assumptions, parameters, and subsequent uncertainties, the calculation of the gustiness (G; Eq. 5) is considered to be among the most uncertain part of the model because it involves spatial and temporal scaling issues in both the driver data and the model formulation.Furthermore, the function that was used to convert the difference between the critical wind speed and wind speed into a damage rate (D β ; Eq. 9) is also thought of as very uncertain.The key parameters in these functions, G adj and R f , are empirical, lack a good observational constraint, and are therefore prime parameters to be tuned for matching the observations.
Given the availability of 60 years of observations of primary damage caused by wind storms between 1951 and 2010 for the region (Gardiner et al., 2010), southern Sweden was selected as the study area for model tuning.Tuning made use of the observed damage volumes for the years 1981 to 2000 because it is the most recent period without major storms.The storm Gudrun (2005) was deliberately excluded from the tuning period so it could be used as an evaluation of model performance.
The simulations used for parameter tuning started in 1981 and therefore had to be forced by a spatially explicit description of the biomass distribution in southern Sweden at the end of the year 1980.In this study, the spatially explicit biomass distribution for 1980 was extracted from a previous study that simulated the effects of forest management and land cover changes in Europe between 1750 and 2010 (McGrath et al., 2015;Naudts et al., 2016).The previous study did not, however, account for the effects of wind storms on woody biomass and is therefore likely to overestimate the biomass in regions with chronic wind stress.This issue was overcome by starting the simulation in 1971 rather than 1981 and by running ORCHIDEE-CAN with the storm damage module between 1971 and 1980 to adjust the biomass to chronic wind stress.Chronic wind stress occurred mainly in western Norway.Subsequently, this adjusted spatially explicit description of the biomass distribution was used in the simulations for parameter tuning.Trail-and-error tuning started on 1 January 1981 and continued until 31 December 2000 by forcing the model with the CRU-NCEP reanalysis to simulate the primary damage caused by storm events including the 1999 Anatol storm.

Critical wind speeds for five tree species
Model implementation and parameterization were tested for all five tree species shared between ForestGALES and ORCHIDEE-CAN (see Sect. 2.3.2) by running ORCHIDEE-CAN for a test pixel.The Fontainebleau forest, which is the closest large forest to the LSCE research institute, was arbitrary selected as the simulation site and ORCHIDEE-CAN was run for 200 years by cycling over the CRU-NCEP reanalysis data from 1901 to 1930.Subsequently, the canopy structure variables simulated by ORCHIDEE-CAN, including h, D, C w , C d , and LAI, were used as the input data for a stand-alone version of ForestGALES (MathCAD version) to estimate the CWS ov and CWS bk .In total four types of CWSs, CWS ov and CWS bk for the inner and outer forest, were simulated by both models for the same canopy structure.

Critical wind speeds over southern Scandinavia
Southern Scandinavia was simulated as a 35 by 35 halfdegree pixel grid.For each pixel, four critical wind speeds, i.e.CWS ov and CWS bk for the inner and outer area of the forest, were simulated making use of the simulated canopy structures for the different tree species and age classes contained in that pixel.To this aim a simulation over southern Scandinavia was set up by using the CRU-NCEP reanalysis for 2010 making use of the parameter values for G adj and R f , www.geosci-model-dev.net/11/771/2018/Geosci.Model Dev., 11, 771-791, 2018 q q q q q q q q q q q 5 10 15 20 0 2 4 6 8 q q q q q q q q q RMSE: 0.   12) with regression coefficients a 0 = −5.299, a 1 = 2.051, a 2 = −0.191,and a 3 = 0.006.This relationship is used to convert CRU-NCEP 6 h mean wind speed to the 30 min maximum wind speed in this study.The RMSE using this regression model to predict the mean value of MWR in each BWS class is 0.48.The open circles (grey colour) show the effect of 6 h temporal aggregation on the MWR from the selected Fluxnet sites in European forest.The grey line is the fitting line of the open circles.
as derived during the tuning phase.The model was restarted from the time point 31 December 2010 of the aforementioned simulation of Naudts et al. (2016), in which forest management and land use were prescribed from the historical reconstructions presented in McGrath et al. (2016).The vegetation over this area mainly consists of coniferous tree species, i.e.Norway spruce (Picea abies (L.) Karst) and Scots pine (Pinus sylvestris L.); however, a species mixture of coniferous tree and broadleaved species such as birch (Betula pendula Roth and B. pubescens Ehrh.) can be found in this region (Drössler, 2010).The simulated spatial distribution of critical wind speeds for each tree species thus reflects the effects on critical wind speeds of age class structure, forest management, and, to a lesser extent, local climate conditions.

Storm damage over Sweden
The capability of ORCHIDEE-CAN to simulate storm damage was evaluated by comparing the observed and simulated storm damage over Sweden from 1951 to 1980 and from 2001 to 2010.A simulation over Sweden was set up by using the CRU-NCEP reanalysis over the last half-century making use of the parameter values for G adj and R f as derived during the tuning phase.The region under study entails a 15 by 15 half-degree pixel grid covering southern Sweden and a section of Norway.The state of the forest in Sweden on 31 December 1940 was described by using the matching year from an existing 250-year-long simulation (see also Sect.3.2).Because the 250-year-long simulation did not account for the effects of wind damage, the first 10 years from 1941 to 1950 of the simulation experiment were intended to reach equilibrium between the vegetation and the mean wind speed.For these 10 years the climate data for 1941 to 1950 were used.Within the study domain, this approaches reduced the standing biomass mainly in western Norway (not shown).Subsequently, the simulation experiment continued from 1951 to 2010, the period for which damage reports are available (Nilsson et al., 2004;Schlyter et al., 2006;Bengtsson and Nilsson, 2007;Gardiner et al., 2010).The years 1981 to 2000, which were used for parameter tuning, were excluded from the evaluation.

Biophysical effects of storm damage in Les Landes
The capability of ORCHIDEE-CAN to simulate the effects of storm damage on the surface albedo in the visible range was evaluated.In January 2009, the storm Klaus passed over Les Landes in France.Simulated albedo and leaf area estimates before and after this storm were compared against observed albedo and leaf area values before and after the storm.Les Landes was selected as a study site because time series for albedo were available for the region (Planque et al., 2017).The region under study is covered by a 6 by 6 halfdegree pixel grid and the simulation is driven by the CRU-NCEP reanalysis between 1991 and 2010.The state of the forest on 31 December 1990 was described by using the matching year from the same existing 250-year-long simulation (see previous section).The first 10 years of the simulation experiment were intended to reach equilibrium between the vegetation and the mean wind speed.For these 10 years the climate data for 1991 to 2000 were used.Within the study domain, this approach had very little effect on the standing biomass (not shown).Subsequently, the simulation experiment continued from 2001 to 2010, the period for which the MODIS albedo product was analysed (Planque et al., 2017).Two ORCHIDEE-CAN simulation pixels, one located at the centre of Les Landes forest with storm disturbance and another located in the east of Les Landes forest without the storm disturbance, were selected for analysing the biophysical responses to the storm damage for the summertime during the 10-year simulation.
Geosci.Model Dev., 11, 771-791, 2018 www.geosci-model-dev.net/11/771/2018/q q q q q q q q q q q q q q q 8 10 12 14 q ORCHIDEE ov ORCHIDEE bk ForestGALES ov ForestGALES bk RMSE,ov: 0.893 RMSE,bk: 1.29 q q q q q q q q q q q q q q q 8 10 12 14 RMSE,ov: 1.44 RMSE,bk: 1.09 q q q q q q q q q q q q q q q 8 10 12 14 RMSE,ov: 1.000 RMSE,bk: 0.634 q q q q q q q q q q q q q q q 8 10 12 14 RMSE,ov: 0.345 RMSE,bk: 0.664 q q q q qqq q q q q q q q q 8 10 12 14 RMSE,ov: 0.722 RMSE,bk: 0.757 q q q q qq q q q q q q q q q 8 10 12 14 RMSE,ov: 0.579 RMSE,bk: 0.789 q q q q q q q q q q q q q q q 8 10 12 14 RMSE,ov: 0.413 RMSE,bk: 0.796 q q q q q q q q q q q q q q q 8 10 12 14 qq qq qq qq q q q q q q q q q q q q q q q q q q q q q 8 10 12 14 qq qq qq qq q q q q q q q q q q q q q q q q q q q q q 8 10 12 14 .Model-calculated critical wind speeds as a function of tree height for five common tree species in Europe.For each tree species the critical wind speed is calculated for overturning and stem breakage for forest located away from a forest edge (outer, a-e) and near a forest edge (inner, f-j).The ORCHIDEE-CAN simulations are shown as symbols and benchmarked against the ForestGALES simulations, which are shown as lines.The scientific names of the tree species are given in Sect.2.2.3.The CWS differences between the ORCHIDEE-CAN and ForestGALES calculation were calculated using Eq. ( 11), for which the CWSs from ForestGALES were treated as the reference.

Downscaling wind fields
The effect of spatial and temporal aggregation on wind fields can be found by the comparison of CRU-NCEP 6 h mean wind speed to the stand-level 30 min maximum wind speed.Due to spatial and temporal aggregation, the CRU-NCEP 6 h mean wind speed reanalysis consistently underestimates the 30 min maximal wind speeds.Stratification by the Beaufort wind scale revealed a clear relationship between this scale and the mean wind ratio for different Beaufort values (see Fig. 2).The low value of 6 h CRU-NCEP wind speed in the BWS10 bin may be due to the fact that all observations in this bin come from a single location wind category (Fig. S3).
A fourth-order polynomial function was fitted through the observed MWR to downscale the CRU-NCEP 6 h mean wind speed to a max value of 30 min mean wind speed within a 6 h time frame.
where U max is the max value of 30 min mean wind speed (m s −1 ) within a 6 h period; a 0 to a 3 are regression param-eters, which have the following values: a 0 = −5.299, a 1 = 2.051, a 2 = −0.191,and a 3 = 0.006.Given the intended use of this analysis to convert the 6 h mean wind speed of the CRU-NCEP reanalysis into a likely wind speed used by ORCHIDEE-CAN on the half-hourly timescale, the mean ratios of the maximum wind speeds were extracted for each Beaufort wind scale for a 6 h averaging period.The value of MWR for a 6 h averaging period increased from 1.0 to 6.0 when the Beaufort wind scale went up from scale 1 to scale 11.

Critical wind speeds for five tree species
Given that ForestGALES has been extensively validated (Hale et al., 2015) and that the wind-throw module in ORCHIDEE-CAN closely followed the physical processes and empirical formulations from ForestGALES, ORCHIDEE-CAN is expected to simulate similar critical wind speeds as estimated by ForestGALES.This is indeed the case for the five tree species sharing the species parameters between ForestGALES and ORCHIDEE-CAN.matched.The difference in CWSs from the two models (RMSE for all CWSs) was limited to a few metres per second ranging from 0.4 to 2.0 for different tree species.Moreover, the estimates show that the critical wind speed close to a forest edge is always lower than the respective value further away from a forest edge.(e.g.Fig. 3f < a, ..., Fig. 3j < Fig. e).Also note that for oak, for example, the critical wind speeds for breakage and overturning are almost identical for small trees, but the difference between both critical wind speeds increases with increasing tree height (Fig. 3e and j).This implies that taller oak stands are more vulnerable to tree overturning than to stem breakage compared to smaller stands.The relationship between tree height and critical wind speed for spruce is very different compared to oak.For spruce the critical wind speeds for breakage and overturning are within several m s −1 from each other and this difference remains more or less constant with increasing tree height (Fig. 3a and f).In other words, both stem breakage and tree overturning may occur simultaneously in tall spruce stands.Furthermore, small beech stands are more vulnerable to breakage (similar to spruce) but tall beech stands are more vulnerable to turnover (similar to oak; Fig. 3c and h).

Critical wind speed over southern Scandinavia
The modelled critical wind speeds for inner and outer forest areas are shown in Fig. 4 and these CWSs were compared with daily maximum wind speed.The contour lines in Fig. 4c and d indicate the spatial distribution of the difference between CWSs and daily maximum wind speed.Across southern Scandinavia, most of the Norway spruce forests can cope with wind speeds exceeding 25 m s −1 for the outer area (away from forest edge; Fig. 4a), but this minimum critical wind speed can be reduced around 5 m s −1 for the inner area (in a forest edge; Fig. 4b).Stands with old spruces of 25 m tall growing on a well-drained soil were simulated to even Geosci.Model Dev., 11, 771-791, 2018 www.geosci-model-dev.net/11/771/2018/q q q q q q q q q q q q q q q q q q q q 1981 1982 1983 q q q q q q q q q q q q q q q q q q q q 0 10 20 30 40 Damaged wood volume ( x 10 m ) q q q q q q q q q q q q q q q q q q q q 1981 1982 1983 sustain wind speeds exceeding 40 m s −1 .The proximity of forest edges is expected to decrease the critical wind speeds (Eqs. 1, 2, and 6).Indeed, a spruce forest close to a forest edge can sustain wind speeds ranging from 10 to 35 m s −1 .As already shown in Fig. 3a and f, the difference in critical wind speeds between overturning and breakage are very small for spruce.This behaviour is sustained over a large spatial domain (compare Fig. 4a and b).We further compare the difference between these CWS maximum daily wind speeds on 9 January 2005.Forests located in the centre of southern Sweden are expected to suffer from storm damage for the Gudrun case (see Fig. 4c, d)

Storm damage over Sweden
The frequency of storm events in the period from 1981 to 2000 used for parameter tuning is about one storm every 5 to 10 years.The observed primary damage of these storms ranges between 2 and 5 × 10 6 m 3 of wood (black line in Fig. 5).The sensitivity of the simulated damage for different values of the relaxation parameter R f is presented in Fig. 5 (grey lines).Using the prior setting, i.e.G adj = 1 and R f = 1, results in an error (RMSE; Hyndman and Koehler, 2006)  q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0 20 40 60 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q −2 − of about 2 × 10 6 m 3 for a calibration period without major storms.Higher parameter values, e.g.G adj = 3.0, make the model overestimate storm damage by up to 10 × 10 6 m 3 .Parameter tuning (G adj = 2.45 and R f = 6.0) resulted in a modest reduction of the deviation between observations and simulations, i.e. an RMSE of 1.35 × 10 6 m 3 , and was therefore used to evaluate the wind-throw model.
The RMSE was largely determined by the simulated damage for 1986, 1987, and 2000.For these years, no damage above the background volume was observed despite the presence of high wind speeds in the CRU-NCEP reanalysis over southern Sweden in December and January for those three years.Regardless of the value used for R f , the high wind speeds in the CRU-NCEP reanalysis resulted in substantial damage volumes compared to the observations.This result suggests that ORCHIDEE-CAN would benefit more from improving the representation of gustiness G adj than from refining the damage relaxation parameter R f .The forest extent has expanded strongly in southern Sweden during the 20th century.In 1940 the forestry system was different from that of today.The clear-felling system that produces the forest edges modelled in this study was not widely adopted until the 1950s.This means that the state of the forest and the exposure to strong wind was very different in the beginning and at the end of the simulation period.ORCHIDEE-CAN partly simulates this transition: in 1940 the average simulated above-ground biomass was 154 m 3 ha −1 , while in 1990 the biomass had increased to 168 m 3 ha −1 .Within this half-century the simulated forest area in Sweden increased from 242 500 to 247 100 km 2 and the share of high stand management increased from 54.8 to 69.2 %.Despite ORCHIDEE-CAN accounting for the observed forest transition, the model mostly underestimated the damaged wood volumes of large storms.
In January 2005, Gudrun caused the biggest recorded storm damage in Swedish history.We extracted the simulation result from the evaluation experiment and calculated the total simulated storm damage by Gudrun.Simulated storm damage was 76.6 × 10 6 m 3 , which compares well to the reported 75.0 × 10 6 m 3 (Gardiner et al., 2010).Moreover, the high-level damage pixels (> 30 m 3 per ha) are located in central Götaland, which demonstrates that when the driver data locate the storm in the right location, the model has the capability to reproduce the spatial distribution for the event (see Fig. 8).
The model, however, underestimated the storm damage in 2007 (Fig. 5).Given the reconstructed wind speed in the driver data, the observed storm damage appears high.The failure of the model to simulate the right order of magnitude in storm damage in 2007 may be for two reasons: (1) the CRU-NCEP reconstruction underestimated the wind speed and/or (2) ORCHIDEE-CAN overestimated the critical wind speed by not accounting for the legacy effects of Gudrun, such as decreased tree stability owing to root damage or increased gap sizes following Gudrun.Additionally, the observed storm damage is more uncertain than for Gudrun because the inventory in 2007 was not made as detailed as for Gudrun.

Biophysical effects of storms
In January 2009, Klaus made landfall in south-western France near Les Landes forest.Following Klaus a decrease in leaf area index and an increase in albedo were reported from SPOT-VGT and MODIS observations, respectively (Planque et al., 2017).Comparison of the simulations against the MODIS observations suggests that the model is able to reproduce the direction of the changes in the surface characteristics following the passing of Klaus in 2009 (Fig. 9a-d).ORCHIDEE-CAN reproduced both the magnitude and direction of the changes in leaf area index (Fig. 9a).Whereas ORCHIDEE-CAN rightly simulated an increase in visible albedo following storm damage (Fig. 9b), the model overestimated the absolute value of albedo and its change.Although this mismatch could be caused by numerous processes and parameters, the use of a constant PFT-specific background albedo in ORCHIDEE-CAN may limit the capability of the model to correctly simulate the growth of an herb layer following stand-replacing disturbances.In line with meteorological theory, large-scale forest disturbances resulted in a decrease in roughness length and transpiration Fig. 9c-d).Decreasing roughness and transpiration could not, however, be confirmed at the site level due to a lack of observations.Geosci.Model Dev., 11, 771-791, 2018 www.geosci-model-dev.net/11/771/2018/ 5 Discussion

Downscaling wind fields
Wind is a highly heterogeneous turbulent flow of air.The flow consists of various sized eddies, and because the energy is conserved, changes in the distribution of the eddy size will result in different flow regimes.Gusts are generated when the flow passes over a heterogeneous surface such as valleys, ridges, objects, or any surface warmer or cooler than its surroundings (Lanquaye-Opoku and Mitchell, 2005;Liu and Weng, 2009).The temporal resolution of a gust is seconds to minutes and its spatial resolution is metres.The highest wind speeds are caused by these gusts and are responsible for most of the storm damage (Usbeck et al., 2010).
The success of a large-scale model in simulating storm damage will thus to a large extent depend on the capability of the model to simulate gusts.At present, large eddy simulations are the most advanced approach to simulate the turbulent flow of air (Moeng, 1984;Dupont and Brunet, 2008).The high computational demand of the method (Yang, 2015) currently prevents it from being used in large-scale models such as ORCHIDEE-CAN.Likewise, advances in regional atmospheric modelling resulted in the capacity to simulate gusts and gust gradients in the wind field such as tornados (Ishihara et al., 2011) and hurricanes (Vickery et al., 2009).The computational requirements of regional models is also too high to consider them for global simulations.High computational needs can be avoided by using statistical downscaling but these methods come at the expense of a poor process representation (Larsén and Mann, 2006;Salameh et al., 2009;Huang et al., 2015;Winstral et al., 2017).Moreover, the statistical relationship used for downscaling the gridded wind field from a coarse to a fine resolution depends on the gridded wind field and its spatial and temporal resolutions.
In this study, we used a statistical approach that builds on the relationship between in situ observations and the gridded reanalysis data.The simplicity of the approach enabled us to focus on implementing, developing, and validating the storm damage model rather than improving the physical representation of gusts in ORCHIDEE-CAN.Downscaling the gridded wind field from a coarse resolution to a fine resolution requires accounting for both spatial and temporal aggregation.Our approach made use of the mean wind ratio to convert the mean wind speed from 6 h CRU-NCEP reanalysis into a mean wind speed at the 30 min time step.By further analysing the contribution of spatial and temporal aggregation to downscaling (see Fig. S3), it was found that the temporal averaging was responsible for about 40 % of the reduction of the variation, whereas spatial averaging was responsible for the remaining 60 %.In other words, simulating storm damage in large-scale models would initially benefit most from increasing the spatial resolution as that would allow us to better account for the extreme wind speeds.
Note that in ORCHIDEE-CAN (revision 4262) the downscaled wind speed is only used to test whether the critical wind speed is exceeded and, if so, to simulate the resulting damage from stem break and overturning.Consequently, the downscaled wind speeds are not used for calculating surface roughness, momentum, heat exchange, and water vapour transfer.The current statistical downscaling is thus not suitable for use in model set-ups in which ORCHIDEE-CAN is coupled to an atmospheric circulation model, for example, LMDZ (Hourdin et al., 2006).Simulating storm damage from coupled land-atmosphere models, e.g.LMDZ-ORCHIDEE-CAN, will thus require an effort to better represent turbulent air flow in the planetary boundary layer.

Storm damage over Sweden
The state of the forest in Sweden on 31 December 1940 was described by using the matching year from an existing simulation, which ran from 1750 to 2010 (Naudts et al., 2016).This existing simulation accounted for land cover changes and changes in forest management following the historical reconstruction of forest cover and management by Mc-Grath et al. (2015), but it did not account for storm damage.However, the forest reconstruction map shows that the forest species over Sweden mainly consist of Scots pine, Norway spruce, and birch, which is comparable with the report made by Helmfrid (1991).
Given its strength, the 1999 storm Anatol resulted in relatively little damage (5 × 10 6 m 3 ), likely because it hit the southernmost part of Sweden where the landscape contains fewer forests and the share of broadleaf forest is much higher than just a few tens of kilometres north of the storm track.Although the good match between the data and the model (Fig. 5) is due to the fact that Anatol occurred within the tuning period, the conditions that are held responsible for the low volume of storm damage were indeed found to be reproduced.The CRU-NCEP reanalysis suggests the Anatol storm track hitting southern Sweden, and ORCHIDEE-CAN uses a forest index of 56 % in the storm track, but 68 % for landscapes 100 km north of the track.Furthermore, in 1999, 67 % of the simulated forests were broadleaved in southern Sweden, whereas roughly one-third were 100 km north of the track.Finally, ORCHIDEE-CAN simulates higher critical wind speeds for broadleaved trees (especially in winter) than for conifers (Fig. 4).
Simulated storm damage between 1951 and 2010, excluding the years 1981 to 2000 which were used for parameter tuning, was used to evaluate the model.Within this period, major wind storms occurred in 1954, 1969, and 2005 and were associated with observed storm damage ranging from 20 to 75 × 10 6 m 3 .Note that the storm damage reported for the evaluation events is about 10 times higher than the damage reported for the events used to parameterize the model.The model error is 5.05 × 10 6 m 3 in the evaluation period.When the large storms in 1954, 1969, and 2005  in the calculation of the RMSE, the model errors decrease to 1.62 × 10 6 m 3 (see Fig. 6).The RMSE is thus largely due to underestimating the damage during big storms.

Subpixel heterogeneity
Gusts and thus storm damage are often generated by surface heterogeneities, for example, forest edges, surface topography, and landscape features, with a very different temperature than the surrounding landscape.Whereas these heterogeneities can be accounted for in large eddy simulations and regional models, large-scale models such as ORCHIDEE-CAN managed to limit their computational costs by, among other simplifications, ignoring these subpixel heterogeneities (Krinner et al., 2005).When implementing processes that are partly driven by these heterogeneities, i.e. storm damage, these models are thus operated at the limit of their design.
In this study, we have tried to overcome this issue by reconstructing subpixel heterogeneity from the wood harvest ag-gregated over the last 5 years and assuming that all gaps are square-shaped and have a surface area of 2 ha (see .This approach enabled ORCHIDEE-CAN to calculate separate critical wind speeds for forest close to a forest edge (inner) and forests away from a forest gap (outer).Although this approach is thought to have contributed to reproducing the observations, it lacks at least three welldocumented subpixel heterogeneities: (1) in ORCHIDEE-CAN the number of gaps varies over time but their surface area is set as constant.Hence, individual gaps do not increase in size and the relationship between gap size and gustiness (Peltola, 2006) is not accounted for.(2) Although each pixel has an altitude, elevation differences within a pixel are not considered.Thus, ORCHIDEE-CAN does not account for the effects of exposure on tree species distribution (Ruel et al., 1997).(3) In ORCHIDEE-CAN, all tree species within a pixel share the same water column, and hence the model does not account for the interactions between tree species and long-term soil water content heterogeneity (Ringeval et al., Geosci. Model Dev., 11, 771-791, 2018 www.geosci-model-dev.net/11/771/2018/2012).All three omissions may have contributed to the failure of the model to simulate the storm damage in 2007 (see Fig. 6).
It is expected that simulating gap size and gap dynamics would improve the performance of the storm damage model.Performance improvements, however, would equally apply to simulating fires because the size, age, and composition of subpixel forest patches were identified as important elements in simulating fire risk and fire behaviour (Keane et al., 2004).Further improvements can be expected for simulating forest regeneration because gap size will determine the growing environment for the recruitment (Rüger et al., 2009).It has even been suggested that to reduce the uncertainty in terrestrial vegetation responses to future climate and atmospheric CO 2 , a change in research priorities away from biomass production and toward structural dynamics and demographic processes should be favoured (Friend et al., 2014).Simulating subpixel heterogeneity is emerging as the next frontier in large-scale models; however, the challenge to do so at low computational costs remains.

Conclusions
The representation of storm damage in ESMs could use empirical models building upon the relationships between total storm damage and predictor variables such as topography, vegetation, soil, historical forest management maps, and wind exposure (Scott and Mitchell, 2005;Kamimura et al., 2008;Lagergren et al., 2012;Moore et al., 2013;Takano et al., 2016).As an alternative that better suits the process-based philosophy of Earth system modelling, more mechanistic approaches could be used.A mechanistic approach, however, requires information on the canopy structure (Peltola et al., 1999;Gardiner et al., 2000).This requirement rules out including mechanistic wind damage models in large-scale land surface models because the latter make use of the big-leaf assumption to represent the canopy.Recently, the big-leaf assumption was replaced by a threedimensional representation of the canopy structure in the large-scale land surface model ORCHIDEE-CAN (Naudts et al., 2015).This change in canopy representation enabled the incorporation of the process formalizations of the standlevel wind risk model called ForestGALES (Gardiner et al., 2000) into ORCHIDEE-CAN.
Incorporating ForestGALES into ORCHIDEE required solving three issues related to differences in spatial scales between the two models.The spatial resolution of Forest-GALES is 10 3 m 2 , whereas the typical resolution used in applications with large-scale land surface models such as ORCHIDEE-CAN is 10 9 m 2 .(1) This difference in spatial scales implies that subpixel heterogeneity had to be accounted for in ORCHIDEE.Although not all sources of subpixel heterogeneity were simulated, vegetation structure at the landscape level was accounted for (see Sect. 2.2.3).
(2) The spatial scale of ORCHIDEE required that wind fields had to be downscaled to account for the occurrence of gusts.A statistical downscaling approach was used in this study (see Sect. 2.3.1).Model performance could, however, benefit from a more mechanistic approach to downscale the wind fields to better account for gusts, regardless of whether the wind fields come from gridded reanalysis data or simulations from an atmospheric circulation model coupled to ORCHIDEE-CAN.A more mechanistic approach to downscale the wind fields may require accounting for subpixel heterogeneity, but it is also needed to enable the use of the storm damage functionality in tandem with atmospheric models, such as the WRF model (Stéfanon et al., 2012) or LMDZ (Hourdin et al., 2006).( 3) When ForestGALES is run on a small spatial scale it is fair to assume that all trees will be damaged when the critical wind speed for that stand is exceeded.This assumption no longer holds for the spatial scales used in ORCHIDEE-CAN.The model was therefore completed by an empirical function to convert the difference between the wind speed and the critical wind speeds into forest damage.
The storm damage functionality of ORCHIDEE-CAN was then parameterized and tested for selected case studies in southern Sweden and south-western France.The model largely captured the 60-year long-term storm damage dynamics over the Swedish forests and simulated a decrease in leaf area and an increase in visible albedo following storm damage in France.The model was thus shown to have the flexibility to reproduce diverse observations, although the validity of the parameters outside Sweden and France still needs to be demonstrated.
In the long term, building the capacity to simulate the impact of wind storms will result in a better understanding of the climate response to the biotic and abiotic disturbance agents in the Earth system, support actionable science to evaluate the effects of changing storm frequency and intensity on global food production, and improve the effectiveness of forest management in mitigating and adapting to climate change.

Figure 1 .
Figure 1.Information flow of this study showing the link between the different elements presented in Sects.2.2 and 2.3 (numbers in the figure refer to section numbers in the text).The diagram shows input data in blue.The dashed box shows critical wind speed calculations according to ForestGALES.

Figure 2 .
Figure2.Distribution of the mean wind ratio (MWR) in each Beaufort wind scale (BWS) and the relationship between the 6 h CRU-NCEP reanalysis wind speed and MWR.The fitting of the relationship (red line) used Eq.(12) with regression coefficients a 0 = −5.299, a 1 = 2.051, a 2 = −0.191,and a 3 = 0.006.This relationship is used to convert CRU-NCEP 6 h mean wind speed to the 30 min maximum wind speed in this study.The RMSE using this regression model to predict the mean value of MWR in each BWS class is 0.48.The open circles (grey colour) show the effect of 6 h temporal aggregation on the MWR from the selected Fluxnet sites in European forest.The grey line is the fitting line of the open circles.
Figure3.Model-calculated critical wind speeds as a function of tree height for five common tree species in Europe.For each tree species the critical wind speed is calculated for overturning and stem breakage for forest located away from a forest edge (outer, a-e) and near a forest edge (inner, f-j).The ORCHIDEE-CAN simulations are shown as symbols and benchmarked against the ForestGALES simulations, which are shown as lines.The scientific names of the tree species are given in Sect.2.2.3.The CWS differences between the ORCHIDEE-CAN and ForestGALES calculation were calculated using Eq.(11), for which the CWSs from ForestGALES were treated as the reference.

Figure 4 .
Figure4.The ORCHIDEE-CAN-calculated lowest critical wind speeds for overturning or stem breakage for forest located near (inner) and away (outer) from a forest edge.When making the display, the critical wind speeds from the three diameter classes and four age groups from Picea species were assessed and the lowest value was compared against the daily maximum wind speed for estimating the damage due to storms.Lowest critical wind speeds in the forest away from a forest edge (outer) (a), lowest critical wind speeds in the forest near a forest edge (inner) (b), lowest critical wind speeds overlaid with the difference between maximum daily wind speed and lowest critical wind speed in outer area on 9 January 2005 (c), and lowest critical wind speeds overlaid with the difference between maximum daily wind speed and lowest critical wind speed in the outer area (d).The contours show the positive wind speed difference in black and the negative wind speed difference in red.Forests within the red contours are expected to suffer from storm damage.

Figure 5 .
Figure 5. Sensitivity of the simulated storm damage over Sweden between 1981 and 2000 for different values for the relaxation parameter (R f ) and the gust factor adjustment G adj (a).Observed storm damage is extracted from Nilsson et al. (2004),Schlyter et al. (2006),Bengtsson and Nilsson (2007), andGardiner et al. (2010).The relative model simulation error (calculated as (estimation−observation) / estimation) for the best-tuned case (b).

Figure 6 .
Figure 6.Comparison of the storm damage simulated by the ORCHIDEE-CAN and the reported wood damage over Sweden from 1951 to 2010.Observed storm damage is extracted from Nilsson et al. (2004),Schlyter et al. (2006),Bengtsson and Nilsson (2007), andGardiner et al. (2010).The dashed-line area is the period from 1981 to 2000, which was selected for parameterization.The RMSE of the estimated storm damage is 1.35 × 10 6 m 3 for the parameterization period and 5.05 × 10 6 m 3 during the evaluation period.The validation period ranges from 1951 to 2010 but excludes the years 1981 to 2000.The relative model simulation error for the validation period from 1981 to 2010 (b).

Figure 7 .
Figure 7.The case study of the storm Edwin (Gudrun) on 8 January 2005.The colour scale shows the damaged woody volume (m 3 ha −1 ) due to this event.

Figure 8 .
Figure 8. Simulated spatial distribution of the damaged wood volume (m 3 ha −1 ) over southern Sweden by Gudrun in January 2005.

Table 1 .
Symbolic notation used throughout the paper.

Table 2 .
Parameter values used in the ORCHIDEE-CAN wind-throw module.The scientific names of the tree species are given in Sect.2.3.2.
(Hale et al., 2015)aining at the bottom of the soil layers.This assumption differs from ForestGALES in which four soil classes with different drainage are distinguished(Hale et al., 2015): freely draining mineral soil, gleyed (i.e.waterlogged and lacking in oxygen) mineral soil, peaty mineral soil, and deep peat.At present, ORCHIDEE-CAN only uses the ForestGALES parameters for freely draining mineral soils.Owing to this assumption ORCHIDEE-CAN is expected to overestimate the critical wind speed, resulting in less damage, for locations with shallow and/or wet soils.