Results

Spatial Patterns of Snow and Forest Structure

Forest structure metrics are highly correlated (Figure A3, Figure 4). LAI’, fVEG, and canopy density tend to increase with elevation and northness. Openness shows the opposite trend, decreasing with elevation and northness (Figure A3). Following expected terrain relationships, absolute SWE is greater at higher elevations on more northern-facing slopes in all three accumulation flights (Figure 5). Lower elevations have steady decreases in SWE in denser canopy, but these relationships become less distinct at higher elevations (Figure A6). The relationships between SWE and terrain are nonlinear and in agreement with previous findings (e.g., Kirchner et al., 2014; Tennant et al., 2015; Varhola et al., 2010) showing that areas with lower SWE tend to be more variable in terms of fVEG and northness, particularly at higher elevations. R-squared values for SWE-elevation regressions are 0.51, 0.68, and 0.43 and slopes are 0.05, 0.19, and 0.09 for the 2008, 2016, and 2022 flights, respectively. R-squared values for SWE-fVEG regressions are 0.12, 0.09, and 0.02 and slopes are -23, -58, and -12 for the 2008, 2016, and 2022 flights, respectively. We calculated that p <0.001 for all models (Figure A9).

“Open Reference Site” Approach

A space-for-structure approach allows us to isolate the effects of vegetation structure by analyzing differences within 30-m grids of up to 900 1-m pixels that we assume experience homogeneous precipitation and the same terrain. We isolate both coarse (fVEG) and fine (LAI’ and openness) forest structures using lidar datasets. This approach is applied across three accumulation flights (2008, 2016, and 2022) and two ablation periods (March to April and April to May, 2016).

Accumulation

DSnowAc patterns from SCB are in general agreement with data used in Varhola et al. (2010). DSnowAc decreases (becomes more negative) with an increase in fVEG, indicating that open areas accumulate more snow relative to average 30-m grid cell values (Figure 7; Figure A7). A linear regression of DSnowAc and fVEG shows greater variability and shallower slopes (-0.23, -0.37, and -0.15 for the 2008, 2016, and 2022 accumulation flights, respectively) than the data found in the Varhola et al. (2010) meta-analysis (-0.37; Table 1).
The RF models have varying levels of skill, with R-squared values of 0.72, 0.62, and 0.16 for the 2008, 2016, and 2022 flights, respectively (Table 2). Forest structure metrics explain a greater amount of variance in the model than terrain metrics (Figure A10). Partial plots of snow accumulation predict DSnowAc values between -4% and -20% per 100% change in fVEG and show that on average the forest sites accumulate less snow than open reference sites (Figure 9). A larger absolute DSnowAc indicates that the discrepancy between open and total grid cell SWE is increasing. In this study, we emphasize patterns consistent across select flights. For the 2008 and 2016 flights, fVEG and LAI’ have the greatest importance on the RF models with about three times as much influence as terrain variables (Figure A10). There are distinct ranges of influence where certain variables predict greater change in DSnowAc. When fVEG < 0.3, we predict consistent decreases of ~5% snow accumulation. At fVEG values > 0.3, these patterns become more variable, and decrease steeply until fVEG ~0.7, where the relationship becomes weaker but continues to predict decreasing DSnowAc of around -15%. Qualitative thresholds can be observed for both LAI’ and openness as well; DSnowAc maintains a maximum around -10% until openness reaches -2.5. A sharp decrease in DSnowAc is predicted until a minimum threshold at around -17% DSnowAc and an openness index of 0 (indicating gap diameter similar to canopy height). DSnowAc shows a steep decreasing trend with LAI’, from DSnowAc around -4% at an LAI’ of 0 to a DSnowAc minimum around -11% where LAI’ is 0.2 (Figure A11). Northness and eastness showed inconsistent results with low importance and no clear trends across the flights (Figure A15; Figure A16). However, the 2008 flight, which best captures accumulation at lower elevations, indicates that southern-facing slopes have a greater difference between open and under-canopy accumulation than other aspects (Figure A14).
Response zones are chosen to quantify snowpack response to canopy removal based on the thresholds identified above (Figure 9). Low response zones are expected to experience the least amount of change post-canopy loss, or post-treatment; high response zones are expected to benefit snow the most. Areas with lower absolute values of DSnowAc (values that approach zero) are typically locations with low density canopy where we would expect greater total accumulation compared with areas with dense canopy and more negative DSnowAc values. Using this framing, we qualitatively identify high (fVEG>0.6 and openness<-2.5), low (fVEG<0.3 and openness>0), and moderate (all remaining fVEG and openness values) snowpack response to canopy loss. Low DSnowAc response areas are primarily at higher elevations, often in the steeper terrain, whereas high response areas tend to be at lower elevations and closer to the valley floor (Figure 6) due to pre-existing patterns in vegetation structure (Figure 4). These response zones show clear differences in average DSnowAc from the space-for-structure approach of ~0, -10, and -30% across low to moderate to high response areas, respectively. Low response areas accumulate more absolute SWE (60 cm in low-response areas vs. 40 cm in high response areas), indicating they would not benefit as much from canopy removal (Figure 10).
The response zones determined through our expanded analysis show clear patterns with previous management and, particularly, forest fires (Figure 11). Areas that have been previously burned are dominated by Ponderosa Pine forest and shrub conversion and have low response classifications. High response areas tend to be forests containing Douglas, Grand, and White Fir. SCB is dominated by moderate response areas, with the proportion of high response areas in planned treatment zones (compared to low-response areas) increasing post-2014.

Ablation

Both DSnowAb metrics show weak increasing linear trends with fVEG and have R2 values of 0.12 and 0.35 for the March-April and April-May flights, respectively (Figure 8). These relationships indicate that an increase in fVEG leads to greater total ablation relative to the open reference sites.
Elevation explains most of the variance in the RF model for the early season ablation (March-April), followed by fVEG and LAI’ (RF r-squared of 0.52; Figure A12). DSnowAb increases as elevation, fVEG, and LAI’ increase. An increase in openness predicts a general decrease in DSnowAb, indicating that sites with more openness show more similarity between ablation in the open and ablation under the canopy. fVEG also controls late season ablation (April-May). fVEG and openness explain most of the variance for late season ablation. High DSnowAb response regions emerge between openness values between -3 (gap diameter 1/20 average tree height) and 0 (gap diameter and tree height are equal) and fVEG values between ~0.5 and 0.9 (Figure A13).