Discussion

Comparison with the Varhola et al. (2010) Framework

Maximum snow accumulation has a negative relationship with increasing fVEG using our “open reference site” lidar-based approach, which aligns with the meta-analysis of Varhola et al. (2010). Our approach takes advantage of the lidar point density to compare large scale (30-m) snow distributions with a nearby open reference location (Figure A5), which controls for factors such as precipitation and terrain. However, a limitation of this approach is that the “open reference” area implicitly depends on the amount and arrangement of vegetation (e.g. not all open 1-m areas are the same). Using the flights that best approximate maximum snow accumulation across three years (i.e. the least ablation at lower elevations), we find a linear trend in relative snow accumulation between fully open and 100% fVEG, with differences of 23% in 2008 and 37% in 2016 (R2 of 0.72 and 0.62, respectively). These slopes are either similar or smaller than those from Varhola et al. (2010), who calculated a slope of 0.37. We see stronger relationships in February 2008 but a lower slope (Figure 7), potentially due to the lack of ablation at lower elevations as compared with March 2016. More lower elevation ablation may accentuate DSnowAc differences because under canopy areas are expected to ablate earlier than open areas in this warmer environment . The correspondence between 2008 and 2016 patterns in DSnowAc with fVEG support the “open reference site” approach and the resulting process inferences we describe below. A random forest model demonstrates that fVEG accounts for the most variation of any predictor variable (fVEG, openness, LAI’, northness, elevation) in two of the accumulation models. The RF results show weaker effects when fVEG is below ~30%. The shallower slopes in the RF partial plots (Figure 9; ~10% change in DSnowAc for 100% change in fVEG) compared with the linear relationships suggest the importance of factors besides fVEG (discussed more in Section 4.2).
Our “open reference site” approach predicts greater ablation with increasing fVEG, and thus reductions in ablation rates with canopy removal, which is the opposite of the relationships found by Varhola et al. (2010). However, our results are consistent with process-based studies in warmer climates, where high levels of longwave radiation from trees may overcome the decreased solar radiation in dense forests, especially on shaded north facing slopes Our results demonstrate a need to refine the Varhola et al. (2010) model of ablation for different climates. Steep north-facing slopes may not receive as much increase in solar radiation after thinning compared to flat areas, which would strengthen the signal of higher ablation with higher fVEG in these areas. Greater under-canopy ablation occurs from a combination of the timing of the ablation season (March-May, when solar irradiance is relatively low) and trees that emit enough longwave radiation to dominate energy fluxes . We find an average of ~15% change in ablation for 100% change in fVEG for the March-April and April-May periods (R2 of 0.13 and 0.34, respectively in Figure 8). Partial plots in the RF model also show positive relationships with fVEG and ablation rate but indicate that there is little effect when fVEG is below ~30% (Figure A13). Ablation rates calculated between two lidar flights are different than season-long ablation metrics used by Varhola et al. (2010). In particular, the lidar-based ablation metrics are affected by snow accumulation during the ablation season, variation in precipitation amount, cold content, energy inputs, and other factors that are highly spatially variable (Figure A3). We attempt to mitigate these issues by using early season (lower elevation) and late season (higher elevation) ablation windows and controlling for initial SWE in each case. However, opportunities exist to improve our method to better utilize ablation information by including distributed snow disappearance estimates with modeling or optical remote sensing. An increasing number of open-source, multi-temporal lidar datasets, collected by ASO for example, in recent years presents an opportunity to repeat these methods across the Sierra Nevada and into domains, like Colorado and Utah, that may show different results due to a variety of energy balance dynamics. Importantly, to be the most useful for this method, these flights must be paired with timely, high point-density, snow-off lidar as well.

Role of Finer-Scale Forest Structure and Terrain on Snow Dynamics

The predictor importance shows that fVEG, LAI’, and openness are more important than terrain metrics (elevation, northness, and eastness) for predicting our DSnowAc metric in the RF model (Figure A9). The importance of these vegetation metrics was clearly shown by our “open reference site” approach since analyses of raw snow accumulation data would fail to disentangle the role of vegetation from the strong influence of elevation and snowfall spatial patterns (Figure 4 and Figure 5). The partial plots from the RF model consistently show a non-linear effect on accumulation from LAI’ and openness, as compared to a more linear relation with fVEG. Overall, our results imply that vegetation has the largest role in reducing snow accumulation when fVEG is high (>0.6), openness is low (vegetation height ≥ canopy gap diameter), and LAI’ is high (>0.25). Inversely, these results indicate that a canopy gap diameter to vegetation height ratio <1 optimizes accumulation, in general support of earlier studies that demonstrated maximum shading across the gap when nearby canopy height is comparable to the gap diameter (i.e. openness ratios <1) . Conversely, partial plots for terrain effects on accumulation and ablation have shallower slopes, with mixed slopes across flights suggesting weaker controls (FigureA14-A18). Our results match Krogh et al. (2020) and Lewis et al., (2023), who used the SnowPALM model to conclude that lower-elevation, south-facing slopes experience post-thinning increases in accumulation.
We also show the importance of fine-scale forest structure controlling ablation rates, and promoting snow retention following thinning, despite limitations in our random forest analysis. Because of differences in energy environments, elevation is important for predicting DSnowAb, whereas fVEG, openness, and LAI’ are more important for the early ablation season metric. Both early and late season ablation estimates indicate a combined impact of fVEG (lowest ablation at <~50% canopy cover) and openness (lowest ablation where gap diameter is greater than tree height) to promote retention. The similar response of snow accumulation and ablation suggests that forest thinning in dense areas – especially if the thinning creates large gaps rather than uniform density reduction – is likely to slow ablation while increasing accumulation, thereby increasing overall snow retention.
Despite the importance of a coarse (simpler) fVEG metric based only on tree spatial coverage, we find that metrics accounting for vertical lidar point density (LAI’) and horizontal and vertical canopy information (openness) help refine predictions of snow-forest interactions. For example, different species of trees often have different LAI’ and height, while in forest patches with similar forest species the spacing and orientation of the trees may be different. Similarly, the land use history or resource limitations of the local area could be important for coarse and fine-scale forest structure information (Stephens et al., 2021). It’s worth noting that these forest structure metrics are highly correlated across the full domain, but their correlation within elevation bands is more mixed due to more homogeneity in forest species, terrain, land use history, and other factors (Figure 4). fVEG is the most important vegetation metric in our analysis, compared with openness and LAI’, is easier to estimate, and matches the linear relationship found in Varhola et al. (2010). Importantly, however, Varhola et al. (2010) did not explore finer scale forest structure information like openness and LAI’. Varhola et al. (2010) also emphasized the inconsistencies between forest structure metrics across studies, which hints at the value of using lidar information to develop repeatable coarse and fine forest structure information.

Management applications and future research directions

While simple linear models capture broad patterns of forest-snow interactions, thresholds at which specific forest structure metrics become important can help refine forest management decisions beyond the current decision-making paradigm that typically does not include benefits to snow water resources. We developed a decision support tool for SCB, and extrapolated that to a landscape scale, by defining qualitative thresholds in fVEG and openness that correspond to high, moderate, and low response of SWE to potential thinning. The Sagehen Project was a mixture of mechanical thinning and prescribed fire that was planned to address wildlife and fire concerns and resulted in the fragmented treatment areas shown in Figure 1. We find near zero DSnowAc (~0%) and more absolute SWE (~50 cm) in the low response areas of SCB, suggesting that our threshold method has identified areas with little response to thinning or other canopy removal. These areas account for a small portion (~7% in 2014) of the total domain, suggesting medium and high response areas were widespread across SCB prior to thinning from the Sagehen Project. The high response areas had a delta SWE of around -25% and lower absolute SWE of ~30 cm. Only 29% of the treated areas as part of the Sagehen Project thinning were high response areas (pre-thinning), while 27% were moderate response. While the resilience of snowmelt-derived water resources was not an explicit priority for the Sagehen Project, it illustrates the potential co-benefits of simultaneously managing forest fire reduction, wildlife habitat, and snow retention . Our decision support tool offers both a retrospective means to assess restoration effectiveness (as done in the Sagehen Project area), as well as a tool for proactive planning efforts.
One of the key advantages of our results for decision support purposes is that it can be extrapolated to larger domains of similar climate and forest conditions or re-developed in new areas with different snow conditions (e.g. colder and windier) that have existing snow-on and snow-off lidar datasets. To demonstrate the scalability of our results for decision support, we extended our low, medium, and high response mapping to an 800 km2 landscape-scale area with planned forest restoration by the USFS. Lidar analyses showed that treated areas had a higher relative proportion of low response zones compared to areas that had not yet been treated, suggesting that treatment planning was optimizing areas for increased snow retention. The increase in the treatments in high response zones after 2014 suggests that this type of lidar-based analysis could have been used proactively to plan treatments for maximum snow effects. Areas near SCB that have burned in recent years (Figure 11b), especially those that have converted to shrub cover, show a much higher proportion of low response areas compared to surrounding unburned and un-thinned forestland. This difference suggests that the fires were effective at changing forest structure in a manner that would increase snow accumulation and retention, but that land cover conversion from forest to shrub may have impacts not estimated with our decision support tool. While wildfires may increase ablation due to black carbon reducing albedo several years post-fire , other studies have shown that reductions in canopy cover from fire can provide net benefits to snowpack and streamflow in the Sierra Nevada . Our method does not consider the effects of fire or other processes that disturb the snow albedo or energy budgets beyond physical changes in forest structure. The simplicity of our methods suggests that replicating our method in different forest types, different forest restoration treatments, and across a range of climates with lidar data coverage could yield both process insight and new decision support tools.