Journal topic
Hydrol. Earth Syst. Sci., 24, 2545–2560, 2020
https://doi.org/10.5194/hess-24-2545-2020
Hydrol. Earth Syst. Sci., 24, 2545–2560, 2020
https://doi.org/10.5194/hess-24-2545-2020

Research article 15 May 2020

Research article | 15 May 2020

# Snow processes in mountain forests: interception modeling for coarse-scale applications

Snow processes in mountain forests: interception modeling for coarse-scale applications
Nora Helbig1, David Moeser2, Michaela Teich3,4, Laure Vincent5, Yves Lejeune5, Jean-Emmanuel Sicart6, and Jean-Matthieu Monnet7 Nora Helbig et al.
• 1WSL Institute for Snow and Avalanche Research SLF, Davos, Switzerland
• 2USGS New Mexico Water Science Center, Albuquerque, NM, USA
• 3Austrian Federal Research Centre for Forests, Natural Hazards and Landscape (BFW), Innsbruck, Austria
• 4Department of Wildland Resources, Utah State University, Logan, UT, USA
• 5Université Grenoble Alpes, Université de Toulouse, Météo-France, CNRS, CNRM, Centre d'Études de la Neige, Grenoble, France
• 6Université Grenoble Alpes, CNRS, IRD, Grenoble INP, Institut des Géosciences de l’Environnement (IGE) – UMR 5001, 38000 Grenoble, France
• 7Université Grenoble Alpes, INRAE, LESSEM, 38402 St-Martin-d'Hères, France

Correspondence: Nora Helbig (norahelbig@gmail.com)

Abstract

Snow interception by the forest canopy controls the spatial heterogeneity of subcanopy snow accumulation leading to significant differences between forested and nonforested areas at a variety of scales. Snow intercepted by the forest canopy can also drastically change the surface albedo. As such, accurately modeling snow interception is of importance for various model applications such as hydrological, weather, and climate predictions. Due to difficulties in the direct measurements of snow interception, previous empirical snow interception models were developed at just the point scale. The lack of spatially extensive data sets has hindered the validation of snow interception models in different snow climates, forest types, and at various spatial scales and has reduced the accurate representation of snow interception in coarse-scale models. We present two novel empirical models for the spatial mean and one for the standard deviation of snow interception derived from an extensive snow interception data set collected in an evergreen coniferous forest in the Swiss Alps. Besides open-site snowfall, subgrid model input parameters include the standard deviation of the DSM (digital surface model) and/or the sky view factor, both of which can be easily precomputed. Validation of both models was performed with snow interception data sets acquired in geographically different locations under disparate weather conditions. Snow interception data sets from the Rocky Mountains, US, and the French Alps compared well to the modeled snow interception with a normalized root mean square error (NRMSE) for the spatial mean of ≤10 % for both models and NRMSE of the standard deviation of ≤13 %. Compared to a previous model for the spatial mean interception of snow water equivalent, the presented models show improved model performances. Our results indicate that the proposed snow interception models can be applied in coarse land surface model grid cells provided that a sufficiently fine-scale DSM is available to derive subgrid forest parameters.

1 Introduction

Snow interception is the amount of snow captured in a forest canopy. As much as 60 % of the cumulative snowfall may be retained in evergreen coniferous forests . In deciduous forests in the southern Andes as much as 24 % of the total annual snowfall may be retained . Due to the sublimation of intercepted snow, a large portion of this snow never reaches the ground , and the interplay of interception and sublimation creates significant below-forest heterogeneity in snow accumulation. estimated that 20 % of the seasonal snow cover in the Northern Hemisphere is located within forested areas. As such, the mass balance of solid precipitation in forested regions, characterized by strong spatial variability of snow accumulation, is a large contributor to the global water budget. Accurately modeling the spatial distribution of snow in forested regions is thus necessary for climate and water resource modeling over a variety of scales . Furthermore, intercepted snow can drastically change land surface albedo values in forested regions. Previous studies observed large albedo differences (a range of 30 %) between snow-free and snow-covered forest stands . Thus, in mountainous areas where forested and alpine regions coexist, accurate estimates of forest albedo play a key role in correctly modeling the surface energy balance. Due to the connectivity between interception and albedo, formulations of surface albedo over forested areas necessitate estimates of intercepted snow .

To date, direct snow interception measurements have only been retrieved from weighing trees. These measurements are limited to the point scale, are resource-intensive sampling, and only allow for the analysis of small- to medium-sized trees or tree elements . However, there are indirect techniques that allow for estimations of interception over larger spatial scales. Indirect measurements that compare snow accumulation between open and forest sites allow for a larger spatial sampling but may be affected by other forest snow processes, such as unloading of the intercepted snow. As such, sample timing of snowstorm conditions needs to be evaluated . Until recently, snow interception could not be characterized over length scales on the order of several tens of meters. However, at these scales snow interception can spatially vary due to canopy heterogeneity. The extensive data set of indirect snow interception measurements in evergreen coniferous forests (further referred to as coniferous forest) in eastern Switzerland collected by is likely the first data set that allows for a thorough spatial analysis of snow interception.

Several statistical models for forest interception of snow water equivalent (ISWE) have been suggested using a variety of canopy metrics and functional dependencies for the rate and amount of storm snowfall . Though these models have been demonstrated to perform well, they often rely on detailed forest canopy density and structure metrics that are either not readily available or cannot easily be upscaled, thus limiting functionality in models where the mean of model grid cells over several hundreds of meters to a few kilometers is required, which potentially reduces validity in large-scale modeling efforts.

Traditional forest metrics used to parameterize snow interception include leaf area index (LAI), canopy closure (CC), and canopy gap fraction (GF) or sky view. These are mainly derived from hemispheric photographs (HP) taken from the forest floor looking upwards. However, these indices can also be estimated from synthetic hemispheric photographs (SP). The SP mimic HP but are generated from aerial light detection and ranging (lidar) data. This requires the inversion of lidar to a ground perspective and conversion from a Cartesian to a polar coordinate system . Prior work has also used return density ratios of lidar, which is computationally faster but less accurate than SP . Canopy structure, or the position of a canopy element relative to the surrounding forest canopy, has also been used to model snow interception. However, as pointed out by , some forest structure metrics such as LAI and CC are highly cross correlated. Therefore, expanded on prior interception models (which mostly rely on the highly cross correlated traditional forest density parameters of LAI and CC) by introducing uncorrelated novel forest structure metrics. Their empirical interception model utilizes the total open area, mean distance to canopy, and CC. While the latter parameter was derived from SP , the first two parameters were directly computed by a digital surface model (DSM). The total open area is defined as the total open area in the canopy around a point, and the mean distance to the canopy defines how far away the edge of the canopy is from a point. Recently extended the mean distance to the canopy vertically by deriving it for 1 m horizontal slices that were normalized with the corresponding elevation above the ground.

Due to the difficulties in measuring snow interception, previous empirical snow interception models were not validated in different snow climates, forest types, or at varying spatial scales. During the Snow Model Intercomparison Project Phase 2 (SnowMIP2) 33 snow models were validated at individual forested and open sites and many models used the snow interception parameterization from . This interception model was one of the first that used canopy metrics (LAI and CC), although a snow interception model for larger scales also requires the greater canopy structure. Overall, SnowMIP2 showed that maximum snow accumulation predictions had large errors compared to observed values in most models, but the snow cover duration was well estimated. Furthermore, a universal best model could not be found because model performances at forest sites varied. This may explain why there is still no common ground with several snow-related variables in land surface models , which led to the current Earth System Model-Snow Model Intercomparison Project (ESM-SnowMIP) showing overall larger errors in simulated snow depth on forest sites than on open sites . Recently validated three snow interception models developed for coniferous forests with observed point snow interception values in a deciduous southern beech (Nothofagus) forest of the southern Andes. All three empirical models required recalibration, with the recalibrated model showing the overall best performance. Similarly, model simulations of largely overestimated the observed accumulated snow depth in a spruce forest at Col de Porte in the southeastern French Alps. They attribute this to errors in the processes linked to the snow interception model based on due to an underestimation of the melt of intercepted snow. In a maritime climate previous snow interception models also failed to accurately model snow interception . While successfully modeled snow interception in a maritime climate, their model consistently underestimated snow interception in a continental climate forest. Overall, this demonstrates the need for more robust parameterizations of the processes affecting snow under forest, which is an important challenge for global snow modeling.

When modeling at resolutions greater than the point scale, accurate implementation of forest snow processes necessitates not just the mean of a grid cell but the standard deviation within a grid cell or model domain. However, to our knowledge, the standard deviation of snow interception has not yet been quantified. In this paper, we propose empirical parameterizations for the spatial mean and standard deviation of snow depth interception (IHS and ${\mathit{\sigma }}_{{I}_{\mathrm{HS}}}$) derived from indirect interception measurements at sites with length scales on the order of several tens of meters. We analyzed an extensive data set consisting of several thousand interception measurements collected immediately after storm events in a discontinuous coniferous forest stand in the eastern Swiss Alps . From a lidar DSM with elevations z , we derived the following two canopy structure metrics: (1) the standard deviation of the DSM (σz) in order to represent the spatial heterogeneity of surface height in a forested model domain; and (2) the spatial mean sky view factor (Fsky), which roughly represents the spatial mean canopy openness but is derived here on the DSM from geometric quantities that describe the received radiative flux fraction emitted by another visible surface patch (i.e., canopy patches) . These two metrics were correlated to spatial means and standard deviation of the indirect interception measurements. We validated the novel models with new indirect snow interception measurements from one site located in the Rocky Mountains of northern Utah, US, and from one site located at Col de Porte in the southeastern French Alps.

2 Data

In this study we only used indirect snow depth interception measurements. Indirect snow interception data were obtained by comparing new snow depth accumulation on the ground between open and forest sites. As such, snow depth interception (further referred to as snow interception) leads to reduced snow depth on the ground at forest sites. This indirect measurement technique allows for a collection of snow interception data over a larger area and also for an investigation of spatial snow interception variability. We used three snow interception data sets as follows: one from the eastern Swiss Alps for the development of snow interception models and two data sets for the independent validation of the developed snow interception models. Of these, one data set was from the Rocky Mountains of northern Utah in the US and one from the southeastern French Alps. In each of the three data sets the snow interception was derived slightly differently, which is described in the following sections.

## 2.1 Eastern Swiss Alps

Indirect interception measurements were collected in seven discontinuous coniferous forest stands near Davos, Switzerland, at elevations between 1511 and 1900 m above sea level (a.s.l.) and primarily consisted of Norway spruce (Picea abies) (Fig. 1a). Mean annual air temperature in Davos (1594 m a.s.l.) is approximately 3.5 C and the average solid precipitation is 469 cm per year (climate normal 1981–2010, https://www.meteoswiss.admin.ch, last access: 12 May 2020). The field sites are maintained and operated by the Snow Hydrology research group of the WSL Institute for Snow and Avalanche Research SLF in Davos, Switzerland. The sites were chosen to limit the influence of slope and topographic shading while capturing as much diversity as possible in elevation, canopy density, and canopy structure (see canopy height models, CHMs, of two field sites in Fig. 2). All seven field sites were equipped in the same manner and consisted of 276 marked and georectified measurement points (about ±50 cm) over a 250 m2 surface area (yellow inlet in Fig. 1a corresponds to each yellow dot). Two nonforested reference sites (open field sites) (see blue dots in Fig. 1a) were equipped with 50 measurements points each to derive the average open-site snowfall (accumulated snowfall).

Figure 1Extent of lidar-derived canopy height model (CHM) with the locations of open field sites (blue points), forested field sites (yellow points), and a SNOTEL site (purple point) (a) close to Davos in the eastern Swiss Alps (∼90 km2; 46.78945 N, 9.79632 E), (b) in the Rocky Mountains of northern Utah, US (∼13 km2; 41.85046 N, 111.52751 W), and (c) in the southeastern French Alps at Col de Porte (∼0.01 km2; 45.29463 N, 5.76597 E). The yellow-framed inlets show the respective snow depth measurement setup at the forested field sites. Underlying orthophotos were provided for the French site by IGN (France) and for the Swiss site by Swisstopo (JA100118). For the site in the US © Google Earth imagery was used.

Figure 2Canopy height models (CHM) for two 50 m×50 m field sites in 1 m grid resolution in the eastern Swiss Alps with (a) high canopy coverage and (b) low canopy coverage (for detailed site descriptions see Moeser et al.2014).

During the winters of 2012/2013 and 2013/2014, snow depth was measured immediately after every storm with a greater than 15 cm depth of snowfall in the open site. In total, nine storm events met the following prestorm and storm conditions that allowed for indirect interception measurements: (1) no snow in the canopy prior to a storm event, (2) defined crust on the underlying snow, and (3) minimal wind redistribution during the storm cycle. New snow was measured down to the prior snow layer crust from the top of the newly fallen snow layer to represent the total snow interception. Total snowfall was measured at the open field sites. Snow interception was obtained by subtracting the total snowfall measured in the forest from the total snowfall measured at the open field site. The extensive measurement data set used in this study is described in high detail in . Preprocessing resulted in 13 994 usable individual measurements from which 60 site-based mean and standard deviation values of snow interception were computed. These 60 values were then utilized to develop the interception parameterizations. For all individual measurements, a mean snow interception efficiency (interception and/or new snowfall open) of 42 % was measured with values ranging from 0 % to 100 %. The probability distribution function (pdf) of all snow interception data can be fitted with a normal distribution (positive part) with a root mean square error (RMSE) of the quantiles between both distributions of 0.6 cm and a Pearson correlation r of 0.99 for the quantiles (Fig. 3). The average storm values of air temperatures covered cold (−12.1C) to mild (−1.9C) conditions.

Figure 3Probability density functions (pdf's) of all individual snow depth interception measurements used for the development (Swiss (CH) data set) and for the validation of the parameterizations (French (F) and US (US) data sets). The dashed lines indicate a theoretical normal pdf for the corresponding data set.

A 1 m resolution gridded lidar DSM was generated from a flyover in the summer of 2010 and encompasses all eastern Swiss Alps field sites (see Fig. 1a for the extent). The initial point cloud had an average density of 36 points m−2 (all returns) and a shot density of 19 points m−2 (last returns only). The 1 m resolution lidar DSM is used for the derivation of the canopy structure metrics, the standard deviation of the DSM (σz), and the spatial mean sky view factor (Fsky) over each 50 m×50 m field site.

## 2.2 Rocky Mountains of northern Utah, US

For the first validation data set, indirect interception measurements were collected at Utah State University's T.W. Daniel Experimental Forest (TWDEF; 41.86 N, 111.50 W), which is located at ∼2700 m a.s.l. in the Rocky Mountains of northern Utah (Fig. 1b). The forest stand is predominantly coniferous and is composed of Engelmann spruce (Picea engelmannii) and subalpine fir (Abies lasiocarpa). However, deciduous quaking aspen (Populus tremuloides) forest stands are also present. Mean annual air temperature is approximately 4 C and mean annual precipitation is approximately 1080 mm . On average 80 % of the precipitation falls as snow. Similar to the sites in the eastern Swiss Alps, two forested sites and one nonforested site were chosen to limit the influences of slope and topographic shading while capturing diversity in canopy density and canopy structure. At both forested sites, measurements were taken along 20 m forested transects every 0.5 m before and after storm events. The after-storm event transect was parallel to the before-storm event transect but displaced by 0.5 m to avoid impacts from the before-storm event transect (yellow inlet in Fig. 1b corresponds to each yellow dot). At one nonforested reference site (open field site) (see blue dots in Fig. 1b) several disordered measurements were conducted within a fenced meadow site (20 m×20 m) (see blue dot in Fig. 1b). Additionally, an automatic weather station nearby provided continuous measurements (Usu Doc Daniel SNOTEL site) (purple dot in Fig. 1b). Because the purpose of the Utah measurement campaigns was not to measure snow interception but rather to investigate the spatial variability of snow characteristics below different forest canopies , the derivation of snow interception differed slightly from the Swiss sites. Accumulated snowfall was first estimated as the difference between prestorm and poststorm total snow depth. Then snow interception was calculated by subtracting the total snowfall derived in the forest from the total snowfall derived at the open field site.

During winter 2015/2016 several measurement campaigns took place. We selected those campaigns that allowed us to reliably derive snow interception from total snow depth measurements before and after storm events. At one of the forested sites we used four parallel 20 m transects (i.e., two storm events) and at a second forested site two parallel 20 m transects (i.e., one storm event). The total snow depth was also measured every time at the nonforested meadow location (open site). Poststorm measurements were made between approximately 1 and 3 d after a recent snowfall, but the total time period between every first and second campaign lasted several days including multiple snowfalls. The storm events were also temporally close, so that the trees may not have been snow free prior to the new snowfall. As such, unloading and snow settling may have influenced these measurements. After parsing the data to further reduce such influences, 95 individual interception measurements remained, resulting in three site-based mean and standard deviation values. For all individual measurements, a mean snow interception efficiency of 33 % was measured, with values ranging from 2 % to 93 %. The pdf of all individual snow interception data can be similarly well fitted with a normal distribution (positive part) with a RMSE of the quantiles between both distributions of 1.3 cm and a Pearson correlation r of 0.98 for the quantiles (Fig. 3). Average storm values of air temperatures covered cold (−7.3C) to mild (−1.4C) conditions.

A 1 m resolution gridded lidar DSM was generated from a flyover in July 2009 and encompasses all field sites (see Fig. 1b for the extent). On average the initial point cloud had 7 returns m−2 and 5 last returns m−2 (shot density). The 1 m resolution lidar DSM is used for the derivation of the canopy structure metrics σz and Fsky over each 20 m transect (field site).

## 2.3 Southeastern French Alps

For the second validation data set, indirect interception measurements were collected in a coniferous forest stand next to the midaltitude experimental site Col de Porte (45.30 N, 5.77 E) at 1325 m a.s.l. in the Chartreuse mountain range in the French Alps (more site details in ). The forest stand is dominated by Norway spruce (Picea abies) with young silver fir (Abies alba) in the understory. Small deciduous trees are present along the northwestern border of the experimental site. Mean annual air temperature is ∼6C and the average solid precipitation at Col de Porte is 644 mm per year. All snow depth measurements were taken by the Snow Research Center (Centre d'Etudes de la Neige (CEN)) in Grenoble, France, as part of the Labex SNOUF project (SNow Under Forest) (Fig. 1c). There were three 8 m transects, each consisting of eight 1 m×0.39 m wooden boxes that were aligned along the north, south, and west axes of the field site. New snow depth was measured inside each box after a storm event, and the box was then cleared of snow. Open-site new snow depth measurements were obtained from snow board measurements at the experimental site. The boards were cleaned after each precipitation event. Interception was then derived as the difference between the open-site and undercanopy new snow box measurements.

During winter 2017/2018 several measurement campaigns were conducted. Four snowstorm events were selected after which new snow depth was measured in all boxes. Snow depth was collected after a major storm event took place. Unloading was visually observed from webcams and had a minimal influence on the measurements. A total of 96 individual interception measurements (4×24 measurements) resulted in four site-based mean and standard deviation values. For the individual measurements, a mean snow interception efficiency of 66 % was measured with values ranging from 1 % to 94 %. The pdf of all snow interception data can be roughly fitted with a normal distribution (positive part) with a RMSE of the quantiles between both distributions of 1.1 cm and a Pearson correlation r of 0.96 for the quantiles (Fig. 3). The average storm values of air temperatures covered mild (−0.9C) to warm (1.7 C) conditions.

A 1 m resolution gridded lidar DSM was generated from flyovers between 30 August and 2 September 2016, which encompassed the entire Col de Porte experimental site (IRSTEA, Grenoble (see Fig. 1c)). The initial lidar point cloud had an average density of 24 points m−2 and a shot density of 17 points m−2 (last return). The initial point cloud right at the transects had an average density of 42 points m−2 and a shot density of 25 points m−2 (last return). The 1 m resolution lidar DSM is used for the derivation of the canopy structure metrics σz and Fsky over the three 8 m transects.

3 Methods

Subgrid parameterizations were derived for site means and standard deviations of the snow interception using forest structure metrics and open-site snowfall. We parameterize mean and spatial variability of snow interception for a model grid cell by accounting for the unresolved underlying forest structure (subgrid parameterization). Forest structure metrics are derived from DSMs to integrate both the terrain elevation and vegetation height.

## 3.1 Forest structure metrics

The sky view factor Fsky describes the proportion of a radiative flux received by an inclined surface patch from the visible part of the sky to that obtained from an unobstructed hemisphere . Fsky is a commonly applied model parameter when computing surface radiation balances and can be easily computed for large areas from DSMs. Fsky integrates previously applied forest structure metrics, such as total open area and mean distance to canopy, because this parameter is able to account for distance, size, and orientation of individual surface (or canopy) patches . We therefore selected Fsky to parameterize the site mean and standard deviation of snow interception (IHS, σHS). Here, we compute Fsky from view factors which are geometrically derived quantities. They can be computed by the numerical methods described within the radiosity approach for the shortwave (SW) radiation balance over complex topography and were originally introduced to describe the radiant energy exchange between surfaces in thermal engineering . Thereby, solve the double area integral using uniform but adaptive area subdivision for surface patches AI, AJ. Fsky for each surface patch AI is one minus the sum over all N view factors FIJ by assuming the sky is one large surface patch. Fsky is computed for each fine-scale grid cell of the DSM as follows:

$\begin{array}{}\text{(1)}& \begin{array}{rl}{F}_{I,\mathrm{sky}}& =\phantom{\rule{0.125em}{0ex}}\mathrm{1}-\sum _{J=\mathrm{1}}^{N}{F}_{IJ}\\ & =\phantom{\rule{0.125em}{0ex}}\mathrm{1}-\sum _{J=\mathrm{1}}^{N}\phantom{\rule{0.125em}{0ex}}\frac{\mathrm{1}}{{A}_{I}}\phantom{\rule{0.33em}{0ex}}\underset{{A}_{I}}{\int }\phantom{\rule{0.125em}{0ex}}\underset{{A}_{J}}{\int }\phantom{\rule{0.125em}{0ex}}\frac{\mathrm{cos}{\mathit{\vartheta }}_{I}\phantom{\rule{0.33em}{0ex}}\mathrm{cos}{\mathit{\vartheta }}_{J}}{\mathit{\pi }\phantom{\rule{0.33em}{0ex}}{r}_{IJ}^{\mathrm{2}}}\phantom{\rule{0.125em}{0ex}}\mathrm{d}{A}_{I}\phantom{\rule{0.125em}{0ex}}\mathrm{d}{A}_{J}\phantom{\rule{0.33em}{0ex}}.\end{array}\end{array}$

Deriving Fsky via Eq. (1) can account for holes in the surface, i.e., small gaps between leaves and branches in forest canopy, provided the DSM is of a high enough resolution to capture this. In this study, the employed DSMs did not resolve small gaps between branches. Common methods to derive Fsky for forested regions are from sine and cosine weighted proportions of sky pixels of HP or SP as suggested, e.g., by or from LAI (e.g., Roesch et al.2001). However, compared to computing Fsky on DSMs these methods rely on extensive field work.

The main advantage of deriving Fsky on DSMs is that Fsky can be derived spatially by averaging all fine-scale Fsky within a coarse grid cell. Here, we use the spatial mean of the sky view factor Fsky Eq. (1) over a field site which is comparable to the spatial mean canopy openness.

The second forest structure metric selected was the standard deviation of the DSM σz of a field site. Though not totally uncorrelated from the spatial mean Fsky (Pearson $r=-\mathrm{0.48}$), we selected σz to serve in coarse-scale models that are not able to rely on computational expensive precomputations of Fsky on fine scales, such as land surface models covering regions of several hundreds to thousands of kilometers. σz is thought to represent the spatial variability of canopy height and terrain elevation of the field site (or model domain).

## 3.2 Subgrid parameterization for forest canopy interception

Modeling the impact of forest canopy on snow accumulation on the ground involves several processes such as interception, unloading, melt and drip, and sublimation. Here, we present novel models for the spatial mean and standard deviation of snow interception. Modeling not only the mean but the standard deviation of snow interception within a grid cell or model domain opens new possibilities for describing the spatially varying snow cover in large grid cells. Empirical parameterizations for site mean and standard deviation of snow interception are derived from the 60 measured mean and standard deviation values from the Swiss data set. Estimates derived using the new models were validated from a comparison to the mean and standard deviation values from the French and US field sites.

Snow interception I was modeled as snow depth HS, i.e., IHS, and not as snow water equivalent SWE, i.e., ISWE. Snow interception models for SWE would be advantageous for model applications because this removes the uncertainties of the consequent empirical snow density parameterization in each model application. However, similar spatial SWE interception measurements comparable to the extensive spatial snow depth interception data set from Switzerland are not available at the moment. The reason similar SWE data sets do not exist is probably that SWE measurements require much more effort and are more time-consuming. We further refrained from deriving a spatial SWE data set from the spatial HS interception data set to avoid any potential errors being introduced when empirically converting measured HS values to SWE. Thus, any future snow density model developments should not affect our snow interception models. Previous interception models estimated new snow density to convert HS into SWE. Models of new snow density typically rely on average storm temperature. Thus, converting HS empirically to SWE and then developing an empirical interception model introduces additional uncertainty. Prior work has shown a standard error of 9.31 kg m−3 when using estimates of density . As such, the snow interception parameterizations developed here are for HS.

From here on, all references will be to site values (mean and standard deviation) without explicitly mentioning the “mean”, unless otherwise stated.

## 3.3 Performance measures

We use a variety of measures to validate the parameterizations, namely the RMSE, normalized root mean square error (NRMSE, normalized by the range of measured data (max–min)), mean absolute error (MAE), the mean absolute percentage error (MAPE, absolute bias with measured minus parameterized normalized with measurements), the mean percentage error (MPE, bias with measured-parameterized normalized with measurements), and the Pearson correlation coefficient r as a measure for correlation. Finally, we evaluate the performance of our parameterizations by analyzing the pdf. We use the two-sample Kolmogorov–Smirnov test (K-S test) statistic values D (Yakir2013) for the pdf (nonparametric method) and compute the NRMSE for quantile–quantile plots (NRMSEquant, normalized by the range of measured quantiles (max–min))) for probabilities with values in the range of [0.1,0.9].

4 Results

## 4.1 Grid cell mean snow interception

### 4.1.1 Model for grid cell mean intercepted snow depth

We parameterized the grid cell mean intercepted snow depth (IHS) by scaling open-site accumulated snowfall PHS using the forest structure metrics Fsky and σz. From these three variables, the interception measurements of the development data set correlated best with PHS (r=0.70). Snow interception efficiency (IHSPHS) correlations were slightly stronger for σz (r=0.71) than for Fsky ($r=-\mathrm{0.69}$).

While it is clear that accumulated snowfall is the key parameter for modeling snow interception by forest canopy and as such regulates its magnitude, the shape of the interception curve is predominantly controlled by forest canopy parameters and the interception model form itself. This behavior of the interception curve has been recently demonstrated by comparing various SWE interception models at single forest sites . To decide on the interception model form we considered previously commonly applied functional relationships with accumulated snowfall such as those from and , as well as simple relationships such as a power law. Together with our observed correlations of the forest structure metrics Fsky and σz with snow interception efficiency, we developed two statistical parameterizations for IHS using two different base functions to scale PHS with either Fsky and σz (Eq. 2) or with only σz (Eq. 3) as follows:

$\begin{array}{}\text{(2)}& {I}_{\mathrm{HS}}\phantom{\rule{0.125em}{0ex}}=\phantom{\rule{0.125em}{0ex}}{P}_{\mathrm{HS}}^{a}\phantom{\rule{0.125em}{0ex}}b\phantom{\rule{0.125em}{0ex}}\frac{\left(\mathrm{1}-{F}_{\mathrm{sky}}{\right)}^{c}\phantom{\rule{0.125em}{0ex}}{\mathit{\sigma }}_{z}^{c}}{\mathrm{1}\phantom{\rule{0.125em}{0ex}}+\phantom{\rule{0.125em}{0ex}}\mathrm{exp}\left(-d\phantom{\rule{0.125em}{0ex}}\left({P}_{\mathrm{HS}}\phantom{\rule{0.125em}{0ex}}-\phantom{\rule{0.125em}{0ex}}f\right)\right)},\end{array}$

with constant parameters: a=0.09 (±1.08), b=0.19 (±0.79), c=0.72 (±0.11), d=0.13 (±0.04) and f=16.44 (±16.33) and

$\begin{array}{}\text{(3)}& {I}_{\mathrm{HS}}\phantom{\rule{0.125em}{0ex}}=\phantom{\rule{0.125em}{0ex}}{P}_{\mathrm{HS}}^{{a}^{\star }}\phantom{\rule{0.125em}{0ex}}{b}^{\star }\phantom{\rule{0.125em}{0ex}}{\mathit{\sigma }}_{z}^{{c}^{\star }},\end{array}$

with constant parameters: ${a}^{\star }=\mathrm{0.82}$ (±0.12), ${b}^{\star }=\mathrm{0.0035}$ (±0.0036), and ${c}^{\star }=\mathrm{0.80}$ (±0.14). The constant parameters result from fitting nonlinear regression models by robust M-estimators using iterated reweighed least squares (see R v3.2.3 statistical programming language and its robustbase v0.92-5 package; ). The 90 % confidence intervals of the parameters are given in parentheses. In both equations PHS and σz are in centimeters.

The accuracy of a derived model between IHS and PHS depended upon the forest structure metrics and the underlying function applied in the potential models. While we investigated previously suggested functional dependencies for the amount of storm snowfall, the best performances were seen when the base function between IHS and PHS was either a power law or a combination of a power law with an exponential dependence. Similar base functions were obtained for fine-scale ISWE models by (exponential) and recently by (power law).

Estimated IHS values from Eq. (2) or (3) increase with increasing PHS, increasing σz or decreasing Fsky. This implies that with increasing forest density (i.e., larger σz), IHS increases faster with increasing PHS. Note that a lower Fsky value denotes more pronounced forest gaps since it is derived from aerial lidar DSM here.

Equations (2) and (3) differ in two ways. First, Eq. (2) incorporates the functional dependency for increasing PHS that snow interception efficiency (interception and/or snowfall) increases with increasing precipitation due to snow bridging between branches until a maximum is reached after which it decreases due to the bending of branches under the load (sigmoid curve as suggested by ). Additionally, a power-law dependency for accumulated open-site storm snowfall is applied to force the sigmoid distribution to zero at very small snowfall events. The sigmoid curve alone is not able to reach zero, potentially breaking the mass balance. In contrast, Eq. (3) solely employs the power-law dependency between IHS and accumulated open-site storm snowfall PHS. The second difference between both equations is that Eq. (2) uses both forest structure metrics (Fsky and σz), whereas Eq. (3) only uses σz. Equation (2) is thus more “complex” and necessitates more time to derive both forest structure parameters, whereas Eq. (3) has a more “compact” form and solely necessitates the estimation of σz.

To evaluate model performances with respect to a simple baseline interception estimate, we linearly fitted the key parameter accumulated snowfall to measured snow interception by assuming constant impact of forest canopy parameters; i.e., IHS= ccPHS with constant fit parameter cc=0.40 (±0.03).

### 4.1.2 Validation of model for grid cell mean intercepted snow depth

Performances of both the newly developed snow interception IHS models (Eqs. 2 and 3) were compared to the IHS measurements from the development data set (Switzerland), as well as the IHS measurements from the combined validation data sets (France and US). In Figs. 4 to 6 we differentiate the validation data set from the development data set by using a black outline around the symbols (validation) instead of colored circles (development). Squares represent the data set from the US and diamonds represent the data set from France.

Figure 4Measured and parameterized site means of intercepted snow depth, i.e., spatially averaged over each site and for each storm date. Parameterized using (a) Eq. (2) and (b) Eq. (3) as a function of the site means of standard deviation of the lidar DSM σz (in color) as well as open-site snowstorm precipitation (size of symbols). Circles represent the development data set from Switzerland, symbols with a black border represent the validation data sets with squares for that from the US, and diamonds for data sets from France.

Figure 4 shows that there is good agreement between IHS and measured interception at all sites for both models. Overall error statistics show good performances for the development and the validation data sets with low absolute errors (e.g., all MAE  1.2 cm), strong correlations (all r≥0.89), and low distribution errors (e.g., all NRMSEquant<8 %) (Table 1). In contrast to the validation data sets the performance statistics for the development data set are slightly reduced for the more compact model (Eq. 3) compared to the more complex model (Eq. 2). Overall, the performance metrics suggest that the simple baseline interception model is worse for both the development and the validation data sets (I(c) in Table 1).

Table 1Performance measures between measurement and parameterization of (I) spatial-mean snow depth interception IHS with (a) Eq. (2), (b) Eq. (3), and (c) a baseline model and of (II) standard deviation of snow depth interception ${\mathit{\sigma }}_{{I}_{\mathrm{HS}}}$ with (a) Eq. (4) and (b) a baseline model. Statistics are shown for the development data set from the eastern Swiss Alps (CH) and for the combined validation data set (US&F).

Figure 5 reveals overall similar performances for both parameterizations as a function of accumulated new snowfall. However, small differences between both parameterizations are visible in the extremes, i.e., for very low and very large IHS and PHS. The bias for the largest PHS (US data set) is larger for the more compact parameterization (Eq. 3), whereas for the smallest PHS (data set from France) the bias is slightly larger for the more complex parameterization (Eq. 2). The bias is more pronounced with regard to the corresponding interception efficiencies, shown in Fig. 5d–f, where the largest bias for the smallest PHS for the complex parameterization (Eq. 2) is −0.24 compared to 0.21 for the more compact parameterization (Eq. 3).

Figure 5Snow depth interception IHS (a, b, c) and interception efficiency IHSPHS (d, e, f) as a function of accumulated open-site snowstorm precipitation PHS and standard deviation of the lidar DSM σz (in color). The y axis of the first column shows measured data, the second column shows model output with Eq. (2), and the third model output with Eq. (3). Site means for each storm event are depicted with colored circles for the development data set from Switzerland, and symbols with a black border depict the validation data sets, with squares for data from the US and diamonds for data from France. Storm means (in PHS bins) are shown in black.

## 4.2 Grid cell standard deviation of snow interception

### 4.2.1 Model for standard deviation of snow depth interception

We parameterized the standard deviation of snow depth interception ${\mathit{\sigma }}_{{I}_{\mathrm{HS}}}$ by scaling PHS using the forest structure metric σz. ${\mathit{\sigma }}_{{I}_{\mathrm{HS}}}$ of the development data set correlated best with PHS (r=0.82). The correlation with mean snow interception IHS was less pronounced (r=0.33). ${\mathit{\sigma }}_{{I}_{\mathrm{HS}}}$ normalized with PHS correlated much better with σz ($r=-\mathrm{0.68}$) than with Fsky (r=0.13).

Building upon the observed power-law functional dependency between mean snow interception IHS and PHS and the observed relationships and correlations for ${\mathit{\sigma }}_{{I}_{\mathrm{HS}}}$, we scaled a power-law function for PHS with the standard deviation of the DSM σz in order to parameterize ${\mathit{\sigma }}_{{I}_{\mathrm{HS}}}$ as follows:

$\begin{array}{}\text{(4)}& {\mathit{\sigma }}_{{I}_{\mathrm{HS}}}\phantom{\rule{0.125em}{0ex}}=\phantom{\rule{0.125em}{0ex}}{P}_{\mathrm{HS}}^{g}\phantom{\rule{0.125em}{0ex}}\frac{h}{\mathrm{1}+{\mathit{\sigma }}_{z}^{j}}\phantom{\rule{0.33em}{0ex}}.\end{array}$

Constant parameters g=0.78 (±0.10), h=13.40 (±11.64), and j=0.53 (±0.12) result from fitting a nonlinear regression model, similar to the derivation of IHS from Eqs. (2) and (3). The 90 % confidence intervals of the parameters are given in parentheses. In Eq. (4) PHS and σz are in centimeters.

${\mathit{\sigma }}_{{I}_{\mathrm{HS}}}$ derived from Eq. (4) increases with increasing PHS or decreasing σz. This implies that with decreasing σz (decreasing forest density), the spatial variability in snow interception increases faster with increasing PHS. The opposite correlation was found between σz and mean snow interception IHS. For a σz converging to zero, modeled ${\mathit{\sigma }}_{{I}_{\mathrm{HS}}}$ via Eq. (4) approaches a constant fraction of precipitation.

Similar to the derivation of a baseline estimate for our IHS models, we linearly fitted the accumulated snowfall to the measured standard deviation of snow interception to evaluate the model performance with respect to a simple baseline estimate for the standard deviation of snow interception. This resulted in ${\mathit{\sigma }}_{{I}_{\mathrm{HS}}}=jj\phantom{\rule{0.25em}{0ex}}{P}_{\mathrm{HS}}$ with constant fit parameter jj=0.20 (±0.01).

### 4.2.2 Validation of model for standard deviation of snow depth interception

Overall, modeled and measured ${\mathit{\sigma }}_{{I}_{\mathrm{HS}}}$ agree well (Fig. 6). Error statistics show good performances for the development and the validation data set with low absolute errors (e.g., all MAE  0.63 cm), strong correlations (all r≥0.92), and low distribution errors (e.g., NRMSEquant<10 %) (Table 1). However, performances are less accurate for the validation data set than for the development data set (e.g., MAE of 0.63 cm as opposed to 0.45 cm and NRMSEquant of 10 % as opposed to 4 %). This was caused by a potential outlier in the validation data set from the US. During one measurement campaign, an open-site accumulated storm snowfall PHS was not available at the same date as the undercanopy measurements. Therefore, this value was estimated from a local automatic weather station (Usu Doc Daniel SNOTEL site; purple dot in Fig. 1b). Additional measurement uncertainty (at the Utah site) was also introduced since interception estimates were integrated values over several snowstorms that occurred during the 13 d between presnowfall and postsnowfall measurement campaigns. When this outlier is removed from the validation data set, performance statistics improve considerably and converge towards the errors of the development data set (cf. MAE decreases to 0.35 cm and the NRMSEquant to 5 %).

Figure 6Measured and parameterized standard deviation of snow depth interception ${\mathit{\sigma }}_{{I}_{\mathrm{HS}}}$ at each site and for each storm date. Parameterized using Eq. (4) as a function of site means of standard deviation of the lidar DSM σz (in color) as well as open-site snowstorm precipitation (size of symbols). Circles represent the development data set from Switzerland, symbols with a black border represent the validation data sets with squares for data from the US and diamonds for data from France.

Overall, the performance of the baseline model for ${\mathit{\sigma }}_{{I}_{\mathrm{HS}}}$ is worse than that of our model performance (II(b) in Table 1). However, because one observed ${\mathit{\sigma }}_{{I}_{\mathrm{HS}}}$ of the validation data set from the US (2.9 cm) was better estimated by the baseline model than our model (4 cm compared to 5.2 cm), the NRMSE and RMSE for the baseline estimates were somewhat better.

To compare modeled (Eqs. 2 and 4) and measured data set mean values from each geographic location (Switzerland, US, and France), we averaged all site values to derive an overall mean of IHS and ${\mathit{\sigma }}_{{I}_{\mathrm{HS}}}$ for each location. The coefficient of variation (CV, description of variability) (CV${}_{{I}_{\mathrm{HS}}}={\mathit{\sigma }}_{{I}_{\mathrm{HS}}}/{I}_{\mathrm{HS}}$) was also calculated for each of the three geographic locations. For the Swiss development data set, the same overall mean, standard deviation, and CV for measured and modeled snow interception was calculated (mean of 9.4 cm, standard deviation of 4.5 cm, and CV of 0.51). For the validation data sets we obtained slightly larger values for modeled IHS (9.3 cm), modeled ${\mathit{\sigma }}_{{I}_{\mathrm{HS}}}$ (3.7 cm), and modeled CV${}_{{I}_{\mathrm{HS}}}$ (0.38) than measured IHS (9.2 cm), measured ${\mathit{\sigma }}_{{I}_{\mathrm{HS}}}$ (3.2 cm), and measured CV${}_{{I}_{\mathrm{HS}}}$ (0.35). If the potential outlying data point from Utah is removed, the same overall modeled and measured mean CV${}_{{I}_{\mathrm{HS}}}$ (0.32) is found along with very close values of the modeled and measured mean IHS (9.8 cm versus 9.9 cm) and modeled and measured ${\mathit{\sigma }}_{{I}_{\mathrm{HS}}}$ values (3.4 cm versus 3.3 cm).

5 Discussion

We proposed two empirical models for spatial mean interception IHS to be employed in hydrological, climate, and weather applications. One model is a more compact model, Eq. (3). This model uses a power-law dependency between IHS and accumulated storm precipitation PHS that is scaled by one forest structure metric: the standard deviation of the DSM σz. The other model, Eq. (2), integrates a more complex parameterization by using a combination of a power law with an exponential dependence similar to the one suggested by for PHS and is scaled by two forest metrics, namely the sky view factor Fsky in combination with σz. For both IHS models, interception increases faster with increasing snowfall when forest density increases (i.e., larger σz). In the more complex model, increasing forest density is implemented by increasing σz and decreasing Fsky. Though Fsky can be precomputed and is temporally valid for many years (unless the forest structure changes due to logging, fires, insect infestations, or other forest disturbances), computing Fsky over large scales and/or with fine resolutions is more computationally demanding than for σz . A subgrid parameterization for the sky view factor of coarse-scale DSMs over the forest canopy would eliminate the precomputation of sky view factors on fine-scale DSMs. Such a subgrid parameterization for sky view factors over forest canopy could be similarly set up as previously done for alpine topography and would lead us towards a global map of sky view factors (cf. Helbig and Löwe2014).

In general, more differences between the compact and more complex modeling approaches only displayed at the extremes. For instance, for small storm precipitation values (PHS=3 cm), the more compact parameterization performs slightly better, whereas for very large storms (PHS=43 cm) the more complex model displayed improved performance. The choice of one of these two models thus depends on the focus range of precipitation values and available computational resources.

Our choice for the functional form of PHS differs from previous parameterizations for snow interception solely using the sigmoid growth $\sim \mathrm{1}/\left(\mathrm{1}+\mathrm{exp}\left(-k\left(P-{P}_{\mathrm{0}}\right)\right)\right)$ or an exponential form $\sim \left(\mathrm{1}-\mathrm{exp}\left(-k\left(P-{P}_{\mathrm{0}}\right)\right)\right)$ (e.g., Aston1993; Hedstrom and Pomeroy1998) with increasing precipitation. While the base function of worked better for , a drawback of this relationship is that interception does not become exactly zero for a zero snowfall amount. To account for this, the model becomes complicated when applied to discrete model time steps . For this reason, selected the relationship proposed by for their parameterization of snow interception. However, the functional form of the model does not account for snow bridging or branch bending, thus modeling interception efficiency as decreasing through time. We also compared means and standard deviations over all the sites as a function of forest metrics and found that the use of storm means can introduce precipitation dependencies that might originate from an insufficient number of sites showing similar forest canopy structure parameter values for a given precipitation (cf. black line compared to colored dots in Fig. 5). Based on the functional dependencies revealed by analyzing our data as a function of PHS and forest structure metrics, a simple power law was able to describe the spatial mean PHS dependency of snow interception (cf. Eq. 3). The equation displayed that with increasing PHS, IHS increases. This is less pronounced with smaller σz or larger Fsky values (Fig. 5). Very recently, a storm event power-law dependency was also found to best describe fine-scale SWE interception in a maritime snow climate . Our base functions for site means and standard deviations thus bear some similarity to previously developed fine-scale snow interception models. Despite an ongoing debate regarding the proper representation of interception, we believe that the interception models presented here have the advantage that they could be applied in various model applications for horizontal grid cell resolutions larger than a few tens of meters. Due to the lack of measurements over larger scales a validation remains impossible at the moment.

We derived just one empirical model for the standard deviation of snow interception ${\mathit{\sigma }}_{{I}_{\mathrm{HS}}}$ that uses a power-law dependency on accumulated storm precipitation PHS scaled by one forest structure metric, namely the standard deviation of the DSM σz. We also tested a more complex model for ${\mathit{\sigma }}_{{I}_{\mathrm{HS}}}$ using both forest metrics (Fsky and σz) that integrates a power-law dependency of PHS. However, model performances for the validation data set did not differ considerably from the ones for the more compact model. Therefore, we propose the more compact parameterization for ${\mathit{\sigma }}_{{I}_{\mathrm{HS}}}$ (Eq. 4) to facilitate broad model applications.

By using Fsky and σz derived from DSMs as forest structure metrics, we focused on the overall shape of the forest. This simplification is similar to the assumption by for solar transmissivity in forests under cloudless sky conditions. They assumed the fraction of solar radiation blocked by the canopy was equal to 1−Vf with Vf being defined as the fraction of the sky visible from beneath the canopy. Our simplification is also in line with previous suggestions. Primarily, to reliably describe interception by forest canopy over larger areas, the larger-scale canopy structure needs to be taken into account instead of only using point-based canopy structure parameters (e.g., Varhola et al.2010; Moeser et al.2016). We proposed to calculate Fsky and σz on DSMs rather than on CHMs to account for terrain and vegetation height. This results from our correlation analysis for measurement data collected in rather flat field forest sites (Sect. 2) and should be verified once spatial snow interception measurements become available in steeper terrain and over larger length scales.

The models for IHS and ${\mathit{\sigma }}_{{I}_{\mathrm{HS}}}$ were statistically derived from measured snow interception data gathered in the eastern Swiss Alps. Naturally, empirically derived parameterizations can only describe data variability covered by the data set. However, even though the parameterizations were developed empirically, we could display that the parameterizations perform well for two disparate, independent snow interception data sets collected in geographically different regions, different snow climates, coniferous tree species, and prevailing weather conditions during collection of the validation data sets (French Alps and Rocky Mountains, US). For instance, in the French Alps, rather warm to mild winter weather conditions predominated, whereas rather mild to cold weather prevailed during the campaigns in the Rocky Mountains of northern Utah in the US. Though snow cohesion and adhesion are clearly temperature dependent, we did not observe decreases in overall performances under these differing weather conditions for our two IHS models, which do not include air temperature. In contrast, in a maritime (warm) snow climate correlations between air temperature and snow interception were recently found by . In addition to the spread in observed temperature conditions, our ranges of accumulated snowstorm PHS values of the development data set are fairly broad (e.g., PHS between 10 and 40 cm). The measurements of the validation data set are well within the range of the development data set values but also cover extremes, such as one very small (PHS=3 cm) and one very large snowfall (PHS=43 cm) (cf. Fig. 3). It is thus reassuring that our models perform sufficiently well in varying climate regions; however, more validation data sets would be advantageous especially in regions experiencing extreme climates such as the cold arctic or warm maritime forests. Despite the existing variability in the data set, more spatial snow interception measurements would clearly help to increase the robustness of our empirical parameterizations.

To date, interception models have been created for SWE instead of snow depth and were mostly point models instead of spatial mean interception parameterizations. As such, a comparative assessment (beyond the independent validation sets in the body of this paper) of our models to absolute performance measures of previous interception models was difficult. However, we calculated relative error estimates for an intermodel comparison with two SWE interception models. We selected the empirical, recently developed SWE model from as well as the 50 m×50 m stratified SWE model (for 50 m×50 m grid cell size) from . The model utilized the same Swiss data set as this study, which is currently the largest set of spatial interception measurements in the world. The model error estimates were calculated for a subset of their data set which included three snowfall events and interception values acquired at three elevations under mild temperature conditions in the McKenzie River basin, Oregon, US (for details on the data see Table 2 in Roth and Nolin2019). We estimated an NRMSE of 28.9 %, MPE of −5.7 %, and MAPE of 31.2 % for the three modeled and measured interception values. The model error estimates were calculated for a subset of the Swiss data set consisting of 34 spatial mean observed interception values (50 m×50 m) and 34 parameterized values. We estimated an NRMSE of 9.3 %, MPE of −16.5 %, and MAPE of 23.5 %. Compared to previous SWE interception models, we obtained improved performances (using means of error estimates over I(a) and I(b), respectively in Table 1). The fairest comparison for our mean error estimates is with the stratified SWE model of . Thus, our mean error estimates reduce as follows for the more complex model (Eq. 2) and the more compact model (Eq. 3): by 9 % and 4 % in the NRMSE, by 60 % and 75 % in the MPE, and by 40 % and 50 % in the MAPE. Our improved model performances, when compared to prior interception models in tandem with a good performance for two distinctly different validation data sets, lend validity to improving coarse-scale climate and hydrologic (watershed and snow) model applications.

Despite the overall good performance of the models, we observed differences between the two validation data sets. The data set collected in France shows improved error statistics for snow interception IHS (e.g., for Eq. (3): RMSE=0.35 cm, NRMSE=4 %, MAE=0.26 cm) when compared to the data set collected in the US (e.g., for Eq. (3): RMSE=1.52 cm, NRMSE=14 %, MAE=1.4 cm). In France, intercepted snowstorm depth was measured as the difference of new snow depth in wooden boxes below trees and open-site new snowstorm depth. This was done in relatively short time intervals after a snowstorm. In the US, intercepted snow was inferred from total snow depth before and after a snowstorm event within forests and in an open site. Derived snow interception was often integrated over several storm events due to longer periods between measurement campaigns. Thus, these measurements were potentially influenced by other snow and forest processes such as snow settling, wind redistribution, sublimation, unloading, and melt and drip. Our interception models, however, only calculate how much snow is intercepted at any point in time, which provides the input for other forest snow process models such as for unloading, sublimation, and melt and drip. We thus assume that these processes will be addressed separately, as in all prior interception models . Despite some uncertainties in the validation data set from the US, it allowed for validation in a different snow climate than the French Alps and also covered a large spread in storm snowfall amounts (Fig. 4).

Differences in model performances between the two validation data sets could also be attributed to the more accurate forest structure metrics for the French data set because of a higher resolution lidar DSM (higher point density of 24 returns m−2 returns and 17 last returns m−2) compared to the lidar flyover from the US (on average 7 returns m−2 and 5 last returns m−2). While it is clear that the higher the point cloud density, the greater the potential detail of derived DSMs, 1 m resolution DSMs computed from point clouds above 5 returns m−2 are usually quite consistent and are suitable for deriving coniferous canopy models that allow tree-level analyses . Current available or scheduled countrywide data sets are now around 1–5 returns m−2 , and these densities can be expected to increase thanks to technical improvements in lidar sensors. Since fine-scale DSMs are the only input required to derive the forest structure metrics Fsky and σz, a global applicability of our snow interception models for coniferous forest would be possible with minimal required information.

To understand if the models would also work in other forest types or in disturbed forests, e.g., due to logging, fires, or insect infestations, more snow interception measurements in deciduous, mixed, and disturbed forests are required. Very recently showed that previously published snow interception models developed for coniferous forests from , , and required recalibration to match observed point snow interception observations in a deciduous southern beech Nothofagus stand of the southern Andes. We investigated the performance of our models for two measurement campaigns in a deciduous quaking aspen (Populous tremuloides) forest in our US field site. The measurement setup (20 m transects) was identical to the ones in the coniferous forest at this location (see Sect. 2.2). Though overall the models compared well with the measurements, the model performance was not as good as for the coniferous forest. Because the lidar DSM was acquired in the summer, i.e., with leaves on the trees, the models naturally overestimated IHS and ${\mathit{\sigma }}_{{I}_{\mathrm{HS}}}$. For instance, using the more complex model for IHS (Eq. 2) we obtained a mean bias of −6 cm, whereas when using the more compact model for IHS (Eq. 3) we obtained a mean bias of −8 cm. For ${\mathit{\sigma }}_{{I}_{\mathrm{HS}}}$, the performance was slightly better overall with a mean bias of −3 cm (Eq. 4). While this shows that the performance is clearly lower in such sites, we assume that the performance would be improved when the lidar is acquired in leaf-off conditions.

The lidar-derived DSM sky view factors do not account for small spaces between leaves or branches, which are well accounted for when sky view factors are derived from HP or LAI. In principle, sky view factors that are computed on DSMs represent, depending on the return signal used to create the DSM, a coarser view of the underlying forest canopy. While this increases the overall fine-scale error, we feel that the ability to calculate both our canopy structure metrics in the Cartesian DSM space, which allows an easy model application, far outweighs fine-scale resolution losses.

6 Conclusion and outlook

The statistical models for spatial mean and standard deviation of snow interception presented here are a first step toward a more robust consideration of snow interception for various coarse-scale model applications. They were built upon a very large data set, validated by two other data sets from different geographic regions and snow climates, and performed well for all three sites and under differing weather conditions. For spatial mean interception all NRMSEs were ≤10 % and for the standard deviation of interception all NRMSEs were ≤13 %. Compared to a previous model for spatial mean SWE at 50 m×50 m the presented models for spatial mean snow interception show improved model performances.

In our observed snow interception data sets, as much as 68 % and on average 43 % of the cumulative snowfall (accumulated snowfall of snowfall event in cm) was retained by coniferous forests (interception efficiency (snow interception/accumulated snowfall) of site means) and as much as 14 % and on average 11 % was retained by deciduous forests. These values compare well to previously observed snow interception in coniferous trees reaching up to 60 % of cumulative snowfall and up to 24 % of total annual snowfall in a deciduous forest in the southern Andes .

The empirical models integrate forest parameters that can be derived from fine-scale DSMs, which can be pregenerated and stored for large regions. One of the presented interception models only relies on the standard deviation of the fine-scale DSM, which is a very efficient way to integrate snow interception in coarse-scale models such as land surface models. This could greatly improve current forest albedo estimates and the subsequent surface energy balance for various model applications such as hydrological, weather, and climate predictions.

However, the presented parameterizations were developed and validated for spatial means and standard deviations over horizontal length scales of a few tens of meters. We can only hypothesize that the parameterizations are also valid at coarser length scales due to the use of nonlocal forest structure parameters. Representative nonlocal forest structure parameters require that a DSM of high enough resolution is available to represent the subgrid variability of forest structure in the coarse-scale model grid cell. However, there was and probably is, to date, no validation data available at large spatial scales. The investigated length scale matches current satellite resolutions (e.g., Sentinel and Landsat), which opens further cross validation and deployment opportunities with satellite-derived parameters such as surface albedos and fractional snow-covered area. With parameterizations for both the mean and standard deviation of snow interception by forest canopy, the distribution of intercepted snow depth in forests can now be derived whenever a sufficiently high-resolution DSM is available.

Data availability
Data availability.

The data that support the findings of this study are available from the corresponding author upon reasonable request (norahelbig@gmail.com).

Author contributions
Author contributions.

NH and DM designed the project idea. NH developed the parameterizations. DM developed the measurement design in Switzerland and performed the snow measurements. YL developed the measurement design in France, and YL and LV performed the snow measurements (as part of the Labex SNOUF project). JES was responsible for the project design in France. JMM was responsible for the lidar measurements in France. MT developed the measurement design in the US and performed the snow measurements. NH wrote the paper, and all coauthors contributed to the paper with data, ideas, and writing.

Competing interests
Competing interests.

The authors declare that they have no conflict of interest.

Acknowledgements
Acknowledgements.

We thank Emmanuel Thibert and Éric Mermin, who performed the GPS and tree measurements of the site setup at Col De Porte, and James A. Lutz, Scott B. Jones, Jobie Carlisle, and David. G. Tarboton for their support and the possibility to work at TWDEF. We also thank Matthieu Lafaysse and Tobias Jonas for the valuable discussions. Any opinions, findings, and conclusions or recommendations expressed are those of the author(s) and do not necessarily reflect the views of the Swiss National Science Foundation.

Financial support
Financial support.

Nora Helbig has been partly funded by the Federal Office of the Environment (FOEN). David Moeser has been partially funded by the USGS South Central Climate Adaptation Science Center. Michaela Teich has been partly funded by the Swiss National Science Foundation (grant nos. P2EZP2_155606 and P300P2_171236), the Utah Agricultural Experiment Station (UAES, as part of the “Forests and snow microstructure: key to water supply in the 21st century?” research project), and the Utah State University (USU) Department of Wildland Resources. The Labex SNOUF project has been supported by a grant from the Université Grenoble Alpes as part of the Labex–OSUG@2020 research (grant no. ANR10 LABX 56). The TWDEF lidar data processing has been supported by the National Science Foundation Established Program to Stimulate Competitive Research (NSF EPSCoR) cooperative agreement (grant no. IIA 1208732), which was awarded to USU as part of the State of Utah Research Infrastructure Improvement Award.

Review statement
Review statement.

This paper was edited by Ryan Teuling and reviewed by two anonymous referees.

References

Andreadis, K. M., Storck, P., and Lettenmaier, D. P.: Modeling snow accumulation and ablation processes in forested environments, Water Resour. Res., 45, W05429, https://doi.org/10.1029/2008WR007042, 2009. a

Aston, A. R.: Rainfall interception by eight small trees, J. Hydrol., 42, 383–396, 1993. a

Bartlett, P. A. and Verseghy, D. L.: Modified treatment of intercepted snow improves the simulated forest albedo in the Canadian Land Surface Scheme, Hydrol. Process., 29, 3208–3226, https://doi.org/10.1002/hyp.10431, 2015. a, b

Bründl, M., Bartelt, P., Schneebeli, M., and Flühler, H.: Measuring branch deflection of spruce branches caused by intercepted snow load, Hydrol. Process., 13, 2357–2369, 1999. a

Dirmeyer, P. A., Gao, X., Zhao, M., Guo, Z., Oki, T., and Hanasaki, N.: GSWP-2: Multimodel analysis and implications for our perception of the land surface, B. Am. Meteorol. Soc., 87, 1381–1398, 2006. a

Essery, R.: Large-scale simulations of snow albedo masking by forests, Geophys. Res. Lett., 40, 5521–5525, https://doi.org/10.1002/grl.51008, 2013. a

Essery, R., Rutter, N., Pomeroy, J., Baxter, R., Stähli, M., Gustafsson, D., Barr, A., Bartlett, P., and Elder, K.: SNOWMIP2: An evaluation of forest snow process simulations, B. Am. Meteorol. Soc., 90, 1120–1136, 2009. a, b

Essery, R. L., Pomeroy, J. W., Parviainen, J., and Storck, P.: Sublimation of snow from coniferous forests in a climate model, J. Climate, 16, 1855–1864, 2003. a

Essery, R. L., Bunting, P., Hardy, J., Link, T., Marks, D., Melloh, R., Pomeroy, J., Rowlands, A., and Rutter, N.: Scaling and parametrization of clear-sky solar radiation over complex topography, J. Hydrometeorol., 9, 228–241, 2008. a

Eysn, L., Hollaus, M., Lindberg, E., Berger, F., Monnet, J.-M., Dalponte, M., Kobal, M., Pellegrini, M., Lingua, E., Mongus, D., and Pfeifer, N.: A Benchmark of Lidar-Based Single Tree Detection Methods Using Heterogeneous Forest Data from the Alpine Space, Forests, 6, 1721–1747, 2015. a

Federal Office of Topography Swisstopo: https://www.swisstopo.admin.ch/en/knowledge-facts/geoinformation/lidar-data.html, last access: 22 November 2019. a

Hedstrom, N. R. and Pomeroy, J. W.: Measurements and modelling of snow interception in the boreal forest, Hydrol. Process., 12, 1611–1625, 1998. a, b, c, d, e, f, g, h, i, j, k, l, m

Helbig, N. and Löwe, H.: Parameterization of the spatially averaged sky view factor in complex topography, J. Geophys. Res., 119, 4616–4625, https://doi.org/10.1002/2013JD020892, 2014. a

Helbig, N., Löwe, H., and Lehning, M.: Radiosity approach for the surface radiation balance in complex terrain, J. Atmos. Sci., 66, 2900–2912, https://doi.org/10.1175/2009JAS2940.1, 2009. a, b, c, d, e, f

Hellström, R. A.: Forest cover algorithms for estimating meteorological forcing in a numerical snow model, Hydrol. Process., 14, 3239–3256, 2000. a

Huerta, M. L., Molotch, N. P., and McPhee, J.: Snowfall interception in a deciduous Nothofagus forest and implications for spatial snowpack distribution, Hydrol. Process., 13, 1818–1834, https://doi.org/10.1002/hyp.13439, 2019. a, b, c, d, e, f

Kaartinen, H., Hyyppä, J., Yu, X., Vastaranta, M., Hyyppä, H., Kukko, A., Holopainen, M., Heipke, C., Hirschmugl, M., Morsdorf, F., Næsset, E., Pitkänen, J., Popescu, S., Solberg, S., Wolf, B. M., and Wu, J.-C.: An International Comparison of Individual Tree Detection and Extraction Using Airborne Laser Scanning, Remote Sens., 4, 950–974, https://doi.org/10.3390/rs4040950, 2012. a

Knowles, N., Dettinger, M. D., and Cayan, D. R.: Trends in snowfall versus rainfall in the western United States, J. Climate, 19, 4545–4559, 2006. a

Krinner, G., Derksen, C., Essery, R., Flanner, M., Hagemann, S., Clark, M., Hall, A., Rott, H., Brutel-Vuilmet, C., Kim, H., Ménard, C. B., Mudryk, L., Thackeray, C., Wang, L., Arduini, G., Balsamo, G., Bartlett, P., Boike, J., Boone, A., Chéruy, F., Colin, J., Cuntz, M., Dai, Y., Decharme, B., Derry, J., Ducharne, A., Dutra, E., Fang, X., Fierz, C., Ghattas, J., Gusev, Y., Haverd, V., Kontu, A., Lafaysse, M., Law, R., Lawrence, D., Li, W., Marke, T., Marks, D., Ménégoz, M., Nasonova, O., Nitta, T., Niwano, M., Pomeroy, J., Raleigh, M. S., Schaedler, G., Semenov, V., Smirnova, T. G., Stacke, T., Strasser, U., Svenson, S., Turkov, D., Wang, T., Wever, N., Yuan, H., Zhou, W., and Zhu, D.: ESM-SnowMIP: assessing snow models and quantifying snow-related climate feedbacks, Geosci. Model Dev., 11, 5027–5049, https://doi.org/10.5194/gmd-11-5027-2018, 2018. a

Latvian Geospatial Information Agency: https://www.lgia.gov.lv/lv/Digit%C4%81lais%20virsmas%20modelis, last access: 22 November 2019. a

Lejeune, Y., Dumont, M., Panel, J.-M., Lafaysse, M., Lapalus, P., Le Gac, E., Lesaffre, B., and Morin, S.: 57 years (1960–2017) of snow and meteorological observations from a mid-altitude mountain site (Col de Porte, France, 1325 m of altitude), Earth Syst. Sci. Data, 11, 71–88, https://doi.org/10.5194/essd-11-71-2019, 2019. a

Lundberg, A., Nakai, Y., Thunehed, H., and Halldin, S.: Snow accumulation in forests from ground and remote-sensing data, Hydrol. Process., 18, 1941–1955, https://doi.org/10.1002/hyp.1459, 2004. a, b

Mahat, V. and Tarboton, D. G.: Canopy radiation transmission for an energy balance snowmelt model, Water Resour. Res., 48, 1–16, https://doi.org/10.1029/2011WR010438, 2012. a

Mahat, V. and Tarboton, D. G.: Representation of canopy snow interception, unloading and melt in a parsimonious snowmelt model, Hydrol. Process., 28, 6320–6336, https://doi.org/10.1002/hyp.10116, 2014. a

Moeser, D., Roubinek, J., Schleppi, P., Morsdorf, F., and Jonas, T.: Canopy closure, LAI and radiation transfer from airborne LiDAR synthetic images, Agr. Forest Meteorol., 197, 158–168, 2014. a, b, c, d, e, f

Moeser, D., Morsdorf, F., and Jonas, T.: Novel forest structure metrics from airborne LiDAR data for improved snow interception estimation, Agr. Forest Meteorol., 208, 40–49, 2015a. a, b

Moeser, D., Staehli, M., and Jonas, T.: Improved snow interception modeling using canopy parameters derived from airborne LiDAR data, Water Resour. Res., 51, 5041–5059, https://doi.org/10.1002/2014WR016724, 2015b. a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q, r

Moeser, D., Mazotti, G., Helbig, N., and Jonas, T.: Representing spatial variability of forest snow: Implementation of a new interception model, Water Resour. Res., 52, 5041–5059, https://doi.org/10.1002/2015WR017961, 2016. a, b, c, d

Morin, S., Lejeune, Y., Lesaffre, B., Panel, J.-M., Poncet, D., David, P., and Sudul, M.: An 18-yr long (1993–2011) snow and meteorological dataset from a mid-altitude mountain site (Col de Porte, France, 1325 m alt.) for driving and evaluating snowpack models, Earth Syst. Sci. Data, 4, 13–21, https://doi.org/10.5194/essd-4-13-2012, 2012. a

Morsdorf, F., Kötz, B., Meier, E., Itten, K. I., and Allgöwer, B.: Estimation of LAI and fractional cover from small footprint airborne laser scanning data based on gap fraction, Remote Sens. Environ., 104, 50–61, 2006. a

Pomeroy, J. W. and Schmidt, R. A.: The use of fractal geometry in modelling intercepted snow accumulation and sublimation, in: Proc. Eastern Snow Conference, vol. 50, pp. 1–10, 1993. a, b

Pomeroy, J. W., Parviainen, J., Hedstrom, N., and Gray, D. M.: Coupled modelling of forest snow interception and sublimation, Hydrol. Process., 12, 2317–2337, 1998. a, b

PRISM Climate Group: 30-Year Normals, Oregon State University, available at: http://www.prism.oregonstate.edu/normals/ (created: 7 November 2012; last access: 28 June 2019), 2012. a

Roesch, A. and Roeckner, E.: Assessment of snow cover and surface albedo in the ECHAM5 general circulation model, J. Clim., 19, 3828–3843, 2006. a

Roesch, A., Wild, M., Gilgen, H., and Ohmura, A.: A new snow cover fraction parameterization for the ECHAM4 GCM, Clim. Dyn., 17, 933–946, 2001. a, b, c, d

Roth, T. R. and Nolin, A. W.: Characterizing maritime snow canopy interception in forested mountains, Water Resour. Res., 55, 4564–4581, https://doi.org/10.1029/2018WR024089, 2019. a, b, c, d, e, f, g, h, i, j, k, l

Rousseeuw, P., Croux, C., Todorov, V., Ruckstuhl, A., Salibian-Barrera, M., Verbeke, T., Koller, M., and Maechler, M.: robustbase: Basic Robust Statistics, available at: http://CRAN.R-project.org/package=robustbase (last access: 12 May 2020), r package version 0.92-5, 2015. a

Rutter, N., Essery, R., Pomeroy, J., Altimir, N., Andreadis, K., Baker, I., Barr, A., Bartlett, P., Boone, A., Deng, H., Douville, H., Dutra, E., Elder, K., Ellis, C., Feng, X., Gelfan, A., Goodbody, A., Gusev, Y., Gustafsson, D., Hellström, R., Hirabayashi, Y., Hirota, T., Jonas, T., Koren, V., Kuragina, A., Lettenmaier, D., Li, W.-P., Luce, C., Martin, E., Nasonova, O., Pumpanen, J., Pyles, R. D., Samuelsson, P., Sandells, M., Schädler, G., Shmakin, A., Smirnova, T. G., Stähli, M., Stöckli, R., Strasser, U., Su, H., Suzuki, K., Takata, K., Tanaka, K., Thompson, E., Vesala, T., Viterbo, P., Wiltshire, A., Xia, K., Xue, Y., and Yamazaki, T.: Evaluation of forest snow processes models (SnowMIP2), J. Geophys. Res.-Atmos., 114, D061111, https://doi.org/10.1029/2008JD011063, 2009. a, b, c

Satterlund, D. R. and Haupt, H. F.: Snow catch by conifer crowns, Water Resour. Res., 3, 1035–1039, https://doi.org/10.1029/WR003i004p01035, 1967. a, b, c, d, e

Schmidt, R. A. and Gluns, D. R.: Snowfall interception on branches of three conifer species, Can. J. For. Res., 21, 1262–1269, 1991. a, b, c, d

Sicart, J. E., Essery, R. L. H., Pomeroy, J. W., Hardy, J., Link, T., and Marks, D.: A sensitivity study of daytime net radiation during snowmelt to forest canopy and atmospheric conditions, J. Hydrometeorol., 5, 774–784, 2004. a

Siegel, R. and Howell, J. R.: Thermal Radiation Heat Transfer, Hemisphere Publishing Corp., 1978. a

Storck, P. and Lettenmaier, D. P.: Measurement of snow interception and canopy effects on snow accumulation and m elt in a mountainous maritime climate, Oregon, United States, Water Resour. Res., 38, 11, https://doi.org/10.1029/2002WR001281, 2002. a, b, c

Suzuki, K. and Nakai, Y.: Canopy snow influence on water and energy balances in a coniferous forest plantation in Northern Japan, J. Hydrol., 352, 126–138, 2008. a

Teich, M. and Tarboton, D. G.: TW Daniels Experimental Forest (TWDEF) Lidar, HydroShare, https://doi.org/10.4211/hs.36f3314971a547bc8bc72dc60d6bd03c, 2016.  a

Teich, M., Giunta, A. D., Hagenmuller, P., Bebi, P., Schneebeli, M., and Jenkins, M. J.: Effects of bark beetle attacks on forest snowpack and avalanche formation – Implications for protection forest management, For. Ecol. Manage., 438, 186–203, https://doi.org/10.1016/j.foreco.2019.01.052, 2019. a

Varhola, A., Coops, N., Weiler, M., and Moore, R. D.: Forest canopy effects on snow accumulation and ablation: An integrative review of empirical results, J. Hydrol., 392, 219–233, 2010. a

Vincent, L., Lejeune, Y., M, L., Boone, A., Gac, E. L., Coulaud, C., Freche, G., and Sicart, J. E.: Interception of snowfall by the trees is the main challenge for snowpack simulations under forests, in: Proceedings of the International Snow Science Workshop (ISSW), Innsbruck, Austria, pp. 705–710, available at: http://arc.lib.montana.edu/snow-science/objects/ISSW2018_O08.4.pdf (last access: 12 May 2020), 2018. a, b, c

Webster, C. and Jonas, T.: Influence of canopy shading and snow coverage on effective albedo in a snow-dominated evergreen needleleaf forest, Remote Sens. Environ., 214, 48–58, 2018. a

Yakir, B.: Nonparametric Tests: Kolmogorov-Smirnov and Peacock, chap. 6, pp. 103–124, John Wiley & Sons, Ltd, https://doi.org/10.1002/9781118720608.ch6, 2013. a