Predicting snow damage in conifer forests using a mechanistic snow damage model and high-resolution snow accumulation data

ABSTRACT

Forest damage caused by heavy wet snow accumulation in the canopy is the second most important abiotic forest disturbance agent in Nordic conifer stands after wind. The extent and frequency of snow damage in the future climate in the Nordic region is a major uncertainty. Few mechanistic models of snow damage risk to trees exist that could support forest management scenario analysis and decision making. We propose a snow damage risk model consisting of a numerical weather prediction-based snow accumulation model for forest canopies and a mechanistic critical snow load model. Snow damage probability predictions were validated on snow breakage data from the winters of 2016 and 2018 covering 3.5 million individual trees in south-eastern Norway derived from pre- and post-damage aerial laser scanning campaigns. The proposed model demonstrated satisfactory damage and no-damage class separation with an AUC of 0.72 and 0.77 in Norway spruce and Scots pine, respectively, and an F1 score of 0.7 in conifers taller than 10 m that suffered moderate stem breakage. The model achieved a classification accuracy that is comparable to that of statistical models but is simpler and requires fewer inputs.

KEYWORDS:

• Snow damage
• mechanistic model
• critical snow load
• snow accumulation

Introduction

Heavy wet snow as a forest disturbance agent

Forest damage in the form of both stem breakage and uprooting caused by heavy wet snow accumulation in the tree crowns has occurred several times in Norway in recent years. For instance, in early 2018 along the south-eastern coast, then in 2020–2021 both in inland areas and in the coastal hills in the south-east, and again as recently as 2022–2023. The winter of 2020–2021 brought the greatest amount of snow damage over the past decade, with a conservative estimate of economic loss at EUR 3 million (NRK Vestfold og Telemark 2021). To add some context, the few existing estimates for the entire European Union range from 1 to 6 million m³ (3% of the total damage) annually over the past 100 years (Nykänen et al. 1997; Schelhaas et al. 2003); in Finland, annual insurance payouts due to snow damage averaged EUR 0.5 million over the past three decades, which constituted 7% of the total insurance compensations to forest owners (Lehtonen et al. 2016). In Finland and Norway alike, snow has become in the past decade the second most important abiotic forest disturbance agent after wind in terms of insurance compensation. At the same time, a decreasing trend was found by Díaz-Yáñez et al. (2016) in the frequency of snow damage events in Norway during the 1995–2014 period. In addition, any forest damage is capable of having indirect cascading effects triggering other natural disturbances, such as making the weakened trees susceptible to bark beetle attacks or increasing fuel loads, even in forests not typically prone to severe wildfires (Schroeder and Eidmann 1993; Klopcic et al. 2009; Machado Nunes Romeiro et al. 2022). From a short-term perspective, the direct effects of snow damage on the economic performance of individual forest estates and the whole forestry industry alike involve reduced timber quality and value, increased harvesting costs, and timber volume loss. The long-term effects of snow damage manifest themselves in reduced forest growth and lower harvest levels.

Reducing snow damage may require adaptations of the prevalent silvicultural practices. However, the existing body of knowledge on the occurrence in space and time of heavy snowfall, the physics of snow accumulation in tree crowns, and the mechanics of tree resistance to breakage appears to be limited. Furthermore, little scientific literature on adaptations to silvicultural practices exists in Norway and, unlike bark beetle attacks, in respect of which a preparedness plan and silvicultural recommendations were developed by the Norwegian Agriculture Agency (Landbruksdirektoratet 2021), nothing similar exists for snow damage. The prevalent practical advice aimed at reducing snow damage is heavy pre-commercial thinning, which is assumed to result in a greater single-tree stability through larger root systems and larger tree stem diameters, as well as more symmetric crowns (Solberg et al. 2022; Skogbrand Forsikringsselskap Gjensidig 2023). Adopting this practice, however, would mean reduced volumes and quality of timber (Mäkinen and Isomäki 2004; Mäkinen and Hein 2006; Moore et al. 2009; Šilinskas et al. 2020). An important direction of research is thus finding silvicultural practices balancing snow damage avoidance and revenue. This requires quantitative biological and economical assessment of any proposed silvicultural practices in terms of snow damage risk.

A risk does exist that the extent and frequency of wet snow damage in Nordic boreal forests might increase in the future as a consequence of the more frequent heavy precipitation events during milder and more humid winters (IPCC 2014). However, a clear-cut conclusion is hard to arrive at due to contradictory evidence and multiple parallel trends that might be at work. While forest disturbance due to snow and ice may be expected to become less frequent and severe globally under a future warmer and either drier or wetter climate, in Europe the prevailing expectations are of more snow damage under a warmer and wetter climate and of less damage under a warmer, but drier climate (Seidl et al. 2017). Subcontinent-scale studies of forest disturbance in a changing climate offer a more nuanced picture: e.g. while Patacca et al. (2023) expect global warming to make snow damage less frequent in Europe, they conclude that the uncertainty associated with the seasonality and severity of extreme wet snow precipitation events is likely to increase. This is in line with Masson-Delmotte et al. (2021), arguing that major interregional and interannual variability is to be expected in precipitation, and the fine-scale heterogeneous pattern of wet snowfall makes the overall direction of change in snow loads difficult to predict (Räisänen 2008, 2016; Strasser 2008; Lehtonen et al. 2016; Hawcroft et al. 2018; Ohba and Sugimoto 2020). Both O’Gorman (2014) and Räisänen (2016) point out that the future trend in northern Europe towards higher mean annual snowfall in colder areas and lower snowfall in warmer areas is more clear and more systematic than the trend in the intensity of extreme daily snowfall, which may remain virtually unchanged over a broad range of climatological monthly air temperatures. Quante et al. (2021) emphasize that extreme snowfall might become less frequent, but hardly less intense in the future in Western Europe. Räisänen (2016) concludes that a greater share of annual snowfall may in the future occur in periods with close-to-zero air temperatures, leading to higher snow accumulation in the forest canopy.

In their study of future snow loads in the built environment in the Nordic countries, Larsson Ivanov et al. (2022) found that the inland and higher-elevation areas on the east side of the Scandinavian mountains may towards the mid-twenty-first century experience an increase in extreme snow loads in all future climate scenarios (the increase is stronger in low-emission scenarios), while the coastal areas may experience a significant decrease. We note that their predicted spatial pattern of higher extreme snowfall is in strong agreement with the forest snow damage observations made in southern Norway over the past decade (Søgaard et al. 2017).

Snow accumulation in tree crowns is controlled by not only the amount of snowfall, air temperature, and relative humidity, but also wind speed. High risk of snow accretion is associated with air temperatures between −3°C and +1°C and wind speeds below 6 m s⁻¹; wind speed is negatively correlated with any type of snow load accumulation due to the fact that higher wind contributes to snow removal from tree crowns (Lehtonen et al. 2014). The effect of climate change on future windiness in northern Europe remains uncertain: Gregow et al. (2020) found little consensus in the literature apart from the expectation that autumn wind speed may slightly increase and that westerly winds may become more common due to the eastward extension of cyclone tracks. This prognosis is shared by a number of other recent studies (Groenemeijer et al. 2016; Kjellström et al. 2018; Ruosteenoja et al. 2019). Feser et al. (2015) note that while the number of storms over northwestern Europe may remain unchanged, their intensities may, similarly to those of snowfall, increase. Future snow accumulation in forest canopies will, thus, be to a greater extent controlled by a combination of mean and extreme snowfall levels and surface air temperature, rather than windiness.

Among the Nordic countries, the future trends in snowfall and wind appear to be studied most extensively for Finland and, to a lesser degree, Sweden. Regarding Finland, the general outlook is that forest canopy snow loads are expected to increase in the east and north of Finland and decrease in the south and west, both along the coast and in the interior (Venäläinen et al. 2020). Lehtonen et al. (2016) offer additional details by breaking down the total crown snow load into components (rime, dry snow, wet snow, frozen snow, see Lehtonen et al. (2014)): in the upcoming two decades, only small changes are predicted with a 10–30% increase in rime, wet snow, and frozen snow loads in northern Finland and no noticeable change elsewhere.

In Sweden, Larsson Ivanov et al. (2022) found similarly large regional differences in the generic case of future snow loads on the built environment; south Sweden and the east coast will experience a decrease in snow loads, while an increase is expected in the north and on the lee (eastern) side of the Scandinavian mountains.

The existing literature indicates that no generic prognosis can be made regarding the future severity, spatial distribution, and frequency of heavy snow precipitation on a subcontinental (Europe), or even regional (the Nordics) scale. The outlook for forest snow damage has necessarily to be nuanced depending on local climate, topography, and forest composition.

Modeling approaches to snow breakage in conifers

Although several bioeconomic scenario analysis models exist in Norway, these lack the capability to predict the snow damage risk. Output from a snow damage risk model should include both probability of damage and its effect (e.g. share of trees damaged) based on the forest structure and/or single-tree variables used in the scenario model, taking into account climatic variables. A recent review of forest disturbance models in Machado Nunes Romeiro et al. (2022) found very few examples of models that could be used in this way. One of the few models developed specifically for Norwegian conditions (Díaz-Yáñez et al. 2019) is based on forest structure variables as well as topographic and climatic variables as input, but predicts the combined risk of wind and snow damage.

There are relatively few attempts to model snow damage risk to trees and branches using a mechanistic modeling approach rather than a statistical one. Some examples include Petty and Worrell (1981), Cannell and Morgan (1989), Peltola et al. (1999), and Nock et al. (2021). Most models of snow damage tend to be statistical models based on observed snow damage (e.g. Valinger and Fridman (1997, 1999), Jalkanen and Mattila (2000), Lehtonen et al. (2014), Díaz-Yáñez et al. (2019), and Suvanto et al. (2021)). Such statistical models are only of use if a large amount of snow damage data is available, along with the key input variables in the models in order to make predictions in other locations or into the future (see Suvanto et al. (2021), Nykänen et al. (1997), and Lehtonen et al. (2014)). They do not allow the investigation of damage to individual trees based on tree characteristics and the level of snowfall from a snow damage event as would be possible with a mechanistic model. Considering existing mechanistic models, the Cannell and Morgan (1989) and the Nock et al. (2021) models only deal with branch breakage, while the Peltola et al. (1999) and the Petty and Worrell (1981) models, although they do deal with snow damage to the whole tree, assume either a displaced crown due to the stem shape or an initial deflection due to the wind. No previous paper, to our knowledge, deals with a mechanistic understanding of snow damage to tree stems without including the effect of wind or initial crown displacement. In this paper, we present a simple mechanistic model that does not require any initial stem displacement and predicts only the snow load required to damage the stems of individual trees. It is very similar in philosophy and physics to the model by Petty and Worrell (1981), but it works on the principle that when the snow load reaches a certain level, the system becomes unstable and even the slightest displacement from any cause (e.g. wind or uneven snow distribution in the crown) results in the collapse of the tree. This means that it is possible to calculate the critical snow load for stem damage with no requirement to include a wind speed or to create an asymmetric distribution of snow within the crown.

In this paper, we assess the ability of this idealized mechanistic model (Section 2.3) based on functions in the ForestGALES model (Hale et al. 2015) to predict stem breakage in conifers at a single-tree level when exposed to heavy wet snow loads. ForestGALES is a hybrid mechanistic model designed to predict the critical wind speed (CWS) leading to tree stem breakage or uprooting at the single-tree or stand level and includes the option to add a snow load in the tree crown.

In November 2016 and January 2018, parts of south-eastern Norway experienced two heavy wet snowfall events. One of the affected areas was the Fritzøe Skoger forest estate located in both the hemiboreal and boreal zones, where high levels of snow accumulation occurred in the forest canopy in the altitude interval of 100–300 m above sea level (m.a.s.l.), resulting in extensive stem breakage, both higher in the tree crown and closer to the ground. Snow damage was mostly concentrated in Norway spruce (Picea abies), with Scots pine (Pinus sylvestris) representing only approximately one quarter of the damaged trees. No evidence exists, however, indicating that the heavy snowfall events also involved elevated wind speeds, making it impossible to directly apply ForestGALES to explicitly estimate the critical snow load (CSL), rather than CWS. In addition to the CSL, a probability assessment of snow breakage would require a reference (actual or estimated) snow load value that the given trees were exposed to, or alternatively a distribution of values characterizing the local “snow climate.” However, the snow precipitation and snow cover data available for the area of interest (Norwegian Meteorological Institute, Norwegian Water Resources and Energy Directorate, and Norwegian Public Roads Administration 2023) was found to be too coarse with a grid size of 1 km and did not adequately capture the observed effect of terrain elevation on snow accumulation. Furthermore, the data was presented as snow load on the ground surface, rather than accumulation of snow in a three-dimensional forest canopy.

We present in this paper (1) a high-resolution (grid size 500 m) numerical weather prediction-based snow accumulation model providing the value of snow accumulation per unit area of forest canopy (Section 2.2), and (2) a mechanistic model to estimate CSL based on a set of tree attributes typically available from a single-tree aerial laser scanning (ALS) or aerial photogrammetric forest inventory (Section 2.3). Thereafter, we compare the predicted snow breakage probabilities to damage observations inferred from the bitemporal forest inventory data collected in 2014 and 2020 and assess the discrimination ability and predictive power of the proposed combination of the snow accumulation and CSL models.

Materials and methods

Study area and single-tree forest inventory data

The Fritzøe Skoger forest is a 73,000-ha private forest estate, where over 60% of the area is productive forest, of which more than 50% is dominated by Norway spruce, according to the forest management plan (Figure 1). For several decades, the primary regeneration method was planting of pure spruce stands, while in recent years this has transformed into establishing a greater share of mixed stands of Norway spruce and Scots pine. In recent decades, regular thinning has been minimized and stand density has mainly been controlled by pre-commercial thinning. The estate lies in a hilly terrain with elevations ranging from sea level up to 600 m.a.s.l. Most of the forest area has soils originating from shallow and discontinuous till sediments. At elevations lower than 100 m.a.s.l., the dominant conifer species is Norway spruce (>85%); the proportion of Scots pine increases to close to 25% at mid-elevations and becomes dominant (>50%) starting at 600 m.a.s.l.

Figure 1. Overview map of the study area.

In 2014 (prior to the snow damage events) and 2020 (after the snow damage events), single-tree forest inventories were conducted by Foran Sverige AB, a forestry mapping and consulting services provider, covering the entire Fritzøe Skoger forest estate, based on two ALS and aerial photography (three natural color and one near infrared bands) campaigns with a nominal point density of 8 and 10 points m⁻², respectively. The ALS point cloud had a horizontal and vertical accuracy of 0.21 and 0.28 m, respectively, based on measurements of 7 control profiles; the orthophoto mosaic resulting from the aerial photography campaign had a resolution of 0.25 m and a horizontal accuracy of 0.07 m, based on measurements of 7 ground control points (both accuracies reported for the 2020 aerial surveys). The two datasets contained such single-tree variables as treetop coordinates in ETRS89 UTM 32N, ALS tree height $h$, tree crown radius $r$ derived by tree crown segmentation, both in m. Treetops were identified by finding local maxima in the ALS point cloud that were at least 2 m higher than their surroundings, and their locations relative to neighboring treetops were used to segment the forest canopy into tree crown polygons and estimate crown radii. Diameter at breast height $d$ in cm estimated by control plot-calibrated allometric models as a function of tree species, $r$, and $h$. Tree species classification from the pre-damage 2014 inventory was used, based on a combination of spectral, such as NDVI (Tucker 1979), and crown geometry properties from a total of 4928 field-identified trees (3314 Norway spruces, 637 Scots pines, 977 deciduous trees) across 152 control plots.

Specifically, the following tree diameter functions were applied.

For Scots pine: (1) $d_B=0.0235796+5.16756r+1.67955h-0.0286432h^2$(1) For Norway spruce: (2) $d_B=0.208008+4.524r+1.67955h-0.0154905h^2$(2) where $d_B$ is diameter at breast height with bark in cm.

The tree diameter functions above were calibrated using a total of 379 control plots established randomly in 45 calibration stands. ALS-based estimates of diameter at breast height $d$ were validated against field measurements and aggregated on a calibration stand level $(R^2=0.65,n=45)$.

Diameter at breast height $d$ was subsequently estimated with correction for double bark thickness as follows: (3) $d=d_B-0.1B$(3) where $B$ is double bark thickness in cm, estimated by the equations below.

For Scots pine (Brantseg 1969): (4) $B=2.9571+1.1499d_B-0.7304{\displaystyle \frac{d_B}{h}}$(4) For Norway spruce (Vestjordet 1967): (5) $B=-0.34+0.831648d_B-0.002832d_B^2-0.010112h^2+0.700203{\displaystyle \frac{d_B^2}{h^2}}$(5) To identify trees that suffered stem breakage as a result of the heavy snowfall in 2016 and 2018, we matched identical trees in the 2014 and 2020 datasets by finding the Euclidean distance between pairs of treetop coordinates, selecting nearest neighbors, and excluding those separated by more than 0.25 m. Heavy accumulation of snow in the tree crowns was the primary cause of stem breakage in the study area in the past decade with no reports on windstorms, forest fires, or other abiotic disturbances causing reduction in treetop heights between 2014 and 2020. The matched dataset was dominated by Norway spruce (63% of the tree count), followed by birch (Betula spp.; 20%), and Scots pine (16%), as classified in the pre-damage 2014 inventory. We filtered out birch and applied a tree height constraint by excluding trees shorter than 10 m and taller than 30 m as measured in 2014. The reason for excluding younger and thus shorter trees was their presumed better ability to withstand high snow loads due to a lower modulus of elasticity and increased flexibility, while trees taller than 30 m were excluded to remove ALS returns potentially misclassified as treetops as, according to the Norwegian forest resource map SR16 (NIBIO 2022), the 99th percentile of dominant tree height in the study area was 29.3 m for Norway spruce and 27.3 m for Scots pine.

A given tree was assumed to have suffered snow breakage if its height as measured in 2014 and 2020 recorded a decrease of ≥0, 1, 2, or 3 m (${\mathrm{\Delta }}_h$). The four thresholds were applied to the single-tree dataset stratified by tree height in 2014 ($h_{2014}$) into three subsets $h_{2014}$ ≥ 10, 15, and 20 m for Norway spruce and Scots pine separately to obtain 24 strata. To balance the individual strata, random samples equal in size to the number of damaged trees were taken from the individual no-damage subsets over 100 iterations. Table 1 presents the summary statistics of the single-tree stratum defined by $h_{2014}$ ≥ 10 m as estimated in the 2014 inventory and Figure 2 shows characteristic examples of snow breakage of various severity in Norway spruce and Scots pine.

Figure 2. Examples of snow breakage of various degrees of severity in (a) Norway spruce and (b) Scots pine. (Photos by Peter Zubkov.)

Table 1. Summary statistics of the 2014 single-tree inventory stratum defined by $h_{2014}$ ≥ 10 m ($n=2,166,717$) by tree species.

Single-tree variable	Mean	Range	Standard deviation
Norway spruce ($n=1,825,651$)
d, cm	20.9	7.7–52.9	5.7
h, m	15.6	10–30	3.9
r, m	1.9	0.15–6	0.5
Taper (cm/cm)	74.5	40.1–130	6.8
Scots pine ($n=341,066$)
d, cm	29.3	14.9–54	4.9
h, m	14.4	10–30	3.2
r, m	2.2	0.15–5.8	0.6
Taper (cm/cm)	49.1	28–100	6.2

High-resolution snow accumulation model

Due to the hilly terrain and the lack of representative weather stations in the near vicinity of the study area, the meteorological dataset used in this study was obtained through numerical weather prediction (NWP). The data was prepared with the Weather Research and Forecasting (WRF) model version WRF-ARW 4.1.2 (Skamarock et al. 2021). It is a state-of-the-art mesoscale atmospheric model, used both for operational forecasting and for simulation of historic weather (i.e. hindcasting). Short-term hindcast simulations were produced by simulating the weather conditions during the entire months of November 2016 and January 2018 for a domain sized 43.5 by 67.5 km containing the Fritzøe Skoger forest estate at a 500 m horizontal grid spacing and at 51 vertical model levels throughout the atmosphere. All simulations were initiated and forced on the lateral boundaries by the global re-analysis ECMWF–ERA5 dataset (Hersbach 2016). In addition to producing standard weather variables, such as wind speed, wind direction, temperature, humidity and precipitation amount, the NWP model outputs detailed information about the precipitation properties, such as mass concentration and liquid water fraction (LWF) of snow and precipitation particles. The level of detail depends on the parameterization scheme of cloud and precipitation processes. In this study, we used the Thompson-Eidhammer aerosol-aware microphysics scheme (Thompson and Eidhammer 2014), with improvements to the melting layer parameterization (Iversen et al. 2021). Further information regarding the model setup is available in Nygaard et al. (2022).

Wet snow accumulation was calculated as a combination of horizontal and vertical precipitation fluxes for an idealized tree assuming a conic crown with a projected area of 1 m² in the horizontal and vertical plane (Lehtonen et al. 2014). Hourly time series of precipitation, wind, temperature, and humidity were extracted from the WRF hindcast, and the accumulated wet snow was calculated for each time step. The time-dependent snow accumulation follows that of Gregow et al. (2008), however our implementation of the accretion model defines wet snow as precipitation with a LWF in the range of 0–0.5. This is based on physical experiments on wet snow accumulation (Sakakibara et al. 2007; Hefny et al. 2009) which show that the adhesion of snow to various surfaces is restricted to a certain range of LWF values.

The snow loss component of the model follows Gregow et al. (2008) and involves melting and wind erosion. The melting rate increases with temperature in the range of 0 to +2.3 °C, with all snow melting away above the upper threshold. A cubic function was applied to account for snow loss due to wind erosion.

Mechanistic critical snow load modeling and snow damage probability estimation

To calculate the critical snow load leading to tree stem breakage, we used a simplification of the deflection loading factor (DLF) function in the ForestGALES model that calculates the additional moment created by a tree crown displaced by the accumulated snow. The form of the DLF function is as follows (Locatelli et al. 2022): (6) $DLF={\displaystyle \frac{1}{1-\left[y\left({\displaystyle \frac{l}{2}}\right)\times \left(W_{crown}+W_{snow}\right)+y\left({\displaystyle \frac{3}{4h}}\right)\times W_{stem}\right]\times {\displaystyle \frac{g}{B_w}}},}$(6) where $l$ is the crown depth in m, $W_{snow}$ is the weight of the snow accumulation in the tree crown in kg, $W_{crown}$ is the crown weight in kg, $h$ is tree height in m, $W_{stem}$ is the stem weight of the tree in kg, $g$ is the gravitational acceleration constant (9.81 m s⁻²), $B_w$ is the bending moment exerted by the wind on the tree in Nm, and $y\left(x\right)$ is a function describing the displacement of the tree at different points $x$ along the tree stem (Gardiner 1989; Quine et al. 2021). DLF does not depend on the wind speed because it is a multiplier applied to the moment due to the wind and is only a function of the characteristics of the tree (stem shape, stem density, stem stiffness, crown size, crown density). A value of DLF = 1 would mean that there is no contribution of the crown and stem to the overturning moment. The average value of DLF with no snow loading is 1.136 (Hale et al. 2015) so that the crown and stem displacement by the wind contributes an additional 13.6% to the overturning moment.

$W_{crown}$ and $l$ are estimated by the equations below: (7) $W_{crown}=\rho \pi r^2{\displaystyle \frac{l}{3},}$(7) (8) $l=b_0+b_1\cdot h$(8) where $\rho$ is the crown density in kg m⁻³, $r$ is the tree crown radius in m, $b_0$ and $b_1$ are the regression coefficients.

The function $y\left(x\right)$ in Equation (6) is given by the following equation, which assumes the stem shape follows a power to the 0.6 function with distance from the treetop (Quine et al. 2021): (9) $\begin{array}{rl}& y\left(x\right)=\frac{F\cdot h^3}{MOE\cdot I_h}[2.5\left(\frac{x}{h}\right)-4.17{\left(\frac{x}{h}\right)}^{0.6}-0.71q\left(\frac{x}{h}\right)\\ & -1.79q{\left(\frac{x}{h}\right)}^{-0.4}+2.5q+1.67],\end{array}$(9) where $x$ is the location along the tree stem measured from the treetop in m, $F$ is the force of the wind, $MOE$ is the modulus of elasticity for the given tree species, $q$ is the ratio of the distance from the treetop to the point where the wind force $F$ is applied to the tree height $h$, and $I_h$ is the area moment of inertia at the base of the tree, approximated by: (10) $I_h\cong {\displaystyle \frac{\pi d^4}{64},}$(10) where $d$ is the tree diameter at breast height in m.

The contribution from the stem weight in Equation (6) is marginal and can be ignored. The bending moment due to the wind $B_w$ can be approximated by $F\cdot h$ and we assume that the wind acts at the center of the tree crown, thus: (11) $x=\frac{l}{2},$(11) (12) $q={\displaystyle \frac{l}{2h}.}$(12) Equation (9) can thus be simplified as follows: (13) $\begin{array}{r}y\left(x\right)={\displaystyle \frac{h^2}{MOE\cdot I_h}\left[5.0\left({\displaystyle \frac{l}{2h}}\right)-5.96{\left({\displaystyle \frac{l}{2h}}\right)}^{0.6}-0.71{\left({\displaystyle \frac{l}{2h}}\right)}^2+1.67\right].}\end{array}$(13) Equation (6) can accordingly be rewritten as: $\begin{array}{r}DLF={\displaystyle \frac{1}{1-\left[{\displaystyle \frac{h^{2}g}{MOE\cdot I_h}\left[5.0\left({\displaystyle \frac{l}{2h}}\right)-5.96{\left({\displaystyle \frac{l}{2h}}\right)}^{0.6}-0.71{\left({\displaystyle \frac{l}{2h}}\right)}^2+1.67\right]\cdot \left(W_{crown}+W_{snow}\right)}\right]}.\left(14\right)}\end{array}$As the denominator in Equation (14) approaches zero, the DLF approaches infinity and the tree collapses under the combined weight of its crown and accumulated snow without any additional wind loading. The CSL in kg can thus be derived as: (15) $\begin{array}{rl}& CSL=\\ & \frac{1}{\frac{64\cdot h^{2}g}{\pi \cdot MOE\cdot d^4}\left[5.0\left(\frac{l}{2h}\right)-5.96{\left(\frac{l}{2h}\right)}^{0.6}-0.71{\left(\frac{l}{2h}\right)}^2+1.67\right]}\\ & -W_{crown}.\end{array}$(15) Model parameter values used in Equations (7) to (15) are given in Table 2.

Table 2. Model parameter values used to calculate critical snow load (CSL) from Equations (7, 8 and 15) (from Locatelli et al. (2022)).

Model parameter	Unit	Tree species
		Scots pine	Norway spruce
$MOE$	N m^−2	$7.3\times {10}^9$	$6.3\times {10}^9$
$b_0$	m	1.3424	0.4206
$b_1$	unitless	0.3156	0.4368
$\rho$	kg m^−3	2.5	2.5

Probability $p$ of tree stem breakage when exposed to a snow load above the CSL is assumed to be a logistic function of the form: (16) $p={\displaystyle \frac{1}{1+e^{{\displaystyle \frac{-{\mathrm{\Delta }}_{SL}}{1000}}}},}$(16) where ${\mathrm{\Delta }}_{SL}$ is the difference between the modeled snow load in kg on the given tree (Section 2.2) and the CSL given by Equation (15). Logistic functions are regularly used in assessing the probability of an event varying between 0 and 1 in a smooth sigmoid shape (e.g. Hosmer et al. (2013), Hale et al. (2015)).

For spatial analysis and visualization of the snow accumulation models and snow breakage predictions, we used ArcGIS Pro version 3.1.

Model discrimination ability and accuracy assessment

To assess the discrimination ability of the CSL model, we calculated the value of the area under the receiver operating characteristic (ROC) curve (AUC) for each of the strata. To make an accuracy assessment of the binary snow breakage predictions, we chose to optimize the cutoff value of $p_{threshold}$ such as to maximize the sum of sensitivity $S_1$ and specificity $S_2$ in the resulting two-by-two confusion matrix, rather than apply a fixed classification threshold of $p_{threshold}=0.5$. For statistical analysis, we used R version 4.2.2 (R Core Team 2022), including the cutpointr package (Thiele and Hirschfeld 2021) for threshold optimization and accuracy assessment, and the ggplot2 package (Wickham 2009) for visualization. The accuracy measures calculated are $S_1$, $S_2$, and F1 score and are given as follows: (17) $S_1={\displaystyle \frac{TP}{TP+FN},}$(17) (18) $S_2={\displaystyle \frac{TN}{FP+TN},}$(18) (19) $F1={\displaystyle \frac{2\cdot TP}{2\cdot TP+FP+FN},}$(19) where $TP$, $TN$, $FP$, and $FN$ are, respectively, the true positives, true negatives, false positives, and false negatives in the corresponding binary confusion matrix.

The reported accuracy measures and $p_{threshold}$ were averaged over 100 resampling iterations on a balanced dataset using a bagging procedure.

Results and discussion

Snow accumulation model

The modeled snow loads as shown in Figure 3 tend to capture well the spatial distribution with higher values in the same elevation intervals as the observed forest snow damage. At higher elevations, the model mainly predicts precipitation as dry snow while in low-lying areas the precipitation is mainly modeled as sleet or rain due to a high LWF. We were prevented from quantifying the accuracy of the snow accumulation model by the absence of weather stations in the study area representing the elevation intervals of interest; the only alternative was comparing estimated snow load values with the spatial distribution of tree height decrease between 2014 and 2020 interpreted as snow breakage.

Figure 3. (a) Maximum snow accumulation maps of the study area in November 2016 and (b) January 2018 (numbers relate to terrain elevation contour intervals), (c) maximum mean snow accumulation (shading indicates +/− 2σ) by elevation in the respective time periods, and (d) snow-damaged trees in 2018 (for ${\mathrm{\Delta }}_h=1m$) by elevation interval and tree species.

We found two patterns in the snow accumulation in different elevation intervals in November 2016 and January 2018 (Figure 3(a–c)). In November 2016, the highest accumulation was modeled in the 100–200 m.a.s.l. interval (mean snow accumulation 86 kg m⁻² with close to half these values in the adjacent 100 m elevation intervals above and below this range). In January 2018, the highest snow accumulation was modeled at 117 kg m⁻² in the same 100–200 m.a.s.l. interval. However, the 2016 maximum was exceeded in a wider belt between 0 and 300 m.a.s.l. with no difference in snow accumulation at higher elevations (Figure 3(c)). Extensive snow breakage in both Norway spruce and Scots pine was found in the 100–300 m.a.s.l. interval, with less damage to Scots pine at elevations below 100 m.a.s.l. where this species is relatively less common and a gradual increase in the proportion of damaged Scots pine starting at 400 m.a.s.l. (Figure 3(d)). The distribution of the observed snow damage by elevation intervals reflects the distribution of the two tree species, with Norway spruce being the dominant species at lower elevations (Figure 1).

The spatial distribution of snow damage demonstrated species-dependent statistical significance ($p<0.05$), with clustering and dispersion patterns in the damage class, as indicated by the Getis-Ord Gi* statistic (Getis and Ord 1992): snow breakage was strongly clustered in the 100–200 m.a.s.l. interval in Norway spruce and in the 0–300 m.a.s.l. interval in Scots pine. A dispersion pattern was found above 400 m.a.s.l., mostly in snow-damaged Scots pine. The no-damage class showed clustering above 300 m.a.s.l. and dispersion below 200 m.a.s.l. in Norway spruce, with a similar, though less pronounced, pattern in Scots pine (Figure 4(a and b)).

Figure 4. Snow accumulation (maximum of November 2016 and January 2018) and Getis-Ord Gi*-based spatial clustering of damage and no-damage observations in (a) Norway spruce and (b) Scots pine. Clustering and dispersion by snow accumulation by tree species in the (c) damage and (d) no-damage class (${\mathrm{\Delta }}_h$≥1 m).

Both clustering and dispersion correlated with snow accumulation: in the damage class, significant clustering ($p<0.1$) was found at median snow accumulation values higher than 115 kg m⁻² in both Norway spruce and Scots pine. Conversely, in the no-damage class clustering was found in Norway spruce at median snow accumulations below 40 kg m⁻², while dispersion was observed at values above 100 kg m⁻² (Figure 4(c and d)).

It should be kept in mind that the snow accumulation model (Section 2.2) was built for a highly simplified conifer crown, taken outside of its spatial and social context, including such factors as forest edge and gap effects, sheltering by neighbors, and snow load sharing with neighboring trees. Furthermore, as noted in Lehtonen et al. (2014), the snow load is greatly affected by individual tree properties and adaptations, such as crown shape. For example, spruce trees at higher elevations tend to develop narrow, candle-shaped crowns that may accumulate higher snow loads than their lowland counterparts due to snow packing and the lower center of gravity preventing shedding of the snow load (Büsgen et al. 1929). This means that our model might underestimate snow accumulation at elevations above 400 m.a.s.l.

To the best of our knowledge, this study is the first ever attempt at predicting single-tree snow breakage using high-resolution NWP methods. In Suvanto et al. (2021), the same objective was dealt with by directly applying the snow accumulation model of Lehtonen et al. (2014) to the spatially more coarse ERA5 global reanalysis data and an “average” tree representing an entire forest stand, while in Lehtonen et al. (2014) no predictions of the risk of tree breakage were made based on the NWP models (such as the meso-scale HIRLAM and the global IFS by ECMWF).

Critical snow loads and snow damage predictions

The CSL model (Section 2.3) combined with the snow accumulation model (Section 2.2) demonstrated satisfactory accuracy of snow damage predictions, with the highest classification accuracy achieved in the moderately restrictive tree height $h_{2014}$ and height decrease ${\mathrm{\Delta }}_h$ strata. Four combinations of the criteria giving the highest and most balanced (in terms of $S_1$ vs. $S_2$) classification performance with their accuracy metrics for Norway spruce and Scots pine are presented in Table 3. Applying a height decrease ${\mathrm{\Delta }}_h$ of at least 2 m to the subset of trees higher than 10 m in 2014 gave the best discrimination ability in terms of AUC (see Figure 5(b) for the ROC curves). The benefits of applying a more restrictive ${\mathrm{\Delta }}_h$ were most expressed in terms of sensitivity and negligible in terms of specificity, irrespective of tree species (Table 3). Applying the model to the strata defined by a higher $h_{2014}$ led to a major drop in $S_2$ across the examined strata with a slight to moderate increase in $S_1$ and negligible effect on AUC and F1 score. We conclude that the proposed model performed best on those conifers that have passed the juvenile development stage (higher than 10 m) at the time of exposure to a heavy snow load and suffered moderately severe stem breakage (height decrease greater than 2 m by the end of the period) as a consequence. In the following, we present and discuss our findings regarding this stratum, unless otherwise indicated.

Figure 5. (a) Density plots of damage probability ($p$) for the no-damage and damage classes by tree species for the entire dataset, showing the optimal $p_{threshold}$ values for the best-performing stratum ($h_{2014}$ ≥ 10 m, ${\mathrm{\Delta }}_h$ ≥ 2 m). (b) ROC curves, incl. $p_{threshold}$, and (c) sensitivity vs. specificity plot for the same stratum for Norway spruce (NS) and Scots pine (SP). Note the different scaling of the vertical axis, reflecting the class imbalance.

Table 3. Mean classification accuracy metrics (after 100 resampling iterations) of the snow breakage model for the four best performing combinations of minimum tree height in 2014 h₂₀₁₄ and tree height decrease Δ_h between 2014 and 2020, indicating the stratum size $n$ of the damage class and the value of $p_{threshold}$ maximizing the F1 score (NS, Norway spruce; SP, Scots pine). Best-performing combination highlighted with bold.

${\mathrm{\Delta }}_h$	Metric	$h_{2014}$
		10 m		15 m
		NS	SP	NS	SP
1 m	$n$	11,497	3,465	7,882	1,670
	AUC	0.71	0.73	0.71	0.69
	$p_{threshold}$	0.55	0.28	0.55	0.29
	$S_1$	0.68	0.76	0.71	0.84
	$S_2$	0.72	0.63	0.68	0.52
	F1	0.69	0.69	0.7	0.68
2 m	$n$	6,214	1,567	4,363	852
	AUC	0.72	0.77	0.72	0.72
	$p_{threshold}$	0.55	0.29	0.55	0.3
	$S_1$	0.72	0.81	0.73	0.87
	$S_2$	0.71	0.66	0.68	0.53
	F1	0.71	0.72	0.71	0.7

It should be borne in mind that the value of ${\mathrm{\Delta }}_h$ is only a proxy for the actual loss of tree height as a result of snow breakage and is not corrected for height growth between 2014 and the time of breakage (late 2016 or early 2018). Depending on the site index and assuming a time interval of three vegetation seasons between 2014 and early 2018, the actual point of stem breakage corresponding to a ${\mathrm{\Delta }}_h$ of 1 m can reasonably be expected to be located 2–3 m below the treetop (Kvaalen et al. 2015; Solberg et al. 2019). Light snow damage confined to <2 m below the treetop is thus outside the scope of this study and may be represented by a close to zero or positive ${\mathrm{\Delta }}_h$.

Two additional sources of uncertainty associated with the single-tree inventory data used are (1) accuracy of treetop identification and tree radius estimation and (2) reliability of the subsequent diameter at breast height estimates. Errors at the treetop identification stage may translate into (a) incorrectly matched treetops between 2014 and 2020, resulting in an imperfect reference dataset, and (b) incorrect diameter at breast height estimates, because crown size is part of the diameter functions (Equation (1 and 2)) and this could lead to incorrect CSL values.

We found that in the best-performing stratum, the damage class had a higher median value of $p$ compared to the no-damage class – 0.65 vs. 0.5 in Norway spruce and 0.4 vs. 0.25 in Scots pine (Figure 6(c)). Across the two classes, Norway spruce had lower CSL values compared to Scots pine; this can traced back primarily to the difference in stem taper between the two species (Table 1), with stem taper values in Norway spruce approaching the zone of risk of snow breakage, i.e. 75 and higher (Petty and Worrell 1981). There was, however, only a slight difference in the CSL between the damage and no-damage classes in Norway spruce with the direction of the difference inconsistent with what was expected. Scots pine had, in contrast, a larger difference in the CSL of 2000kg vs. 1500 kg in the damage and no-damage classes, respectively, and a correct direction of the relationship (Figure 6(b)). For the snow accumulation values, we found a major difference (1200 kg vs. 500 kg in the damage and no-damage classes, respectively) in Norway spruce and a smaller difference in Scots pine (Figure 6(a)). The probability values ($p$) were also higher in the elevation interval 0–300 m.a.s.l., which is the one with the highest snow accumulation levels, especially in Norway spruce (Figure 7). The $p$ values correlated stronger with tree snow loads rather than the CSL (Table 4), suggesting that the occurrence of snow breakage was predominantly controlled by snow accumulation in both tree species. This is unsurprising considering the finding in Päätalo et al. (1999) that snow accumulation above 60 kg m⁻² makes conifers susceptible to stem damage and uprooting and the modeled snow accumulation greater than 60 kg m⁻² across the elevation intervals of interest as shown in Figure 3(c). A snow load of 60 kg m⁻² can be conceptually interpreted as 60 cm of fresh snowfall, however Nykänen et al. (1997) offer a word of caution indicating the oversimplified nature of this conversion ignoring LWF in the snow and side wind; besides, such conversion would imply a strong assumption of zero pre-existing snow load in the forest canopy.

Figure 6. Boxplots of (a) snow accumulation in the tree crowns, (b) CSL, and (c) $p$ for the no-damage and damage classes by tree species ($h_{2014}$ ≥ 10 m, ${\mathrm{\Delta }}_h$ ≥ 2 m) (NS, Norway spruce; SP, Scots pine).

Figure 7. Boxplot of $p$ by elevation, class, and tree species ($h_{2014}$ ≥ 10 m, ${\mathrm{\Delta }}_h$ ≥ 2 m) (NS, Norway spruce; SP, Scots pine).

Table 4. Correlation matrix of the damage probability ($p$) versus tree crown snow accumulation (SL) and CSL by tree species (NS, Norway spruce; SP, Scots pine) ($h_{2014}$ ≥ 10 m).

		Pearson correlation coefficient
		SL	CSL
$p$	NS	0.78	0.05
	SP	0.69	−0.12

We note that the detected snow breakage cannot be reliably attributed to a specific snow accumulation event, i.e. either November 2016 or January 2018, because of the unavailability of remotely sensed data of sufficiently high spatial resolution collected between 2014 and 2020. The observed snow damage pattern can arguably be interpreted a combined effect of the two events, i.e. that the trees were first weakened by sustained bending in 2016 (Petty and Worrell 1981) and suffered extensive stem breakage when exposed to even greater loads in 2018 (Nykänen et al. 1997).

Despite its simplicity, the CSL model (Section 2.3) performed reasonably well and demonstrated a satisfactory discrimination ability with a clear, though imperfect separation of classes (Figure 5(a)). Similarly to the snow accumulation model in Section 2.2, we assumed a highly idealized conical tree crown shape without regard to the effect of tree species and elevation-dependent adaptations, such as candle-shaped vs. broad crowns in Norway spruce at higher vs. lower elevations, respectively (Lehtonen et al. 2014). We further assumed that the snow load acts symmetrically at the tree crown’s estimated center point; this may not hold in the situations of asymmetric snow accumulation due to wind or the shielding effect of the neighboring trees leading to an upward shift along the tree stem of the point where the snow load acts. The latter would cause stem breakage close to the treetop, a situation our model is not expressly designed to handle.

Model performance in the context of other existing snow damage models

In Suvanto et al. (2021), the same winter of 2017–2018 marked by heavy snowfall was used to validate their model’s snow damage risk predictions for Finland, and similarly to us, they concluded that crown snow load was the main driver of snow breakage during that winter. Our snow damage model achieved a classification performance that is comparable to that reported in Valinger and Fridman (1999) and Suvanto et al. (2019), both following a statistical approach, but calculated for individual trees, rather than stand or plot level values and using a simpler, physics-based approach. A clear advantage of our proposed model is that it requires as inputs only the most basic ALS or photogrammetric single-tree characteristics and reference snow load values (modeled or historical), making it suitable for applications in forest management planning and scenario analysis.

Similarly to Peltola et al. (1999), we saw that snow accumulations above 60 kg m⁻² made both Norway spruce and Scots pine susceptible to stem breakage, however without any additional wind loading – since wind contributes to shedding of the crown snow accumulation (Nykänen et al. 1997; Lehtonen et al. 2014) – and irrespective of the tree’s location relative to the forest edge.

The performance of our mechanistic model is in line with the earlier finding reported in e.g. Nykänen et al. (1997) that the main biotic factor controlling the CSL in conifers is stem diameter, followed by tree height, MOE, and crown length (see Equation 15), with the former two represented by taper, earlier identified as the single most important factor of tree stability under a snow load (Petty and Worrell 1981).

Further work

Obvious directions for further research into the risk of snow damage in Nordic conifer stands include:

1. Improvements to snow accumulation modeling, including making the model aware of the spatial and social status of trees (e.g. forest edge effects, sheltering by and snow load sharing between neighboring trees, and local topography) and making snow accumulation dependent on the tree species and crown shape.

2. Validating the snow accumulation model based on actual observations of snowfall composition and snow accretion to forest canopies.

3. Improvements to the CSL model, such as implementing more advanced stem taper functions, differentiating between tree species and their adaptations in terms of crown shapes and sizes, and importantly, adding functionality to predict the CSL at different heights along the tree stem (see Nykänen et al. (1997)).

Conclusion

In this paper, we present (1) a high-resolution NWP-based snow accumulation model for forest canopies and (2) a mechanistic CSL model for predicting the risk of snow damage in common Nordic conifer tree species and we validate the predictions using snow damage observations collected from single-tree ALS surveys. Despite the simplified and idealized formulation of the CSL model, we found the combined snow damage model to be both intuitive and capable of discriminating between damage and no-damage classes in the event of moderately severe stem breakage and of correctly representing the relative contribution of CSL and snow accumulation to the probability of snow damage. Our accuracy assessment indicates that the proposed model is well suited for applications in scenario analysis and snow damage risk assessment.

Acknowledgements

The authors are grateful to Christer Sandum and Fritzøe Skoger AS for making available the ALS and forest inventory data and to Eva Solbjørg Flo Heggem for pre-processing the forest inventory data.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Data availability statement

The data that support the findings of this study are available from the corresponding author, PZ, upon reasonable request. Restrictions may apply to the availability of the source data, which were used for this study and are available from the corresponding author with the permission of Fritzøe Skoger AS.

Additional information

Funding

This work was supported by the Research Council of Norway under Grant 301745 as part of the MARCSMAN (Maximizing the Resilience and Carbon Sequestration in Managed Norway Spruce Forests) project.