Canopy Effects on Snow Accumulation: Observations from Lidar, Canonical-View Photos, and Continuous Ground Measurements from Sensor Networks

Abstract

A variety of canopy metrics were extracted from the snow-off airborne light detection and ranging (lidar) measurements over three study areas in the central and southern Sierra Nevada. Two of the sites, Providence and Wolverton, had wireless snow-depth sensors since 2008, with the third site, Pinecrest having sensors since 2014. At Wolverton and Pinecrest, images were captured and the sky-view factors were derived from hemispherical-view photos. We found the variation of snow accumulation across the landscape to be significantly related to canopy-cover conditions. Using a regularized regression model Elastic Net to model the normalized snow accumulation with canopy metrics as independent variables, we found that about $50\%$ of snow accumulation variability at each site can be explained by the canopy metrics from lidar.

1. Introduction

The snowpack in California’s Sierra Nevada has long served as the primary water resource for agricultural and urban uses [1]. For seasonal forecasts following the onset of snow melt, the estimation methods are turning from statistical estimates that use historical records to spatio-temporal water-balance estimates with integrated time series and spatial data [2,3]. Quantifying the spatio-temporal distribution of snow accumulation enables more-accurate forecasts of snow melt and streamflow, especially in a warming climate as statistical relations broke down; and it is also a long-standing challenge in snow hydrology [4,5]. In the high Sierra, orographic effects drive solid-phase precipitation falling over mid-to-high elevations [6,7]. During the snow-accumulation period, vegetation intercepts snowfall, causing the snowpack to distribute unevenly under the canopy. In boreal forests, as much as 60% of cumulative snowfall may be intercepted by forest canopy in mid-winter, and annual sublimation losses can be 30–40% of annual snowfall [8,9]. In mixed-conifer forest, both observations and modeling results have shown that canopies can reduce 10–25% of total snow accumulation in surrounding open areas; and the reduction can be 30–40% for snow accumulation under canopy [10]. The canopy effects are also found to vary with local landscapes [11]. Accurately quantifying canopy interception of snowfall is the foundation to estimate the total accumulation and snow melt with higher accuracy and precision during the spring season.

The canopy interception of snowfall can be quantified as the snow-storage capacity of the canopy and interception efficiency (interception/snowfall). The snow-storage capacity is the maximum amount of snowfall that can be intercepted by the canopy. It is determined by the leaf area, tree species, and initial canopy snow load [12]. The interception efficiency is found to decrease with increasing snowfall, initial canopy snow load and temperature. It has been observed to increase with leaf-area index (LAI) and canopy coverage [8,13].

The coniferous canopies interception of snowfall is challenging to directly measure and quantify, which has led to interception being inferred from snowfall patterns on the ground, as noted above. Previous studies designed special weighing devices such that the weight of the intercepted snow and accumulated snow can be measured at the same time. The total snow interception is found to be correlated with the season or event snowfall [8,12]. Similar findings were observed from simulating the snow evolution by using energy balance models with canopy metrics [14,15,16]. Thus, several process models have incorporated this statistical finding to estimate canopy-cover effects on snow accumulation [17,18,19,20].

A common approach to calculate canopy interception uses canopy metrics that are highly correlated with snow accumulation. Retrieving canopy metrics has advanced in recent years. The technology has been advancing from the traditional plant-canopy analyzer [18,21,22,23,24], to hemispherical-view camera [25,26], and recently, to lidar [27,28]. The plant canopy analyzer was commonly used for retrieving the LAI in the forest. By using the hemispherical-view camera, the pixels of an image can be classified as either canopy-covered or clear; and thus the percentage of clear view for each zenith angle can be quantified as the sky-view factor, which was also found to be a statistically significant predictors for parameterizing snowfall interception in process models [25,29]. The point-cloud data collected using lidar can be used for reconstructing the 3-dimensional canopy structures if the point-cloud is sufficiently dense. Algorithms have been developed for deriving LAI from lidar point clouds, and additional canopy metrics from lidar often get derived for quantifying the snowfall interception [30].

In addition to canopy-metric retrieval from lidar, canopy effects can also be quantified using statistical models, with dense spatial measurements of snow depth or snow water equivalent (SWE) [28,29]. Most previous studies were conducted using lidar measurements, either airborne or terrestrial. Both the airborne and terrestrial lidar can provide dense spatial snow-depth measurements (>10 $\mathrm{pts}/\mathrm{m}$). With the extensive footprint provided by airborne lidar scans, the canopy effect on snowpack spatial distribution can be quantified with large samples. Terrestrial lidar has a much smaller footprint compared to airborne lidar [31], however, it is able to provide multiple scans per season. Thus the temporal variation in canopy effects can also be determined.

One short coming in using lidar is that it can lack temporal completeness, especially during the precipitation season, when cloud cover limits taking measurements. Lidar also requires clear-sky conditions to take measurements to prevent the laser pulse intensity from attenuating because of rain drops and snow flakes [32]. A dense cluster of snow-depth sensors can compensate the weakness of lidar in terms of temporal consistency. Combining the vegetation structures derived from lidar measurements and continuous snow-depth measurements, there is potential for the spatial variation of snow accumulation to be accurately quantified. In our study, we used long-term spatially dense snow measurements in the Sierra Nevada, together with the lidar-derived canopy metrics, to study the canopy effect on seasonal snow accumulation.

The overall aim of the work reported here is to explain and predict the spatial variability of snow accumulation using lidar-derived canopy metrics and clustered snow-depth sensor measurements. We address two major question. First, to what extent can one use lidar-derived canopy metrics to predict the snow accumulation spatially? Second, what is the relative importance of various lidar-derived canopy metrics?

2. Methods

2.1. Study Areas and Snow-Depth Sensor Data

The study was conducted over three areas in California’s central and southern Sierra Nevada (Figure 1a). For each study area, snow-depth sensors (Judd Communications) were deployed in clusters (Figure 1b,c;

Table 1. Elevation information of each site and time-series data availability.

Site	Sub-Site	Elevation, m	Data Availability, Water-Year a
Pinecrest	Upper	1808–1834	2014–2017
	Lower	1748–1778	2014–2017
Providence	Upper	1975–1984	2008–2016
	Lower	1730–1740	2008–2016
Wolverton	Site1	2225–2227	2008–2016
	Site2	2250–2266	2008–2016
	Site3	2590–2602	2008–2016
	Site4	2630–2648	2008–2016

^a Water year begins on 1 October of the previous calendar year, e.g., water year 2014 begins on 1 October 2013.

), with topographic characteristics (elevation, aspect) varying between clusters and canopy-cover conditions varying within each cluster. The snow-depth sensors have an accuracy of ±1 $\mathrm{cm}$, with measurement range of 0.5 to 10 m and beam width of 22 degrees. All study sites encounter wet winter for most of the years because of their proximity to the Pacific. Due to the orographic effect, most precipitation falls as rain at elevations below 1500 m and as snow above 2000 m [28]. The date of peak snow-accumulation is commonly recognized as 1 April [1]. Pinecrest, located in the Stanislaus-Tuolumne Experimental Forest, is the lowest in elevation and is also relatively flat. The lower site of Providence has a similar elevation range as Pinecrest, and the upper site is 200 m higher. Wolverton, the highest of the three study areas, receives mainly snow, whereas Providence and Pinecrest receive precipitation as both rain and snow. There were 12 sensors installed at Pinecrest in 2014, 27 sensors at Providence, and 25 sensors at Wolverton, installed in 2008. All time-series data are shown in Figure S1.

2.2. Lidar Data

The lidar survey was performed in August 2010. The flight parameters and sensor settings are described in

Table 2. Flight parameters and sensor settings ^a.

Flight Parameters		Equipment Settings
Flight altitude	600 m	Wavelength	1047 nm
Flight speed	65 $\mathrm{m}{\mathrm{s}}^{-1}$	Beam divergence	0.25 mrad
Swath width	233.26 m	Laser PRF	100 kHz
Swath overlap	50%	Scan Frequency	55 Hz
Point density	10.27 ${\mathrm{m}}^{-2}$	Scan angle	$\pm {14}^{°}$
Cross-track resolution	0.233 m	Scan cutoff	3${}^{°}$
Down-track resolution	0.418 m	Scan offset	$0^{°}$

^a Lidar-survey flights were operated by NCALM, and point-cloud data were produced based on their processing software.

. The point-cloud lidar data were used for generating raster data sets. The raw point-cloud files, we divided them into 250 × 250-m tiles using LAStools lidar-processing software [33]. For all data points in the point-cloud data set, each one of them has a classification attribute and the data point was a laser return from the ground if the classification attribute equals 2. We extracted the ground points from each tile and interpolated them into a 0.5-m resolution digital elevation model (DEM) using a simple kriging model with a spherical covariance function [34]. The 250 × 250-m DEM tiles were mosaic-ed together to form a single DEM of the study area. A digital surface model (DSM) was generated from all first returns of the lidar point cloud. Subtracting the DEM from the DSM produces the canopy-height model (CHM). Individual trees were segmented out from the CHM using a watershed-segmentation algorithm implemented in SAGA GIS software [35]. The segmentation was conducted after applying a Gaussian filter with a 3 × 3 moving window to fill pits in CHM and suppress irrelevant local maximas [36,37]. Over each snow sensor location, canopy metrics, including mean canopy height, standard deviation of canopy height, and canopy-cover fractions, were extracted at search radii from 2 to 40 m at 1-m increments. These canopy metrics were calculated by masking out the circular area with certain radius around the snow-depth sensor location. The canopy-height statistics were calculated with the CHM pixels that are within the circular area. The canopy-cover fractions were determined by calculating the percentage of canopy pixels within the circular area. The canopy pixel is classified from the CHM if a pixel value is greater than 2 m [28]. The distance from the sensor location to the closest tree trunk was also calculated by using the segmented-tree information. In total we derived more than 100 canopy metrics from the above methods.

2.3. Canonical-View Images

The canonical-view images, are the hemispherical photos taken at each snow-depth sensor. The sky-view factors (f, SVF) at each individual zenith angle $\theta$ were derived from the raw image, with each pixel classified into binary representation of sky and non-sky elements. Within each zenith angle, the SVF was estimated by the ratio of the number of sky pixels and the number of total pixels. An example of the canonical-view image and the estimated SVF is shown in Figure 2. The total sky-view factor at each sensor node was also estimated using the equation below [29,38],

$$f_{tot}=\frac{{\int }_{\theta =1}^{\theta =90}sin\left(\theta \right)f\left(\theta \right)}{{\int }_{\theta =1}^{\theta =90}sin\left(\theta \right)}$$

(1)

The sky-view-factor data are available for Wolverton and Pinecrest. Considering that only 6 sensor locations at Pinecrest have the images, the data available are not enough to build robust statistical models, so Pinecrest is not included in analysis with the sky-view factors. We included the total sky-view factor and the sky-view factors at each zenith angle as independent variables for modeling the snow accumulation observed at each sensor location. The results are compared with the modeled results that use lidar-derived canopy metrics as predictors.

2.4. Snow Accumulation Events Detection

The data availability over time for each site is shown in Table 1. We studied the canopy effect on snow accumulation using the following procedure (Figure 3).

• Get the moving average of each snow-depth time series with a window size of 2 days. Then calculate the 1st order gradient of the time series. This made estimates less vulnerable to high-frequency noise in the snow-depth data.
• The 1st-order gradients over all sensors, were used to calculate the $x\%$ quantile of the gradient. The quantile statistic was then compared with a pre-configured threshold to determine if most sensors observed snow accumulation. Neighboring accumulating days were then grouped together to form a single event.
• For snow-accumulation event detection, we set the quantile for snow accumulation as 30%. It means that if 30% of sensors show an ascending trend in one day, we can classify this day as an accumulation day.
• The daily gradient thresholds were also need to be optimized, along with the gap length between two adjacent snow accumulation dates. The optimized threshold for snow accumulation events is 0.1 cm. If two snow accumulation events were temporally close, we used the following rule to determine if the two neighboring events can be merged together or not.

For snow-accumulation events, the optimized way to combine two neighboring events is to first judge if the length of gaps between two events are shorter than one third of the sum of length of two events. Then, if the snow depth data for most sensors does not show a descending trend during the gap period, the two close events were combined into one.

The parameters that were used in the above steps were derived from a sensitivity analysis. The analysis was performed based on the snow-depth time series collected from Wolverton in water year 2010. We manually extracted the snow-accumulation periods from the multivariate time-series data and used the extracted periods as the standard to be compared with. The analysis was conducted in a Monte Carlo manner, with parameters randomly drawn from their predefined sets. The algorithm classifies each day of the water year as either an accumulation day or non-accumulation day. Comparing the results produced from the algorithm to the manually extracted periods, we calculated precision and recall as the classification accuracy for the entire year. We used these two metrics to determine the near-optimal parameter set that was used.

2.5. Statistical Analysis

All extracted accumulation events were used for statistical modeling, with features derived from the lidar data and the sky-view factors derived from the canonical-view camera images. We conducted regression analysis to study if canopy metrics can be used as predictors for estimating snow accumulation at various canopy-covered conditions. For each individual accumulation event, the total snow accumulation at each sensor node was estimated as $\Delta H=H_k-H_0$ where $H_k$ is the snow-depth at the last time step and $H_0$ is the snow depth when the precipitation event started. Considering topographic effects on precipitation along the elevation gradient, we offset the total solid precipitation for each individual event at each site using elevation. For all study sites, the sensors were clustered in a lower-elevation site and a upper-elevation site. We calculated the mean snow accumulation at both clusters and we subtracted the absolute difference of the two averages from the snow accumulations at the upper site to derive the offset accumulations. The offset results were standardized to the range of 0–1. The detrended target values are regressed using Elastic Net (implemented in Python Scikit-Learn [39]), which is a regularized regression method that linearly combines both $L1$ and $L2$ penalties in the regression model. We assumed a linear regression problem defined as:

$$y=X\beta +\epsilon$$

(2)

where y is the target value and X is the matrix of all covariates. The regression coefficient $\widehat{\beta }$ can be estimated with:

$$\widehat{\beta }=\underset{\beta }{arg\; min}{(\parallel y-X\beta \parallel }^2+{\lambda }_2{\parallel \beta \parallel }^2+{\lambda }_1{\parallel \beta \parallel }_1)$$

(3)

where ${\lambda }_1$ and ${\lambda }_2$ are the weightings for the $L1$ and $L2$ regularizations, respectively. The Elastic Net was chosen over other regularized regression approaches for its ability to address correlated covariates and high numbers of covariates [40,41]. We have also considered nonlinear and ensemble models, such as the Random Forest regression, but the performance was not as good as the Elastic Net regression (results not shown). In our case, the canopy metrics can be highly correlated when the search radii are close and the number of covariates included in our analysis is more than 100.

In order to have representative estimates of how much variability can be explained by the Elastic Net model we used bootstrap to resample the data for 20 iterations and estimated the cross-validated coefficient of determination ($R^2$) within each iteration. The distribution of the $R^2$ was estimated from multiple bootstrapping results.

We also applied correlation analysis to explore the most-informative radius of lidar-derived canopy features and the most-informative zenith angle of the sky-view factors from the canonical view images. We computed the correlation between the snow accumulation from each individual event and the lidar-derived mean canopy height at various search radii and at various zenith angles. The correlation coefficients (R) were compared at various radii and angles to select the optimal radius and zenith angle. Considering that Pinecrest has a relatively short record, most of which is during a recent drought, we did not conduct the analysis for Pinecrest. Also, camera images are not available for Providence, thus we only conducted radius dependency analysis at that site. At Wolverton, we selected a few near-optimal searching radii and zenith angles. We used these selected variables and conducted a step-wise linear regression with forward selection [42] for exploring the relative importance between variables. We used step-wise linear regression rather than other embedded feature selection approaches in order to show the results from s perspective of correlations.

3. Results

3.1. Snow Accumulation Events Extracted from Snow-Depth Time Series

Before running the snow-accumulation detection algorithm on all time series, we performed a sensitivity analysis on the 2010 time-series data from Wolverton. We found that as long as the gap length is about $1/3$ of the total length of the neighboring events, with the snow-accumulation rate threshold being more than 0.1 cm and less than 1 cm, and the quantile of the snow-accumulation sensors being in between 30% to 50%, both precision and recall are higher than other parameter sets (see Figure S2 for analysis results.)

We applied our snow-accumulation events-extraction algorithm on all snow-depth sensors for all time periods when cleaned snow-depth data were available. The performance of the detection algorithm is similar to manual extraction, done previously for some years (data not shown). As shown in Figure 4, the algorithm is able to detect most major snow-accumulation periods. A summary of accumulation events for each site is shown in Table S1 and the distribution of the magnitude of the accumulation at each study area is shown in Figure S3. Wolverton has more events because it is at higher elevations. At the time of this analysis, 2016–17 data were not available for Wolverton.

3.2. Statistical Modeling Results

The variability that the Elastic Net model can explain over the three sites is shown as in Figure 5. The uncertainty range of the variability that can be explained by the Elastic Net model is much larger for the Pinecrest (0.55 on average and $\pm 0.1$ of variation) analysis than the other two areas (0.45 on average at Providence and 0.25 on average at Wolverton, and $\pm 0.05$ of variation for both). Note that Pinecrest is the lowest and Wolverton the highest elevation of the 3 sites.

At Providence and Wolverton, excluding the minor accumulation events (≤15 cm) can increase the variability that can be explained by the Elastic Net model, with more than 50% explained at Providence and 40–50% explained at Wolverton (Figure 6a,b). Due to the fact that for minor accumulation events the signal strengths are not greater than the uncertainty range of the snow-depth sensors, including these data points will degrade the performance of the Elastic Net model. At Wolverton, the spatial variability of snow accumulation that can be explained reaches the maximum when the mean snow accumulation is between 15 $\mathrm{cm}$ and 30 $\mathrm{cm}$ (Figure 6b). At Pinecrest, no particular trends can be observed, as the number of data points is limited. When including most of the data points, the variability explained stabilized around 40–60% (Figure 6c).

For Wolverton, the only study area that both SVF and lidar are available, the trends in Figure 6b show that $R^2$ is the maximum when the threshold is above 15 cm. We thus filtered the data with mean precipitation that are below 15 cm. We conducted three sets of analysis, including using lidar-derived canopy metrics as the predictors, using SVF as the predictors, and using both lidar and SVF as the predictors in the Elastic Net model. As shown in Figure 7, the coefficient of determination calculated from 20 bootstrapping runs of predicting the snow accumulation suggests that the trained models are more accurate in predicting the total snow accumulation at the unobserved sensor locations if using lidar-derived variables because the variability of the $R^2$ is larger than that using sky-view factor. Also, the third box-plot of this figure suggests that using both lidar-derived variables and sky-view factors is only as good as just using lidar, indicating that the sky-view factors do not provide more information than lidar in terms of predicting snow accumulation.

The canopy-related variables derived from lidar and canonical-view images are compared based on the correlation coefficients between the total solid precipitation and the canopy attributes over each individual event across sensors, as well as the step-wise linear regression. The correlation analysis (Figure 8a) shows that the surrounding canopies have a stronger effect on the snow accumulation on the ground than does the canopy right above. The canopy mean height within a 15-m radius at Providence is the most effective distance, while the optimal radius is about 8 m at Wolverton. For sky-view factor, the optimal zenith angle is about ${21}^{°}$ at Wolverton. In Figure 8a, we identify each individual precipitation event by the transparency of each curve, from which we can see that heavier storms have more dominant weights for characterizing the canopy effects at different search radii from lidar data and zenith angles from canonical-view photos. In addition, the step-wise regression analysis (Figure 8b) conducted on the selected optimal variables of lidar-derived and canonical-view imagery features suggest that the canonical-view imagery features are more important and the marginal information that lidar provides is limited comparing to the first canonical-view imagery feature selected.

The canopy-cover information is slightly more important than canopy-height metrics, which include tree heights and tree-height standard deviations. This was verified by correlation coefficients in the regression analysis between snow accumulation and both tree height at incremental search radii and SVF at incremental zenith angles. In general, SVFs are more correlated with snow accumulation than is tree height. The step-wise regression analysis also suggests that the sky-view factor at the optimal zenith angle, is more important than is tree height at the optimal search radius. Although the tree height is an important metric characterizing trees in the forest, it does not necessarily represent the density and interception capacity of the canopy. Even though sky-view factor only represents the canopy-cover condition at the lowest layer of canopy, it is still a direct index representing the interception capacity of part of the tree crown. In contrast, the tree height does not necessarily represent the interception capacity, which is the reason that the lidar-derived variables are not as important as the sky-view factors. Please note that having a few important variables does not guarantee more-accurate estimates in snow accumulation. Combining all lidar-derived variables together, we have made more-accurate estimates of snow accumulation than using the sky-view factors (Figure 7).

In addition, we compared the correlation coefficient between different types of lidar-derived canopy-related features and the snow accumulation over different sensor nodes. The features include mean canopy height over the search radii, standard deviation of the height, maximum canopy height, and canopy coverage. As is shown in Figure 9, the amount of data at Pinecrest is not enough to draw solid conclusions. At Providence and Wolverton the correlation coefficient is a concave shaped function of both canopy-height mean and canopy coverage at various searching radii. The maximum canopy height at the smallest search radius correlates the most with the snow accumulation. The standard deviations of the canopy heights at various search radii show contrast trends at Providence and Wolverton.

4. Discussion

4.1. Canopy Effect at Different Elevations

Among the three sites studied, the variability of snow accumulation that the canopy-related variables can explain varies from site to site and also depends on the mean cumulative precipitation over the entire event. The difference between sites can be attributed to different elevations, as the canopy-cover density decreases as elevation increases and the solid-phase precipitation increases with elevation. For example, the instrumentation locations at Pinecrest and Providence are at much lower elevations compared to Wolverton. About 40% of snow-accumulation variability is attributable to the canopy effects at these two sites, however only 25% can be explained over Wolverton. This suggests that at higher elevations, where precipitations is heavier, the canopy effects can be diluted by the heavy snowfall, which is similar to previous findings [13], where total interception of snowfall saturated when the total precipitation reached certain thresholds for different tree species. In addition, we observed some noise introduced by the low-precipitation events in the regression analysis [11]. Figure 6 suggests that the spatial variability of precipitation is less explainable by the canopy-related variables when the total precipitation is small.

The difference of the canopy effect on snow accumulation at the lower elevations is also affected by the rain-on-snow events. Analyzing rain-on-snow events would need more meteorological data, especially temperature and relative humidity, for estimating the rain-snow partition along the elevation gradient [43]. A rain event can accelerate snowmelt by altering the energy balance of the snowpack [44,45]. As such, the snow accumulation during mixed rain and snow events reflects the net of snow precipitation and snowmelt. The spatial patterns of total snow accumulation at each location need to be analyzed with a full surface-energy-balance model such that canopy effects on accumulations can be described.

4.2. Forest Thinning Effects on Snow Accumulations

At Pinecrest, considering that not all sensor locations have sky-view factors, we used multi-year snow-depth data from surveys and the sky-view factors for all survey locations, comparing correlation patterns observed from this site with Wolverton. Similar to Figure 8a, in Figure 10 we see that the adjusted $R^2$ peaked around a zenith angle of ${30}^{°}$. Over the snow survey locations, the forest above has been thinned the same amount, but with either even density or heterogeneous density. It seems that even thinning has a greater impact on canopy effects, as the correlation coefficient is close to zero for most zenith angles except when the angle approaches 90 degrees. The thinning with variable forest density holds a similar pattern comparing to the control group, with adjusted $R^2$ peaked around a zenith angle of 30–${50}^{°}$. But the magnitudes of correlations are smaller during the dry years (2013–2014).

Comparing lidar-derived canopy metrics at increment radii, Figure 9 suggests that the most important canopy structures may not be the canopy layers right above the measured locations. The canopy surrounding within a few meters could be even more important, as the interception capacity can be larger when the trees are clustered together versus for a single tree. Our findings are aligned with what was found previously [11,12]: multi-layered forest cover increases canopy-interception efficiency, resulting in significant reductions in subcanopy snow accumulation. Even thinning will make the trees evenly spaced [46], and as such the number of forest layers will be greatly reduced, with commensurate reduction in canopy interception of snow. This is the reason that the snow depth is not correlated with the SVF at any zenith angles in the evenly thinned forest at Pinecrest, while the control group and the forest with variable density still retain multi-layered canopies.

4.3. Potentials to Extend the Analysis

The integrated sensing system, including lidar, canonical-view images, and ground measurements, allowed us to study canopy effects on snow accumulations. The snow-accumulation detection algorithm can be applied to most mountainous areas that have seasonal snow, with minor tuning of parameters to fit the local snow-accumulation patterns. However, the current statistical analysis is limited to only canopy effects, and did not address the interactions between canopy and terrain or canopy and meteorological conditions. Although the current statistical model trained with canopy metrics can explain about 50% of the snow-accumulation variability, more may be explained when local topography and meteorology are considered. For canopy-terrain interactions, a spatially denser deployment of snow-depth sensors is needed in order to build relationships between terrain and snow that are statistically significant. For meteorological data, other than air temperature and relative humidity mentioned in Section 4.1, wind distribution is also important as it affects canopy interception, snow unloading from the canopy, and snow redistribution in the forest [8,12]. We did not include wind in our analysis, as there is limited wind monitoring on site and the wind effect is negligible in the Sierra Nevada [30], but wind can be important for regions like the Alps and the Rocky Mountains.

5. Conclusions

We found correlations between the lidar-derived canopy attributes and snow accumulation extracted from the multi-year time-series snow-depth measurements. The correlation is stronger when the precipitation event has more than 15 cm of snow accumulation. And the correlation is also much stronger at a lower elevation (<2000 m) because of denser vegetation. Although the lidar-derived canopy attributes are complementary to sky-view factors in explaining the snow-accumulation variability, the sky-view factors at certain zenith angles are more important than lidar-derived variables in terms of correlation. According to the correlation analysis using lidar-derived canopy metrics at various search radii, the canopy surrounding the snow surface within 8-m radius is more important than canopy structures within either smaller radius or larger radius, indicating that a clustered canopy effect is stronger than is canopy of a single tree. The above findings suggest great potential of using lidar and ground measurements for studying canopy effects on mountain snowpack.