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.