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² 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 km²; 46.78945^∘ N, 9.79632^∘ E), (b) in the Rocky Mountains of northern Utah, US (∼13 km²; 41.85046^∘ N, 111.52751^∘ W), and (c) in the southeastern French Alps at Col de Porte (∼0.01 km²; 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).

Download

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.

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.

Download

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⁻² (all returns) and a shot density of 19 points m⁻² (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.

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.

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}+\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}-\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⁻² returns and 17 last returns m⁻²) compared to the lidar flyover from the US (on average 7 returns m⁻² and 5 last returns m⁻²). 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⁻² 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⁻² (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 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.

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.