Methods

Study Site

Sagehen Creek Basin (SCB), located in the eastern Sierra Nevada, spans 28 km2 and roughly 1800-2650 meters in elevation. The study domain is about 7 by 8 km, extending beyond the watershed (Figure 1). SCB has a unique record of research spanning more than 50 years . SCB sits within the Sierran Steppe-Mixed Forest-Coniferous Forest-Alpine Meadow Province, and dominant vegetation types include Jeffrey Pine (P. jeffreyi), Lodgepole Pine (Pinus contorta), Ponderosa Pine (P. ponderosa), Red Fir (A. magnifica), Sugar Pine (P. lambertiana), Western White Pine (P. monticola), and White Fir (Abies concolor) . The Mediterranean climate in this region has warm, dry summers and cold, wet winters. Precipitation is dominated by snow . There are three snow telemetry (SNOTEL) sites within or directly surrounding the basin (http://www.wcc.nrcs.usda.gov/snow/) which measure depth and SWE using snow pillows. Median peak SWE ranges from 37 cm at the lowest SNOTEL site (Station name: Independence Creek; Station ID: 540, elevation 1961 m) to 113 cm at the highest site (Station name: Independence Lake, Station ID: 541, elevation 2541 m). The mean 30-year winter (DJF) temperature ranges from -1.7°C to -1.8°C from the lowest to highest SNOTEL site (Figure 1).
The ecohydrology of SCB since ~1200 CE has been shaped by anthropogenic influences. Indigenous forest management (e.g., prescribed burns), primarily influenced by the presence of the Washoe (Wá∙šiw) Tribe, formed a heterogeneous ecosystem from ~1200 CE until the mid-19th century with fire return intervals around two years . Subsequent settlement led to a period of timber harvesting until the early 20thcentury and fire suppression thereafter. This suppression resulted in dense, homogenous forests susceptible to high-intensity disturbance (e.g. wildfire, drought; . The “Sagehen Fuels Reduction Project” (Sagehen Project), a collaborative effort aimed at restoring and enhancing ecological heterogeneity and resilience to wildfires, has had a substantial impact on the landscape of SCB. Fuel treatments, including prescribed burns and forest thinning, were implemented at target sites in the late 2010s (Figure 1c). The goal of these management initiatives was to restore ecosystem function and enhance species conservation. However, less is known about how these treatments, or landscape disturbances, explicitly influence hydrological processes driven by snowmelt .

Lidar Data and Processing

Lidar is an active remote sensing tool that emits pulses of light at near infrared wavelengths (typically 1064 nm). At the target scales for this study (individual tree to hillslope), lidar data is collected using instruments aboard fixed-wing aircrafts. Full waveform pulses are emitted from lidar scanner and reflected off objects and the stronger returns are discretized into points based on the return time of the pulse. Lidar resolution is quantified as a function of ground and non-ground (e.g. vegetation) returns per square-meter (See Table A4 for resolution specifications of the data used in this study).
Major limitations of airborne lidar data are the temporal resolution and cost of acquisition, which can inhibit multiple flights. While peak SWE (historically around April 1) is typically targeted for lidar data acquisitions, SNOTEL data often reveal ablation at lower elevation sites before that date . Therefore, early lidar flights can be valuable to investigate accumulation dynamics. However, climate variability complicates the prediction of peak SWE timing, and later flights can provide more useful information to stakeholders interested in the magnitude of spring melt.
The study period covers 2008 to 2022 (Table A4), with specific dates of inquiry driven by lidar data availability. Later flights were impacted by significant ablation, particularly at lower elevations (Figure 2). Despite their temporal limitations, these data present a unique opportunity because a series of multi-temporal flights in March, April, and May 2016 captured both early and late season ablation, allowing us to take advantage of spatially distributed SWE data. Lidar data was sourced from ASO and the National Center for Airborne Laser Mapping (NCALM). Snow depth was calculated by differencing snow-on and snow-off flights. Snow-off lidar was flown by NCALM in the summer of 2014 with an average point density of 9 pts/m2 . Snow-on lidar was flown in the early accumulation season of 2008 (on February 10) by NCALM, with an average point density 5 pts/m2 and again in 2016 by ASO at near peak SWE (March 26), early in the ablation season (April 17), and late in the ablation season (May 18) with average point densities ~2 pts/m2 . A second snow-off lidar was acquired by NCALM in the fall of 2020 (11/20) with an average point density of 21 pts/m2 . Snow-on lidar post-2020 was also flown by NCALM. This flight was targeted to capture peak SWE in the spring of 2022. Early season accumulation followed by a long period of ablation shifted the peak SWE ~10 days earlier than April 1. Therefore, the data was collected on March 21 to target peak SWE with an average point density of 28 pts/m2 .
Data processing was modified from the workflow of Kostadinov et al., 2019 (See Appendix A for more details on lidar processing) to optimize reproducibility and to enable analyses in diverse vegetation types. Vertical biases were corrected across the watershed based on available ground-truth points (Table A4; Figure A4).
After the vertical bias was corrected, we classified the snow-off vegetation using a refined classification logic (Table A1) to produce four canopy strata classes: open, short, understory, and tall (Figure 3). In addition, an open reference class represents open areas with a 1 m buffer around the canopy to mitigate edge effects for our analyses (Figure A5). Snow depth was calculated by taking the height above ground and filtering based on the canopy strata classifications (Table A2). Snow density was modeled using the Snow Physics and Lidar Mapping (SnowPALM) model, which calculates compaction as a function of albedo and snow density decay using values calibrated from SNOTEL sites, as well as the effect of canopy interception and sublimation of snowfall, shortwave transmittance through the canopy, and longwave radiation emitted by the forest, and the full energy balance of the snowpack . Where data was available, SnowPALM was validated against field and SNOTEL data . The median difference between SnowPALM and ground-truth snow density was <1 kg/m3. We did not correct for this discrepancy in our final SWE calculations, since this difference is smaller than the range of spatial variability in density. The SnowPALM model calculated snow density for each water year at 1-m spatial resolution. From lidar snow depth and SnowPALM density, we calculated SWE at 1 m (pixel) and averaged that to 30 m (grid) resolution.
We chose several key terrain variables to characterize the landscape as proxies for the energy and mass-balance dynamics that control snowpack accumulation and ablation . Northness (cos(aspect) * sin(slope)) primarily represents solar radiation while eastness (sin(aspect) * sin(slope)) represents wind loading that alters mass redistribution, since the predominant wind direction is from the west. Elevation correlates with both precipitation and temperature.
Forest structure metrics were chosen to represent the wide range of classifications that are used in other studies consistent with coarse metrics (i.e. Varhola et al. 2010) and newer, lidar-based metrics of fine forest structure . See Table A3 for forest structure classification descriptions. We chose to investigate the fraction of vegetation (fVEG), a modified leaf area index (LAI’), and openness index. fVEG is a coarser-scale forest structure metric that indicates the portion of a grid cell covered by pixels with maximum canopy height greater than three meters, impacting both the incoming radiation dynamics and interception . LAI is the most consistent with previous studies used by Varhola et al. (2010) that used full “forest cover” to define dense, interacting canopies. LAI’, calculated using a combination of a point density ratio and canopy height model, serves as a proxy index for the total surface area of the leaves or needles on the canopy, acting as a limiting factor for canopy interception. This metric was then scaled to the maximum forest height . Openness represents gap dynamics influencing incoming solar radiation and breaks in potential interception. We defined openness as the ratio of the gap diameter (or twice the maximum distance to nearest canopy) to the average forest height in each grid cell. We took the natural log of the openness variable so that negative values represent forest gap to average tree height ratios less than one and positive values represent forest gap to average tree height ratios greater than one .

“Open Reference Site” Analysis

The space-for-structure approach aims to control for terrain and isolate secondary vegetation impacts on snow. The primary assumption associated with this method is that the dynamics dictated by the forest structure at one location are representative of dynamics that would occur if the structure at another location shifted towards the former.
Following the methodology of Varhola et al. (2010), SWE was normalized based on a “reference” value using Equation 1, below. This normalization was calculated to control for precipitation and terrain differences across the study domain.
\(\Delta y\,=\,100\frac{x-R}{R}\) (1)
Here, ∆y represents the percent difference between the target value (x) and a reference value (R). These percent change values were calculated for several inputs. The general form of Equation 1 was used to determine how accumulation changed with increasing canopy cover. Open reference areas were defined as areas 1 m from the canopy edge with no canopy (Figure A5).
For the lidar flights during the accumulation season, the average SWE of all 1-m pixels within a 30-m grid cell was taken to represent the target SWE (x), and the average SWE of the open reference 1-m pixels within that same 30-m grid cell was taken to represent the reference SWE (R). The details of the open reference approach are explained in Appendix A and match the Varhola et al. (2010) framework. Low elevation areas (<2,100 m) were excluded from the 2016 and 2022 analyses due to ablation; ablation removes the accumulation signal. We refer to the final metric from Equation 1 as delta accumulation (DSnowAc).
To calculate ablation, we took the SWE loss between each ablation season flight to produce early season (March-April) and late season (April-May) ablation metrics at 1-m resolution. We limited early season ablation data to lower elevations (<2,100 masl) and late season ablation data to higher elevations (>2,100 masl). Because of the limited temporal resolution of these data, if a pixel had no SWE, we could not be sure when that pixel ablated completely. Thus, we excluded all snow-free data, in addition to any areas that experienced accumulation.
For the ablation flights, we modified Equation 1 to create a normalized difference ablation metric (delta ablation, or DSnowAb).
\(y=\ 100\cdot\left(\,\frac{x_{i}-x_{f}\,}{x_{i}}-\,\frac{R_{i}-R_{f}}{R_{i}}\right)\ \)(2)
Where “i” and “f” identify the SWE from the initial and final lidar flights, respectively, within the ablation period.

Random Forest Analysis

We used the Random Forest (RF) machine learning tool implemented through the Ranger package in R to calculate the controls on snowpack accumulation and ablation. RF classifiers use a set of regression trees created using bagged samples from a training dataset. We chose to use a RF because it is an ensemble classifier that can be used with multi-modal data and accounts for covariance in the input variables, and it has been previously used to understand the impact of canopy cover on snowpacks (Krogh et al., 2020, Lewis et al., 2023).
We experimented with several forest structure metrics based on height and lidar point density, and narrowed target metrics based on both preliminary RF runs and knowledge about specific density metrics that have been found to be important in this region . See Table A3 and Figure A3 for details about tested metrics.
We performed RF analyses with ∆y determined from Equation 1 and Equation 2 as the response variables (See Appendix A for more details on RF hyperparameters). Predictor variables for the space-for-structure approach remained the same: elevation, northness, eastness, fVEG, openness, and LAI’.

Decision Support Tool

Based on the results of the RF models, we identified a series of “response” thresholds from slope changes in plots of select predictor variables. Thresholds represent areas where thinning impacts snowpack response (specifically DSnowAc response) either positively (more accumulation) or negatively (less accumulation). We incorporated the results into a simple decision support tool for stakeholders to make informed decisions about forest management. The goal was to capture forest structure thresholds that could be used to target forest structures for thinning. Using the Sagehen Project as a natural experiment, we tested these thresholds against historically disturbed areas to evaluate the effectiveness of thinning on snowpack resilience. We expanded response criteria classification to a larger area covering the eastern extent of the 2014 NCALM flight with similar existing vegetation types (EVT) to SCB (from https://www.landfire.gov/evt.php) to demonstrate how this method could be scaled. The study domain was dominated by coniferous vegetation, aggregated based on the collapsed vegetation type name and physiognomic classification. This allows us to qualitatively assess how broad-scale patterns align with finer-scale patterns and investigate how prior management and disturbance correlates with response classification.