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

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


ception 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. (2015bMoeser et al. ( , 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) 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 Hydrol. Earth Syst. Sci., 24, 2545-2560, 2020 https://doi.org/10.5194/hess-24-2545-2020 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) 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 σ I 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(Moeser et al., , 2015a. 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 north-ern Utah, US, and from one site located at Col de Porte in the southeastern French Alps.

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

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. (2014Moeser et al. ( , 2015a. 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.

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 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 SNO-TEL 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).

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. Opensite 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.

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

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 I J 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 = −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).

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 . 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.

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 (maxmin)), 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].

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 = −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: (1 − F sky ) c σ c z 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 = 0.82 (±0.12), b = 0.0035 (±0.0036), and c = 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.
(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 = cc P HS with constant fit parameter cc = 0.40 (±0.03). as well as open-site snowstorm precipitation (size of symbols). Circles represent the development data set from Switzerland, symbols with a black border represent the validation data sets with squares for that from the US, and diamonds for data sets from France.

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 Hydrol. Earth Syst. Sci., 24, 2545-2560, 2020 https://doi.org/10.5194/hess-24-2545-2020 Table 1. Performance measures between measurement and parameterization of (I) spatial-mean snow depth interception I HS with (a) Eq. (2), (b) Eq. (3), and (c) a baseline model and of (II) standard deviation of snow depth interception σ I HS with (a) Eq. (4)

Model for standard deviation of snow depth interception
We parameterized the standard deviation of snow depth interception σ I HS by scaling P HS using the forest structure metric σ z . σ I 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). σ I HS normalized with P HS correlated much better with σ z (r = −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 σ I HS , we scaled a power-law function for P HS with the standard deviation of the DSM σ z in order to parameterize σ I HS as follows: (4) Constant parameters g = 0.78 (±0.10), h = 13.40 (±11.64), and j = 0.53 (±0.12) result from fitting a nonlinear regres-sion 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. σ I 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 σ I 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 σ I HS = jj P HS with constant fit parameter jj = 0.20 (±0.01).

Validation of model for standard deviation of snow depth interception
Overall, modeled and measured σ I 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 σ I HS is worse than that of our model performance (II(b) in Table 1). However, because one observed σ I 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 σ I HS for each location. . Measured and parameterized standard deviation of snow depth interception σ I HS at each site and for each storm date. Parameterized using Eq. (4) as a function of site means of standard deviation of the lidar DSM σ z (in color) as well as open-site snowstorm precipitation (size of symbols). Circles represent the development data set from Switzerland, symbols with a black border represent the validation data sets with squares for data from the US and diamonds for data from France.
The coefficient of variation (CV, description of variability) (CV I HS = σ I HS /I 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 val-ues for modeled I HS (9.3 cm), modeled σ I HS (3.7 cm), and modeled CV I HS (0.38) than measured I HS (9.2 cm), measured σ I HS (3.2 cm), and measured CV I HS (0.35). If the potential outlying data point from Utah is removed, the same overall modeled and measured mean CV I 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 σ I HS values (3.4 cm versus 3.3 cm).

Discussion
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 ∼ 1/(1 + exp(−k(P − P 0 ))) (e.g., Satterlund and Haupt, 1967;Schmidt and Gluns, 1991;Moeser et al., 2015b) or an exponential form ∼ (1 − exp(−k(P − P 0 ))) (e.g., Aston, 1993;Hedstrom and Pomeroy, 1998) with increasing precipitation. While the base function of Sat-terlund 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 σ I 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 σ I 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 σ I 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 pointbased 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 σ I 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 strati-fied 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 Hydrol. Earth Syst. Sci., 24, 2545-2560, 2020 https://doi.org/10.5194/hess-24-2545-2020 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 σ I 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 σ I 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.

Conclusion and outlook
The statistical models for spatial mean and standard deviation of snow interception presented here are a first step toward a more robust consideration of snow interception for various coarse-scale model applications. They were built upon a very large data set, validated by two other data sets from different geographic regions and snow climates, and performed well for all three sites and under differing weather conditions. For spatial mean interception all NRMSEs were ≤ 10 % and for the standard deviation of interception all NRMSEs were ≤ 13 %. Compared to a previous model for spatial mean SWE at 50 m × 50 m the presented models for spatial mean snow interception show improved model performances.
In our observed snow interception data sets, as much as 68 % and on average 43 % of the cumulative snowfall (accumulated snowfall of snowfall event in cm) was retained by coniferous forests (interception efficiency (snow interception/accumulated snowfall) of site means) and as much as 14 % and on average 11 % was retained by deciduous forests. These values compare well to previously observed snow interception in coniferous trees reaching up to 60 % of cumulative snowfall (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 highresolution DSM is available.
Data availability. The data that support the findings of this study are available from the corresponding author upon reasonable request (norahelbig@gmail.com).
Author contributions. 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.