the Creative Commons Attribution 4.0 License.

the Creative Commons Attribution 4.0 License.

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

### David Moeser

### Michaela Teich

### Laure Vincent

### Yves Lejeune

### Jean-Emmanuel Sicart

### Jean-Matthieu Monnet

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.

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 (Pomeroy and Schmidt, 1993; Pomeroy et al., 1998; Storck and Lettenmaier, 2002). In deciduous forests in the southern Andes as much as 24 % of the total annual snowfall may be retained (Huerta et al., 2019). Due to the sublimation of intercepted snow, a large portion of this snow never reaches the ground (Essery et al., 2003), and the interplay of interception and sublimation creates significant below-forest heterogeneity in snow accumulation. Rutter et al. (2009) 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 (see Essery et al., 2009; Rutter et al., 2009). 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 (e.g., Roesch et al., 2001; Bartlett and Verseghy, 2015; Webster and Jonas, 2018). 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 (e.g., Roesch et al., 2001; Roesch and Roeckner, 2006; Essery, 2013; Bartlett and Verseghy, 2015).

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 (Schmidt and Gluns, 1991; Hedstrom and Pomeroy, 1998; Bründl et al., 1999; Storck and Lettenmaier, 2002; Knowles et al., 2006; Suzuki and Nakai, 2008). 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 (e.g., Satterlund and Haupt, 1967; Schmidt and Gluns, 1991; Hedstrom and Pomeroy, 1998; Moeser et al., 2015b; Vincent et al., 2018). 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 Moeser et al. (2015b) 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 (*I*_{SWE}) have been suggested using a variety of canopy metrics and functional dependencies for the rate and amount of storm snowfall (e.g., Satterlund and Haupt, 1967; Schmidt and Gluns, 1991; Hedstrom and Pomeroy, 1998; Hellström, 2000; Lundberg et al., 2004; Andreadis et al., 2009; Moeser et al., 2015b; Huerta et al., 2019; Roth and Nolin, 2019). 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 (Moeser et al., 2014). Prior work has also used return density ratios of lidar, which is computationally faster but less accurate than SP (Morsdorf et al., 2006). 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 Moeser et al. (2015b), some forest structure metrics such as LAI and CC are highly cross correlated. Therefore, Moeser et al. (2015b, 2016) 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 (Moeser et al., 2014), 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 Roth and Nolin (2019) 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) (Essery et al., 2009; Rutter et al., 2009) 33 snow models were validated at individual forested and open sites and many models used the snow interception parameterization from Hedstrom and Pomeroy (1998). 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 (Dirmeyer et al., 2006), 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 (Krinner et al., 2018). Recently Huerta et al. (2019) 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 Hedstrom and Pomeroy (1998) model showing the overall best performance. Similarly, model simulations of Vincent et al. (2018) 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 Hedstrom and Pomeroy (1998) 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 (Roth and Nolin, 2019). While Roth and Nolin (2019) 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 (*I*_{HS} 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 (Moeser et al., 2014, 2015a, b). From a lidar DSM with elevations *z* (Moeser et al., 2014), 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 (*F*_{sky}), 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) (Helbig et al., 2009). 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.

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 m^{2} 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).

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 Moeser et al. (2014, 2015a, b). 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.1 ^{∘}C) to mild (−1.9 ^{∘}C) conditions.

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 (*F*_{sky}) 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 (PRISM Climate Group, 2012). 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 (Teich et al., 2019), 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.3 ^{∘}C) to mild (−1.4 ^{∘}C) conditions.

A 1 m resolution gridded lidar DSM was generated from a flyover in July 2009 and encompasses all field sites (Mahat and Tarboton, 2012; Teich and Tarboton, 2016) (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 *F*_{sky} 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 Morin et al., 2012; Lejeune et al., 2019). 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 ∼6 ^{∘}C 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) (Vincent et al., 2018) (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.9 ^{∘}C) 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 *F*_{sky} over the three 8 m transects.

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 *F*_{sky} 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 (Helbig et al., 2009). *F*_{sky} is a commonly applied model parameter when computing surface radiation balances and can be easily computed for large areas from DSMs. *F*_{sky} 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 (Helbig et al., 2009). We therefore selected *F*_{sky} to parameterize the site mean and standard deviation of snow interception (*I*_{HS}, *σ*_{HS}). Here, we compute *F*_{sky} 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 (Helbig et al., 2009) and were originally introduced to describe the radiant energy exchange between surfaces in thermal engineering (Siegel and Howell, 1978). Thereby, Helbig et al. (2009) solve the double area integral using uniform but adaptive area subdivision for surface patches *A*_{I}, *A*_{J}. *F*_{sky} for each surface patch *A*_{I} is one minus the sum over all *N* view factors *F*_{IJ} by assuming the sky is one large surface patch. *F*_{sky} is computed for each fine-scale grid cell of the DSM as follows:

Deriving *F*_{sky} 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 *F*_{sky} for forested regions are from sine and cosine weighted proportions of sky pixels of HP or SP as suggested, e.g., by Essery et al. (2008) or from LAI (e.g., Roesch et al., 2001). However, compared to computing *F*_{sky} on DSMs these methods rely on extensive field work.

The main advantage of deriving *F*_{sky} on DSMs is that *F*_{sky} can be derived spatially by averaging all fine-scale *F*_{sky} within a coarse grid cell. Here, we use the spatial mean of the sky view factor *F*_{sky} 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 *F*_{sky} (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 *F*_{sky} 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., *I*_{HS}, and not as snow water equivalent SWE, i.e., *I*_{SWE}. 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 (e.g., Hedstrom and Pomeroy, 1998; Moeser et al., 2015b; Roth and Nolin, 2019; Huerta et al., 2019) 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 (Hedstrom and Pomeroy, 1998). 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* (Yakir, 2013) for the pdf (nonparametric method) and compute the NRMSE for quantile–quantile plots (NRMSE_{quant}, normalized by the range of measured quantiles (max–min))) for probabilities with values in the range of [0.1,0.9].

## 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 (*I*_{HS}) by scaling open-site accumulated snowfall *P*_{HS} using the forest structure metrics *F*_{sky} and *σ*_{z}. From these three variables, the interception measurements of the development data set correlated best with *P*_{HS} (*r*=0.70). Snow interception efficiency (*I*_{HS}∕*P*_{HS}) correlations were slightly stronger for *σ*_{z} (*r*=0.71) than for *F*_{sky} ($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 (Roth and Nolin, 2019). To decide on the interception model form we considered previously commonly applied functional relationships with accumulated snowfall such as those from Hedstrom and Pomeroy (1998) and Moeser et al. (2015b), as well as simple relationships such as a power law. Together with our observed correlations of the forest structure metrics *F*_{sky} and *σ*_{z} with snow interception efficiency, we developed two statistical parameterizations for *I*_{HS} using two different base functions to scale *P*_{HS} with either *F*_{sky} and *σ*_{z} (Eq. 2) or with only *σ*_{z} (Eq. 3) as follows:

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

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; Rousseeuw et al., 2015). The 90 % confidence intervals of the parameters are given in parentheses. In both equations *P*_{HS} and *σ*_{z} are in centimeters.

The accuracy of a derived model between *I*_{HS} and *P*_{HS} 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 *I*_{HS} and *P*_{HS} was either a power law or a combination of a power law with an exponential dependence. Similar base functions were obtained for fine-scale *I*_{SWE} models by Moeser et al. (2015b) (exponential) and recently by Roth and Nolin (2019) (power law).

Estimated *I*_{HS} values from Eq. (2) or (3) increase with increasing *P*_{HS}, increasing *σ*_{z} or decreasing *F*_{sky}. This implies that with increasing forest density (i.e., larger *σ*_{z}), *I*_{HS} increases faster with increasing *P*_{HS}. Note that a lower *F*_{sky} 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 *P*_{HS} 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 Satterlund and Haupt, 1967; Moeser et al., 2015b). 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 *I*_{HS} and accumulated open-site storm snowfall *P*_{HS}. The second difference between both equations is that Eq. (2) uses both forest structure metrics (*F*_{sky} 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., *I*_{HS}= *c**c* *P*_{HS} with constant fit parameter *c**c*=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 *I*_{HS} models (Eqs. 2 and 3) were compared to the *I*_{HS} measurements from the development data set (Switzerland), as well as the *I*_{HS} 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 4 shows that there is good agreement between *I*_{HS} 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 NRMSE_{quant}<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).

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 *I*_{HS} and *P*_{HS}. The bias for the largest *P*_{HS} (US data set) is larger for the more compact parameterization (Eq. 3), whereas for the smallest *P*_{HS} (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 *P*_{HS} for the complex parameterization (Eq. 2) is −0.24 compared to 0.21 for the more compact parameterization (Eq. 3).

## 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 *P*_{HS} using the forest structure metric *σ*_{z}. ${\mathit{\sigma}}_{{I}_{\mathrm{HS}}}$ of the development data set correlated best with *P*_{HS} (*r*=0.82). The correlation with mean snow interception *I*_{HS} was less pronounced (*r*=0.33). ${\mathit{\sigma}}_{{I}_{\mathrm{HS}}}$ normalized with *P*_{HS} correlated much better with *σ*_{z} ($r=-\mathrm{0.68}$) than with *F*_{sky} (*r*=0.13).

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

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 *I*_{HS} from Eqs. (2) and (3). The 90 % confidence intervals of the parameters are given in parentheses. In Eq. (4) *P*_{HS} and *σ*_{z} are in centimeters.

${\mathit{\sigma}}_{{I}_{\mathrm{HS}}}$ derived from Eq. (4) increases with increasing *P*_{HS} or decreasing *σ*_{z}. This implies that with decreasing *σ*_{z} (decreasing forest density), the spatial variability in snow interception increases faster with increasing *P*_{HS}. The opposite correlation was found between *σ*_{z} and mean snow interception *I*_{HS}. 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 *I*_{HS} 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 *j**j*=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., NRMSE_{quant}<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 NRMSE_{quant} 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 *P*_{HS} 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 NRMSE_{quant} to 5 %).

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 *I*_{HS} 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 *I*_{HS} (9.3 cm), modeled ${\mathit{\sigma}}_{{I}_{\mathrm{HS}}}$ (3.7 cm), and modeled CV${}_{{I}_{\mathrm{HS}}}$ (0.38) than measured *I*_{HS} (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 *I*_{HS} (9.8 cm versus 9.9 cm) and modeled and measured ${\mathit{\sigma}}_{{I}_{\mathrm{HS}}}$ values (3.4 cm versus 3.3 cm).

We proposed two empirical models for spatial mean interception *I*_{HS} 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 *I*_{HS} and accumulated storm precipitation *P*_{HS} 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 Moeser et al. (2015b) for *P*_{HS} and is scaled by two forest metrics, namely the sky view factor *F*_{sky} in combination with *σ*_{z}. For both *I*_{HS} 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 *F*_{sky}. Though *F*_{sky} 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 *F*_{sky} over large scales and/or with fine resolutions is more computationally demanding than for *σ*_{z} (Helbig et al., 2009). 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öwe, 2014).

In general, more differences between the compact and more complex modeling approaches only displayed at the extremes. For instance, for small storm precipitation values (*P*_{HS}=3 cm), the more compact parameterization performs slightly better, whereas for very large storms (*P*_{HS}=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 *P*_{HS} differs from previous parameterizations for snow interception solely using the sigmoid growth $\sim \mathrm{1}/(\mathrm{1}+\mathrm{exp}(-k(P-{P}_{\mathrm{0}})\left)\right)$ (e.g., Satterlund and Haupt, 1967; Schmidt and Gluns, 1991; Moeser et al., 2015b) or an exponential form $\sim (\mathrm{1}-\mathrm{exp}(-k(P-{P}_{\mathrm{0}})\left)\right)$ (e.g., Aston, 1993; Hedstrom and Pomeroy, 1998) with increasing precipitation. While the base function of Satterlund and Haupt (1967) worked better for Moeser et al. (2015b), 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 (Moeser et al., 2016). For this reason, Mahat and Tarboton (2014) selected the relationship proposed by Hedstrom and Pomeroy (1998) for their parameterization of snow interception. However, the functional form of the Hedstrom and Pomeroy (1998) 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 *P*_{HS} and forest structure metrics, a simple power law was able to describe the spatial mean *P*_{HS} dependency of snow interception (cf. Eq. 3). The equation displayed that with increasing *P*_{HS}, *I*_{HS} increases. This is less pronounced with smaller *σ*_{z} or larger *F*_{sky} 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 (Roth and Nolin, 2019). 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 *P*_{HS} 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 (*F*_{sky} and *σ*_{z}) that integrates a power-law dependency of *P*_{HS}. 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 *F*_{sky} 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 Sicart et al. (2004) for solar transmissivity in forests under cloudless sky conditions. They assumed the fraction of solar radiation blocked by the canopy was equal to 1−*V*_{f} with *V*_{f} 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 *F*_{sky} 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 *I*_{HS} 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 *I*_{HS} 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 Roth and Nolin (2019). In addition to the spread in observed temperature conditions, our ranges of accumulated snowstorm *P*_{HS} values of the development data set are fairly broad (e.g., *P*_{HS} 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 (*P*_{HS}=3 cm) and one very large snowfall (*P*_{HS}=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 Roth and Nolin (2019) as well as the 50 m×50 m stratified SWE model (for 50 m×50 m grid cell size) from Moeser et al. (2015b). The Moeser et al. (2015b) model utilized the same Swiss data set as this study, which is currently the largest set of spatial interception measurements in the world. The Roth and Nolin (2019) 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 Nolin, 2019). 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 Moeser et al. (2015b) 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 Moeser et al. (2015b). 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 *I*_{HS} (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 (Roesch et al., 2001). 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 (Kaartinen et al., 2012; Eysn et al., 2015). Current available or scheduled countrywide data sets are now around 1–5 returns m^{−2} (e.g., Federal Office of Topography Swisstopo, 2019; Danish Geodata Agency, 2019; Latvian Geospatial Information Agency, 2019), 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 *F*_{sky} 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 Huerta et al. (2019) showed that previously published snow interception models developed for coniferous forests from Hedstrom and Pomeroy (1998), Lundberg et al. (2004), and Moeser et al. (2016) 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 *I*_{HS} and ${\mathit{\sigma}}_{{I}_{\mathrm{HS}}}$. For instance, using the more complex model for *I*_{HS} (Eq. 2) we obtained a mean bias of −6 cm, whereas when using the more compact model for *I*_{HS} (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.

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 (Pomeroy and Schmidt, 1993; Pomeroy et al., 1998; Storck and Lettenmaier, 2002) and up to 24 % of total annual snowfall in a deciduous forest in the southern Andes (Huerta et al., 2019).

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.

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

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.

The authors declare that they have no conflict of interest.

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.

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.

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

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

Danish Geodata Agency: https://download.kortforsyningen.dk/content/dhmpunktsky, last access: 22 November 2019. 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