Toward a novel laser-based approach for validating snow interception estimates

,


Introduction
The hydrology of snow dominated forests is controlled by interactions of mass and energy fluxes between snow and forest structural elements. As forest cover increases, snow accumulation on the ground is typically reduced because of canopy snowfall interception and subsequent sublimation, which can account for as much as 60% of the cumulative snowfall depending on forest type, duration of snow storage in the canopy, and seasonal hydrometeorological conditions (Hedstrom & Pomeroy 1998;Molotch et al., 2007). The sensitive connection between forest structure and snow interception therefore has important implications for the hydrology in any region around the globe where the major proportion of total water input comes from snow. Understanding this relationship is increasingly important with widespread observed and projected shifts from snow to rain (Klos et al., 2014), changes in the frequency of winter rain-on-snow events (Floyd & Weiler, 2008;Musselman et al., 2018), and changes in forest vegetation due to fire (Westerling, 2016), drought (Allen et al., 2010), insects (Bentz et al., 2010;Frank, Massman, Ewers, & Williams, 2019) and other disturbance processes that might be altered by a changing climate and/or forest management.
While the importance of snow interception has long been acknowledged, it is also difficult to measure, map, and model. Direct measurement has typically been limited to resource-intensive point measurements derived from weighed trees, which are generally limited to small or medium trees (Hedstrom & Pomeroy, 1998;Knowles, Dettinger, & Cayan, 2006;Storck et al., 2002;Suzuki & Nakai, 2008) or tree branches (Brundl, Bartelt, Schneebeli, & Fluhler, 1999;Schmidt & Gluns, 1991). Indirect measurements have compared snow accumulation between forested sites and nearby clearings. Although indirect measurements have advantages (e.g. estimating spatial variance), the accuracy has long been questioned (e.g., Miller 1964) and the measurements can be confounded by wind redistribution of snow, branch unloading, and reference site size (Moeser, Stähli, & Jonas, 2015).
A potential novel method could use terrestrial lidar (TLS) because three-dimensional lidar point clouds, based on the laser return-time/distance relationship, can be generated for any size tree. The point clouds can be transformed into an α-shape (convex hull with a polyhedral surface approximating the shape and volume the tree), from which volumes can be calculated (Edelsbrunner, 1994;Lafarge, Pateiro-Lopez, Possol & Dunkers, 2014). The α-shape convexity (i.e. α-convexity) (Rodriguez-Casal, 2007) is tunable according to TLS spatial resolution and study objectives -be it calculating volumes of fine-scale structures like branches or smoothed-over structures like snow-laden trees. Intercepted snow volume can be estimated by subtracting snow-free tree volume from snow-on tree volume. Furthermore, novel autonomous terrestrial laser scanning (ATLS) systems (Adams et al., 2014;Eitel et al., 2013) could enable timeseries characterization of seasonal dynamics associated with snow interception and unloading. Ultimately, TLS based snow interception may have a number of advantages over traditional methods, including: spatially explicit estimation of snow interception for different aspects or portions of tree canopies, data collection in difficult to access terrain (Adams, Bauer, & Paar, 2014) known to be important contributors to water budgets (Bühler, Adams, Bösch, & Stoffel, 2016;Hood & Hayashi, 2010), and providing time-efficient data for calibration of emerging aerial lidar (ALS)-based snow interception models to specific forest types (Deems et al., 2013;Moeser et al., 2015;Painter et al., 2016).
The objectives of this study were twofold: (1) Test the feasibility of using ATLS to directly map intercepted snow mass and (2) Determine the effect of α-convexity on ATLS-derived snow mass estimates. In doing so, this study provides a preliminary feasibility assessment for estimating snow interception mass with terrestrial laser scanning, providing information on challenges and opportunities for future research.

Study Site
Two model hanging trees measuring 1.83 meters (m) in height (hereafter referred to "left tree" and "right tree" -see Figure 1) were installed prior to winter 2017 following established approaches outlined in Hedstrom and Pomeroy, 1998. The trees were off-the-shelf Christmas trees, bilaterally symmetrical and not representing a specific species. Artificial trees were utilized to avoid desiccation and interception estimates that may be affected by progressive needle drop. The trees were installed at the University of Idaho McCall Field Campus (44.9353472°, -116.0820167°) in the mountains of west-central Idaho, which receives an average of 3.4 m of snow fall (depth) per year at 1528 m elevation (Western Regional Climate Center, 2016).

Field Measurements
Load cells measured strain gauge output (mV/V), an electrical signal which is proportional to the applied excitation voltage, from the hanging trees in one-minute intervals. Known masses were hung from the load cells to verify measurement accuracy and to develop a calibration equation (mass = 30.18 * mV/V -4.6917) which converted strain gauge output to kg. The originating mass of a snow-free tree at the beginning of each winter was subtracted from subsequent measurements to calculate snow masses for each scan.
An ATLS scanned one side of the trees at a distance of 6.2 m and produced two high resolution point clouds per day (1.12 cm spot size and 0.20 cm point spacing at 10 m) ( Figure 1b). The ATLS employs a rugged time-of-flight laser rangefinder (optoNCDT ILR 1191 with 905 nm near infrared laser and 1.7 mrad beam divergence; Micro-Epsilon Messtechnik GmbH & Co. KG, Ortenburg, Germany) designed for harsh environments (see Eitel et al., 2013 for more detail). The ATLS completed one scan in 13 hours.
Assuming each tree canopy to be bilaterally symmetrical and circular in shape, and given a known canopy diameter, distance from scanner to canopy perimeter, and location of the ATLS, trigonometric calculations yielded 47.7% for the percent of tree canopy perimeter "seen" by the ATLS. Rounded to 50% for the analyses, snow masses obtained from the load cells were therefore divided by two and averaged across each ATLS scan duration to allow for comparison with ATLS data.

Lidar volume estimates
The ATLS point clouds were transformed into "alpha-shape" structures in R using the alphashape3d R package (Lafarge, Pateiro-Lopez, Possol and Dunkers, 2014). An alpha-shape structure is a three-dimensional convex hull generated from a point cloud, approximating the shape and volume of the scanned structure. The convex hull is fitted with Delauney Triangulation (drawing triangles between points so that there is no overlap between triangles). The α-convexity parameter (α) (Edelsbrunner & Mucke, 1994) is selected by the user, corresponding to data resolution and units of the input data. An α-convexity parameter of 1 corresponds to the convex hull; the closer α is to zero, the more triangle borders are deleted to make a better fitting, flexible and concave hull (see Figure 3). A sensitivity analysis was conducted by comparing measured vs. ATLS derived snow mass (see section 2.5) with α = 0.010, 0.015, 0.020, 0.025, and 0.030 m -a range deemed appropriate given the ATLS spatial resolution (see section 2.2). This permitted the selection of α that resulted in the best agreement between observed and TLS based estimates of intercepted snow mass. The originating volume of a snow-free tree at the beginning of each winter was subtracted from subsequent measurements to calculate snow volumes for each scan.

Statistical analysis
The precision and accuracy of laser derived estimates of snow interception (objective 1) was determined by fitting a simple linear regression model in R (R Core Team, 2013) between ATLS derived intercepted snow mass in kg (independent variable) and load cell derived intercepted snow mass in kg (dependent variable) (Piñeiro, Perelman, Guerschman, & Paruelo, 2008). Following common practice when using convex hull approaches (e.g., Lafarge, Pateiro-Lopez, Possol & Dunkers, 2014), a sensitivity analysis to determine the effect of α-convexity on ATLS derived snow mass estimates (objective 2) was performed by comparing R2 (goodness-of-fit), Root Mean Square Error (RMSE) in kg, and regression intercept/slope estimates (indicative of model bias) for each model combining alternative α levels (see section 2.3) and snow density estimation methods (see section 2.4). The optimal model was determined by selecting, first, the best performing snow density estimation method and, secondly, the α level that produced the least under/over estimation in ATLS based mass predictions (i.e. closest to 1:1 line) and lowest RMSE.

Results
A total of 115 complete ATLS scans were recorded for the left tree in the first winter between January and April; a total of 69 complete ATLS scans were recorded for the left tree in the second winter between November and March. A total of 83 and 69 scans were recorded for the right tree over the same time periods.
Discrepancies in sample size between the left and right tree were related to incomplete ATLS scans in which scanner malfunction truncated a portion of the scene(e.g., see Figure 1b).
Data exploration using results from the left tree revealed that measured snow interception (averaged across multiple complete scans) in the first winter averaged 3.19 kg (1.0% of season total), with a standard deviation of 4.59 kg and a maximum of 16.95 kg (4.8% of season total). During the second winter, measured snow interception for the left tree averaged 3.31 kg (2.3% of season total), with a standard deviation of 2.76 kg and a maximum of 10.19 kg (4.5% of season total). Estimated mean fresh snow densities using equations [1][2][3][4] for the left tree in both winters are summarized in Table 1. Density values are similar to Mair et al. (2016), in which estimates produced by equation [2] approximated the Brazenec (2005) constant, estimates produced by equation [3] were higher than the constant, and estimates produced by equation [4] displayed a wider range.
Sensitivity analyses (Table 2) using data spanning both winters for both trees demonstrated that the fresh snow density constant [1] (Brazenec, 2005) consistently produced higher R2 (model fit) and lower RMSE (unexplained variance) than empirical variable-density equations [2] (Diamond & Lowry,1953), [3] (LaChapelle, 1962), and [4] (Hedstrom & Pomeroy), in that order. An α-convexity parameter of α = 0.025 was selected for further analysis because, at this level, regressions utilizing the density constant produced slopes close to 1.00 (0.97 for the left tree and 1.07 for the right tree), calibrating the two methods (see Figure 2). Underestimation in TLS based mass estimates is evident with α = 0.010, 0.015, and 0.020 m; overestimation in TLS based mass estimates is evident with α = 0.025 and α = 0.030 m (see slopes in Table 2). Simple linear regression utilizing the density constant and α-shape resolution of 0.025 m yielded R2 = 0.71 / RMSE = 1.06 kg for the left tree and R2 = 0.69 / RMSE = 0.91 kg for the right tree. Intercepts of -1.39 for the left tree and -1.34 for the right tree further illustrate model bias.

The effect of snow density and scan duration on model performance
To our knowledge, this is the first study that explores the suitability of high resolution, automated terrestrial lidar to directly estimate canopy snow interception, with direct comparison to the established hanging tree method (Hedstrom & Pomeroy, 1998). The most precise proxies for measured snow interception mass (R2 [?] 0.69), with the least variation in unexplained variance (RMSE [?] 0.91 kg), were obtained by multiplying the ATLS-derived snow interception volume estimates by a regionally specific fresh snow density constant or 100 kg/m3 (Brazenec, 2005). A density constant of 120-125 kg/m3, a reasonable estimation for the region, would have further improved model performance.
In contrast, ATLS-derived snow interception volume estimates in conjunction with dynamic fresh snow density estimation equations based on air temperature reduced R2 and increased RMSE model estimates. It may be that because equations [2][3][4] are empirical, they represent relationships between fresh snow density and air temperature specific to their respective experimental locales: Sierra Nevada Mountains in California (Diamond & Lowry, 1953); mid-continental Canadian boreal forest (Hedstrom & Pomeroy, 1998); and, a variety of avalanche monitoring sites across the western United States (LaChapelle, 1962). Unexplained model variance in ATLS based mass predictions may also be partially explained by the long ATLS scan duration (13 hours). subsequent changes to the density of freshly intercepted snow (i.e. metamorphosis) during the course of one scan, or retention of metamorphosed snow between snow events, may have resulted in unexplained model variance. Furthermore, this study was not designed to account for losses of intercepted snow due to wind-driven sublimation and/or unloading, processes that also may have resulted in unexplained model variance.
The exclusive reliance on fresh snow density equations that don't take snow metamorphosis into account (Brazenec, 2005;Diamond & Lowry, 1953;Hedstrom & Pomeroy, 1998;LaChapelle 1962) likely contributed to model noise and/or underestimation of predicted density/mass, and experimental development of adjustment factors for these equations was not in the scope of this study. Although under-and overestimation were minimized by conducting a sensitivity analysis to select the degree of α-convexity (see section 4.3), these conclusions point to the need for shorter ATLS scan durations. Recent advances in lidar technology might help to address the slow scan duration time with the relatively new availability of rugged, relatively low-cost (< $10,000), fast scanning lidar instruments (Condliffe, 2018). Further, laser return intensity obtainable from TLS might provide valuable insights on snow properties that affect density, such as changes in grain size/shape and overall wetness (e.g., Eitel et al., 2016, Kaasalainen, Kaartinen, & Kukko, 2008. Although this study emphasized automated data collection, model performance could also have been improved with in situ density measurements of intercepted snow. Further research is of sampling snow density on trees.

The effect of weather and snow properties on model performance
TLS can detect minute changes (<10cm) in snow depth, but several factors can decrease the number and intensity of received signals (Deems, 2013;Prokop, 2008). Atmospheric occlusion from heavy snow or fog can interfere with lidar returns, and wet snow surfaces can lead to adsorption of lidar pulses on the target itself (Deems, 2013;Kaasalainen, Kaartinen, & Kukko, 2008;Prokop, 2008). Snow is also strongly forwardscattering; the proportion of forward scattering of lidar pulses increases with scan angle (Deems, 2013), potentially leading to missed returns on the edges of targets. Despite these issues, the short distance to target (6.2 m) and small scan angle (4.22°) should have minimized the variance between measured and ATLS derived snow mass due to atmospheric occlusion or forward scattering. Likewise, sampling during the coldest months should have minimized unexplained model variance or underestimation of snow volume due to lidar pulse adsorption by snow with high water content. Although these confounding factors were minimized, they were not entirely eliminated in this pilot study.

Conclusion
This study provides valuable insights into the use of TLS for estimating snow interception mass with terrestrial lidar data. Initial results indicate agreement between predicted and measured values of intercepted snow mass (R2 [?] 0.69 and RMSE [?] 0.91 kg) when utilizing an α-shape convexity parameter (α) of 0.025 m and a constant snow density estimate (100 kg/m3). To further improve TLS derived snow interception estimates, future research is needed to develop improved approaches to estimate in situ canopy intercepted snow density, explore the sensitivity of TLS snow volume estimates to changing snow conditions and quantities within the canopies of a variety of live trees of different sizes and for a range of temperatures that affect branch flexibility, and/or reduce beam divergence and reconstruct occluded structural elements. Snow interception is challenging to measure or model, but our findings highlight the potential of lidar technology to efficiently and accurately estimate intercepted snow mass. This is a potentially useful development for the collection of interception data in remote terrain, as well as the calibration of aerial lidar-based snow interception models to distinct forest types around the globe.

Data Availability Statement
The data that support the findings of this study are available from the corresponding author upon reasonable request.