The prediction of uneven snowpack response to forest thinning informs forest restoration in the central Sierra Nevada

The Sierra Nevada has experienced unprecedented wildfires and reduced snowmelt runoff in recent decades, due partially to anthropogenic climate change and over a century of fire suppression. To address these challenges, public land agencies are planning forest restoration treatments, which have the potential to both increase water availability and reduce the likelihood of uncontrollable wildfires. However, the impact of forest restoration on snowpack is site specific and not well understood across gradients of climate and topography. To improve our understanding of how forest restoration might impact snowpack across diverse conditions in the central Sierra Nevada, we run the high‐resolution (1 m) energy and mass balance Snow Physics and Lidar Mapping (SnowPALM) model across five 23–75 km2 subdomains in the region where forest thinning is planned or recently completed. We conduct two virtual thinning experiments by removing all trees shorter than 10 or 20 m tall and rerunning SnowPALM to calculate the change in meltwater input. Our results indicate heterogeneous responses to thinning due to differences in climate and wind across our five central Sierra Nevada subdomains. We also predict the largest increases in snow retention when thinning forests with tall (7–20 m) and dense (40–70% canopy cover) trees, highlighting the importance of pre‐thinning vegetation structure. We develop a decision support tool using a random forests model to determine which regions would most benefit from thinning. In many locations, we expect major forest restoration to increase snow accumulation, while other areas with short and sparse canopies, as well as sunny and windy climates, are more likely to see decreased snowpack following thinning. Our decision support tool provides stand‐scale (30 m) information to land managers across the central Sierra Nevada region to best take advantage of climate and existing forest structure to obtain the greatest snowpack benefits from forest restoration.


| INTRODUCTION
California's Sierra Nevada snowpack is the primary water resource for agricultural and urban use, with economic importance in the tens of billions of dollars annually (Stewart et al., 2004;Zheng et al., 2018).
The future of Sierra Nevada snowpack, however, is threatened by anthropogenic climate change (Riddlebarger et al., 2021) and wildfire (Smoot & Gleason, 2021), as decreasing snowpacks are melting earlier in the spring across much of western North America (e.g., Knowles et al., 2006;Mote et al., 2005).The Sierra Nevada currently provides over 60% of California's developed water supply, much of which has historically fallen as snow during the colder months (Edelson & Hertslet, 2019).As average winter and spring temperatures have increased across western North America, spring snowmelt-driven streamflow has shifted earlier and decreased in magnitude, bolstering the importance of conserving existing snowpacks to provide ample water supply for tens of millions of end users and the fourth largest economy in the world (Mote et al., 2005).
Early studies in the Sierra Nevada have shown the significant differences between snow accumulation under the canopy and forest clearings (Church, 1933), endorsing the possibility of increasing snow retention through proper forest management.One technique to increase snow accumulation and melt is through selective forest thinning, which can minimize the primary drivers of snow loss, specifically canopy-intercepted snowfall, longwave radiation, and wind (Tague et al., 2019;Varhola et al., 2010).Previous studies in cold regions have demonstrated that up to 50% of cumulative snowfall can be intercepted by, and up to 60% of snowfall can be sublimated from, the canopy during midwinter conditions (Hedstrom & Pomeroy, 1998).Likewise, canopies reduce wind speed thereby lowering turbulent sensible and latent heat exchange (Conway et al., 2018) and incoming direct shortwave radiation from shading of the snow surface (Musselman et al., 2013(Musselman et al., , 2015)).Wind controls snow redistribution and enhanced wintertime sublimation, while also increasing snow ablation during the melt season, when sensible energy warms the snowpack (Varhola et al., 2010).Prior research highlights the significant regional variability and complicated interactions between snowpack response and vegetation removal from differences in climate, topography, forest type, and the spatial distribution of trees on the landscape (Varhola et al., 2010).Since trees have a limited snow loading capacity, the influence of forest cover decreases in importance as the total snowfall increases beyond a certain threshold (Boon, 2009).
Vegetation disturbance, whether natural or human-caused management actions, can alter both the timing and quantity of snowmeltdriven streamflow.There is, however, substantial variation in both the magnitude and direction of hydrologic response depending on physiographic features and forest structural conditions (Goeking & Tarboton, 2022;Ren et al., 2021).For example, snow storage duration in forest stands relative to adjacent open areas has been observed to range from lasting more than 3 weeks longer (Koivusalo & Kokkonen, 2002;Lundquist et al., 2013) to disappearing more than 12 weeks earlier (Dickerson-Lange et al., 2017).Consequently, predictions for the net effects of forest cover change on snow storage from management, such as forest thinning or wildfire, are typically site specific (Gleason et al., 2013;Harpold et al., 2020;Schneider et al., 2019), reflecting differences in climate and vegetation structure, with limited ability to make inferences about unstudied areas (Dickerson-Lange et al., 2021).
Thinning and prescribed fire are being implemented at unprecedented levels across the western United States to reduce wildfire hazard, increase forest resilience, and improve forest health.An example of these landscape-scale forest restoration efforts is the 6300 km2 Tahoe-Central Sierra Initiative (TCSI; see Figure 1; Manley et al., 2020Manley et al., , 2022) ) that spans across the American, Truckee, and Yuba River watersheds, as well as the Lake Tahoe basin.The TCSI is a multiagency coalition aimed to accelerate landscape forest restoration to "improve the health and resilience of the Sierra Nevada" (Manley et al., 2022).Recent snow dynamics modelling along the west shore of Lake Tahoe showed that reducing the density of shorter trees (<10 and <20 m tall) resulted in the most benefit to the snowpack in terms of maximum depth and timing of melt (Harpold et al., 2020;Krogh et al., 2020).However, these studies focussed only on the higher elevations that were characterized by dense and fire-suppressed forests and may not represent the much larger and more climatically and Remote sensing-based tools offer the potential to account for local forest structure, which has successfully been used to estimate post-thinning evapotranspiration (ET) response in the Sierra Nevada using multispectral imagery (Roche et al., 2018) and three-dimensional forest structure derived from lidar (Harpold et al., 2020).However, broadly available decision support tools to enhance the hydrologic outcomes of management do not yet incorporate criteria based on forest structure and operational constraints associated with different tree harvest methods.For example, heavy machinery can harvest a wide range of tree diameters but cannot operate on slopes steeper than 35% or in areas over 300 m from established roads.Conversely, hand thinning can be conducted in any terrain but cannot remove vegetation over $35 cm diameter at breast height (roughly corresponding to a 10 to 20 m tall tree; Crowfoot & Crowell, 2021;Tahoe Conservancy, 2021).Additionally, forest restoration projects must account for threatened forest-dependent species habitats, such as the California spotted owl and northern goshawk (Cabiyo et al., 2021;Crowfoot & Crowell, 2021).Silvicultural prescriptions commonly provide target basal area, canopy cover, and maximum tree diameters for removal, but the selection of individual trees for removal is commonly made in the field at the stand scale, and it is rarely informed by landscape-scale ecosystem benefits, such as water conservation and yield.A combination of stand-scale information with fine-scale, physically based snow modelling into decision support tools would greatly enhance the ability of managers to design silvicultural prescriptions to enhance water-related benefits.
Currently, there are only a handful of tree-scale snow energy and mass balance models capable of simulating the effect of heterogenous silvicultural prescriptions on snow water availability, such as the Flexible Snow Model (FSM2; Mazzotti et al., 2020)

and Snow
Physics and Lidar Mapping model (SnowPALM; Broxton et al., 2015), which have historically had limited translation into management support tools (Krogh et al., 2020;Mazzotti et al., 2021).These models must be run over larger and more heterogeneous forests and climates to best inform management solutions and are important for simulated response to forest restoration because stand-scale snowpack response requires accounting for effects of local radiation at the tree to gap scales (Broxton et al., 2021;Mazzotti et al., 2023).
Here, we use the high-resolution (1 m) SnowPALM model on five unique 23-75 km 2 subdomains to include a wide range of topography, climate, and vegetation.SnowPALM was developed to quantify the effects of vegetation structure on stand-scale snow accumulation and ablation in New Mexico and Colorado.One important advantage of SnowPALM is that it parameterizes the 1 m vegetation height and vertical canopy density based on airborne lidar, which is available at large spatial scales.This model has since been used to calculate the effect of forest restoration on snow retention and water yield in the Sierra Nevada (Harpold et al., 2020;Krogh et al., 2020), as well as to quantify the effect of forest structure change from a high-severity fire on snowpack dynamics in New Mexico (Moeser et al., 2020)
Each of these subdomains occur within the broader TCSI landscape.
The TCSI spans a wide elevation range from 1500 to 3311 m.a.s.l, a Mediterranean climate, and mean annual precipitation of $1500-2000 mm, roughly 60% to 90% of which occurs as snowfall, depending on elevation.The forests across the TCSI primarily consist of Sierran mixed conifer, including Jeffrey pine (Pinus jeffreyi), ponderosa pine (Pinus ponderosa), lodgepole pine (Pinus contorta), Douglas fir (Pseudotsuga menziesii), and white fir (Abies concolor), along with upper elevation subalpine conifer, commonly dominated by red fir (Abies magnifica), with scattered stands of aspen (Populus tremuloides) and variable understory species (Krogh et al., 2020).We selected the five subdomains to capture a wide range of meteorological, topographic, and vegetation characteristics to better characterize the variability of the central Sierra Nevada and to include locations where forest restoration is either planned or has been recently completed.
All of the planned restoration and fuels reduction projects aim to address critical barriers to increasing the pace and scale of forest restoration in the Sierra Nevada.For example, the French Meadows project proposes to treat 4800 ha with mechanical thinning, mastication, hand thinning, reforestation, and aspen and meadow restoration, as well as 2800 ha of prescribed burning, in collaboration with numerous state and federal agencies (Edelson & Hertslet, 2019) Lastly, the Lake Tahoe West Restoration Partnership (LTW) is working collaboratively to restore the resilience of 24,000 ha of forests, watersheds, and communities across the west shore of the Lake Tahoe basin.One of the stated goals of the LTW project is to use thinning to help forests recover from fire, drought, and insect outbreaks, as well as to ensure that fires burn primarily at low to moderate severity for ecological benefits (Tahoe Conservancy, 2021).

| Snow modelling
The Snow Physics and Lidar Mapping (SnowPALM) model (Broxton et al., 2015), used to simulate the impact of forest thinning on snowpacks in this study, calculates snow mass and energy balance at a high spatial (1 m) and temporal (hourly) resolution, enabling detailed snowpack simulation in areas with complex terrain and heterogeneous vegetation.Input data for the model include elevation and vegetation (canopy height and density) layers derived from lidar, as well as a standard suite of meteorological data (precipitation, air temperature, pressure, wind speed and direction, and incoming shortwave and longwave radiation).SnowPALM simulates the mass and energy balance of the snowpack and soil layers separately, in addition to calculating the energy exchange with the atmosphere.The model accounts for attenuation and enhancement of shortwave and longwave radiation based on the location of tree shadows, sky view factor, and parameterized longwave enhancement from tree trunks (Mahat & Tarboton, 2012), interception and sublimation (Hedstrom & Pomeroy, 1998;Liston & Elder, 2006), and wind distribution of snowfall (Winstral & Marks, 2002).The model simulates the impact of precipitation events, including rain-on-snow, as incoming snowfall adds directly to the frozen water in the snowpack, but rainfall can either freeze to the snowpack (if the snowpack is below freezing) or percolate through the snow into the soil (if the snowpack is isothermal).We use the model to simulate the total water emanating from the snowpack or atmosphere that is available for infiltration and runoff, which we will refer to as RAIM, or the sum of snowmelt (including percolated rain-on-snow) and rainfall on bare ground.Further details about SnowPALM can be found in Broxton et al. (2015).
Following the methods of Krogh et al. (2020), we use hourly precipitation, air temperature, wind speed and direction, air pressure, downward shortwave and longwave radiation, and specific humidity at 1/8 degree resolution from the North American Land Data Assimilation System (NLDAS-2; Xia et al., 2012).These data are downscaled to 1 m resolution using lapse rates (for temperature and precipitation) from the Parameter-elevation Regressions on Independent Slopes Model (PRISM) monthly climate data, as in Krogh et al. (2020) and Moeser et al. (2020).To speed up computation time, SnowPALM uses parallel processing for 500 m tiles within the subdomain.We used a high-performance computer with 32 cores at the University of Nevada, Reno, to execute the model using vegetation height and canopy density generated for the two thinning scenarios and for the current forest conditions.We ran SnowPALM from 1 October to 31 July for water year 2015-2016 (hereafter WY16), 2016-2017 (WY17), and 2017-2018 (WY18) for each subdomain.These three years spanned a range of annual meteorological conditions from relatively dry (WY18; 84% of normal precipitation) to average (WY16; 98% of normal) to relatively wet (WY17; 176% of normal).
The model parameterization used in this study is largely the same as Harpold et al. (2020) and Krogh et al. (2020).However, we finetuned several SnowPALM variables using snow depth and snowwater equivalent (SWE) measurements from Snow Telemetry (SNOTEL) and other weather stations within our study areas: East Shore (Marlette Lake, #615), West Shore (Ward #3, #848), Yuba (Station ID: GOL), Sagehen (weather stations within SCEF; Kirchner et al., 2020), and French Meadows (Dolly Rice stations and Talbot Camp data from Bales et al., 2020).We removed biases for air temperature and precipitation, adjusted albedo parameters to match melt rates from peak SWE to snow disappearance, and used maximum leaf area index (LAI) from local in situ measurements (van Gunst, 2012).

| Lidar datasets and virtual thinning experiments
We used a suite of airborne lidar collections to calculate the vegetation height, canopy cover, and digital elevation model for our subdomains, including data from the National Center for Airborne Laser Mapping's (NCALM) 2013-2014 USFS Tahoe National Forest campaign (for French Meadows, Sagehen, and Yuba; https://doi.org/10.5069/G9V122Q1) and Watershed Science's 2010 Lake Tahoe Basin campaign (East and West shores; https://doi.org/10.5069/G9PN93H2), as well as the 2018 USGS 3DEP campaign (for TCSI; USGS, 2015).
We used the FUSION/LDV LIDAR analysis software (McGaughey, 2014) to calculate the vegetation height, which we define as the difference between the canopy height model and the bare-earth elevation at 1 m resolution.We also used the FUSION software to calculate the percent tree canopy cover, defined as the ratio between tree cover (>2 m above the ground surface) and total returns per square metre.
Published tree segmentation data from Xu et al. (2018) was applied for the West Shore and East Shore, while we calculated our own tree segmentation data for French Meadows, Sagehen, and Yuba using similar methods.The Xu et al. (2018) algorithm uses the lidar point cloud to delineate each tree and create shapefiles outlining the approximate footprint and maximum height of individual trees.Using the footprint and maximum tree height, we create the two virtual thinning scenarios for the SnowPALM modelling by removing the vegetation when the maximum tree height was less than 10 or 20 m tall.Specifically, we set the corresponding vegetation height and canopy cover to zero for all pixels within the footprint of each tree we remove.These two height thresholds mimic light and heavy forest thinning prescriptions from the USFS, respectively, across the Sierra Nevada designed to improve drought tolerance and reduce forest fuels.These treatments aim to lower the risk of severe forest fires, while maintaining other ecosystem benefits, such as water conservation and yields, wildlife habitat, and recreational values.

| Random forests decision support tool
We developed random forests (RF; Breiman, 2001)  Next, we used the results from the SnowPALM modelling within the five subdomains to create a decision support tool to inform land managers on how to selectively thin vegetation to maximize RAIM both within and outside of our five subdomains.We used an RF algorithm, this time trained on all the data from the five SnowPALM subdomains, to predict the most productive locations to thin within the TCSI region.The RF model was developed to predict how changes in vegetation height and canopy cover affect snow accumulation and ablation under a wide range of topographic (elevation, slope, northness, and eastness) and climatic (air temperature, precipitation, and wind) conditions, as well as the initial vegetation height.
An "out-of-bag" importance metric (Breiman, 2001) for each predictor was calculated by the random forests and compared with the original RF out-of-bag dataset, following the methodology from Krogh et al. (2020).Predictor variables with higher errors after the random permutation were the most sensitive to change and thus are the most important in this analysis.

| SnowPALM model results
After calibrating our SnowPALM models to SNOTEL sites located within each subdomain, the resulting root-mean-square error (RMSE) for WY16-18 SWE for the West Shore is 106 mm, Sagehen is 209 mm, Yuba is 210 mm, and East Shore is 72.8 mm, while the WY15-17 snow depth RMSE for French Meadows is 133.8 mm (see Figure 2).
We successfully ran the SnowPALM model over each of the five subdomains for all thinning scenarios for three water years, producing daily outputs for all desired variables.An example histogram of the vegetation height from the Sagehen subdomain before and after both light (removing vegetation <10 m tall) and heavy (removing vegetation <20 m tall) thinning can be seen in Figure 3a.To give a better idea of where the USFS could enact vegetation removal, the shaded region in Figure 3b shows a 300 m buffer from all established roads across terrain with slopes less than 35%, superimposed upon the vegetation height in Sagehen.Over 75% of the Sagehen subdomain, and nearly 65% of the TCSI, lie within 300 m of an established road.
Histograms of ΔRAIM show that the magnitude of increased or decreased RAIM varied between the control and the two restoration scenarios, light thinning (Figure 4a) and heavy thinning (Figure 4b).RAIM increased across an average of 96% and 92% of the West Shore and Sagehen subdomains, respectively, for all three water years and both thinning scenarios (Figure 4).However, RAIM only increased across 36% and 28% of the Yuba and East Shore subdomains, respectively.RAIM increased across 58% of the French Meadows subdomain, but the median (±standard deviation) RAIM increase was negligible at 0.06 ± 1.03%.

| Snowpack response to thinning in each subdomain
Most subdomains show similar patterns between the light and heavy thinning scenarios, with the heavy thinning scenario resulting in larger magnitude changes in RAIM than the light thinning scenario for regions with both increasing and decreasing RAIM.Averaged over the three water years, we see a 2.8% (light thinning) and 6.7% (heavy thinning) median RAIM increase in the West Shore, 2.4% and 4.5% increase in Sagehen, 0.1% increase and 7.0% decrease in the Yuba, 0.5% and 4.8% decrease in the East Shore, and 0.1% and 0.1% increase in French Meadows, respectively.However, there are different patterns in the Yuba between thinning scenarios that highlight the processes controlling ΔRAIM.Specifically, in WY17, there is an average ΔRAIM of 1 ± 37 mm for the light thinning scenario, but an average ΔRAIM of À217 mm in the heavy thinning scenario.We believe that this is largely because WY17 was such a wet year in the Yuba (total precipitation of 2376 mm, as compared with a WY16 total of 1346 mm and WY18 total of 1127 mm), of which 92% of the precipitation fell as snow, as opposed to WY16 (82%) and WY18 (72% as snow).During this abnormally cold and wet year, the average wind speed was also 15-20% higher than WY16 and WY18.We found that removing vegetation in the heavy thinning scenario created large enough clearings to increase radiation and sublimation, in addition to blowing snow out of the Yuba subdomain entirely.These effects reduced the RAIM more for the heavy thinning scenario than for the light thinning scenario.We calculated median wind redistribution values of 0.9 ± 6.1% during the control scenario in the Yuba, 1.1 ± 6.7% during the light thinning scenario, and over 11.3 ± 11.2% during the heavy thinning scenario (Table 2).We believe that since WY17 was wet, cold, and windy, a much larger percentage of the falling snow was blown out of the Yuba subdomain than during other water years or study areas.Other subdomains had smaller initial clearings and lower wind speeds, which probably buffered this effect and led to less snow being blown out of the study areas.
We see a clear relationship between ΔRAIM with vegetation height, vegetation density, and LAI in some subdomains, but not others (Figure 5).Specifically, we see large positive ΔRAIM in Sagehen and the West Shore for most vegetation heights over 7 m tall and densities of at least 40% canopy cover.However, in the Yuba, we see positive ΔRAIM for only denser (>50% canopy cover) and taller (>14 m) vegetation and negative ΔRAIM for shorter and less dense vegetation.Again, we observe little ΔRAIM in response to thinning in French Meadows or the East Shore where the vegetation is shorter and the snowpack is smaller and less responsive to thinning than in other subdomains.Our results build upon previous explanations of Differences in climate, topography, and land use history dictate location-specific forest species and structure that cannot be captured with a single relationship.Our results are identical to Harpold et al. (2020) and Krogh et al. (2020) for the West Shore subdomain where the response of RAIM across forest height and density was relatively consistent (Figure 4).The results for the West Shore appear to be transferrable to Sagehen, which is a similar elevation and climate.However, we observed more variable RAIM responses to thinning in dense and tall vegetation in the Yuba, with negative ΔRAIM in lowdensity forest patches across the East Shore, Yuba, and French Meadows (Figure 4).The results show that observations from a single subdomain may not be transferrable to the others, particularly when the sites vary significantly in either their climate or pre-thinning forest structure (Figure 5).
We partially repeated the analysis of Varhola et al. (2010) by correlating the changes in fractional forest cover with resulting changes in SWE accumulation and ablation due to thinning (Figure 6a).Previous cross-site synthesis efforts (Varhola et al., 2010) suggested that snow accumulation should increase by 4% for each 10% decrease in forest fraction.Our results from the West Shore showed only a 1.1% increase in snow accumulation (Krogh et al., 2020) for each 10% decrease in forest cover (Figure 6a).We observed similar trends for Sagehen but a lower gain of 0.1% per 10% reduction in forest cover.
However, the Yuba and East Shore showed the opposite trend, where removal of forest cover resulted in decreases in average snow accumulation in both subdomains.The trend for French Meadows was statistically insignificant, which corresponded with the small magnitudes of ΔRAIM (Figure 4).The ablation relationships with forest cover removal were generally weaker than accumulation and showed large variability.Sagehen and the West Shore showed increased ablation with forest removal, which is opposite of the general trend from Varhola et al. (2010).However, the East Shore, Yuba, and French Meadows showed decreasing ablation with forest cover removal.The East Shore had the largest negative response of SWE and ΔRAIM to forest restoration compared with the other subdomains.The East Shore had much higher shortwave radiation (133.8%),shorter average vegetation height (26.6%), windier conditions (130%), and thinner canopy density (37.4%) compared with the average of the other subdomains.These losses of available water from forest restoration highlight the strong influence of snow evaporation and sublimation (which dominate the overall mass balance) on snow and water dynamics, despite the influence of canopy shading in the energy balance.
F I G U R E 5 Relationships between ΔRAIM, initial vegetation height, and vegetation density for each subdomain based on WY18 (the driest year), as compared with changes in RAIM for low, medium, and high LAI value pixels during the heavy thinning scenario.Symbols are coloured by the percent change in RAIM, and symbol size represents the control LAI.Note how ΔRAIM is the highest for different vegetation structures for each subdomain, implying different management recommendations across sites.
F I G U R E 6 Relative changes in (a) snow accumulation and (b) snow ablation from the East Shore (red), Sagehen (blue), West Shore (green), and Yuba (cyan) subdomains, with the 10th and 90th percentiles indicated in matching dashed lines.Note that removing vegetation in the West Shore and Sagehen leads to higher snow accumulation and ablation, while it leads to lower accumulation and ablation in the East Shore, French Meadows, and Yuba.

| Using RF results for a landscape-scale decision support tool
We used the ΔRAIM calculations to calibrate 30 m scale random forests (RF) models for each subdomain, independently, to inform the most important variables for a landscape-scale analysis.We trained each RF using topographic and climatic variables to make the models more transferable to other locations: canopy density, ΔLAI, elevation, slope, northness, eastness, vegetation height, and mean wintertime (December-February) air temperature, precipitation, and wind speed.
The out-of-bag RMSEs for site-level models were quite small, suggesting that the RFs achieve sufficient model performance to accurately capture the variability and heterogeneity of ΔRAIM across the three water years and two thinning scenarios.We observed overall RMSE of 3.77% (Sagehen), 7.82% (Yuba), 4.99% (West Shore), 9.66% (East Shore), and 0.95% (French Meadows).Not surprisingly, ΔLAI was the most important variable for all subdomains, on average, while incoming shortwave radiation was the second most important variable everywhere except French Meadows, where we did not see a clear pattern of variable importance.Overall, topographic variables had relatively little effect on the RF models, since they were highly correlated with air temperature, precipitation, and incoming radiation, which capture some of the same processes.
We expand on the RF-type approach from Krogh et al. (2020) using the complex SnowPALM outputs for each domain to scale to the 6300 km 2 TCSI restoration area.However, the variability in RAIM in response to thinning across subdomains (Figure 4) necessitates an approach that can leverage both climate and forest structure information beyond that of Krogh et al. (2020).We developed another RF model that leveraged the large response of sites and thinning scenarios using the results from all five subdomains and all three water years (Figure 7).We calculated a TCSI-scale out-of-bag RMSE of 6.1%, which is larger than the mean RMSE of the four site-specific RMSEs (3.4%), but low enough to provide confidence for prediction.To determine the efficacy of our TCSI-scale RF, we performed a leave-one-out cross validation, where we trained the RF model on four of the five subdomains and calculated the ΔRAIM at the fifth subdomain, as compared with the SnowPALM ΔRAIM generated at that subdomain (Figure 8).We calculated median percent errors of 24.2% (West shore), À1.4% (Yuba), À2.4% (Sagehen), 5.0% (East Shore), and À18.7% (French Meadows).Note that since ΔRAIM is so small in parts of the East Shore and French Meadows, the percent error appears more inflated.Since we calculate reasonable leave-one-out cross validation errors and low RMSEs, we use the calibrated TCSIscale RF to predict the effects of vegetation thinning on RAIM across the TCSI region.The RF model at the TCSI landscape scale gives information on the importance of each predictor variable.The most important controls on landscape-scale response were wind (27.8) and ΔLAI (20.4), followed by precipitation (12.3) and pre-thinning LAI (9.3), while the remaining variables were relatively unimportant (Figure 7a).
Note that precipitation, wind, and air temperature alone account for over 46.5% of the importance.Average wintertime windspeed was a stronger predictor of ΔRAIM than ΔLAI, despite a relatively small range of average wind speeds across the subdomains (3.53 to 5.09 m/s; Table 1).Subdomains that had the highest ΔRAIM also had the lowest wind scour, while those that had the most negative ΔRAIM also had the highest amount of wind scour.Specifically, we see the largest RAIM increases in the West Shore and Sagehen, where the wind scour index only increases by 1.2% and 2.4% between the control and heavy thinning scenario, respectively, while we see the largest RAIM decreases in the East Shore, where the wind scour increased by 9.5% due to vegetation removal (Table 2).
The limitations of SnowPALM modelling are worth noting for future studies.Lidar data to run the model are limited in spatial and temporal extent, while climate forcing data are less certain over mountainous regions far from weather stations used for calibration.
heavy thinning scenarios, respectively.ΔRAIM is negative in response to thinning on the western slope of the Sierra Nevada, especially at lower elevations, as well as along the east shore of Lake Tahoe and in the northwest corner of the TCSI region.On the other hand, we observe large positive ΔRAIM along the crest and leeward side of the Sierra Nevada and along the west and south shores of Lake Tahoe for both thinning scenarios (Figure 7b,c).The largest positive ΔRAIM are concentrated immediately east of the Sierra Crest, coincident with the highest snow levels and overstocked, fire-suppressed forests.Both positive and negative ΔRAIM are increased when more forest cover is removed in the heavy thinning scenario.
Our results differ with one of the few existing landscape-scale

1
Map showing subdomains of Yuba, Sagehen, French Meadows (FM), East Shore (ES), and West Shore (WS) in the Tahoe Central Sierra Initiative (black outline) region, as well as the towns of Truckee and South Lake Tahoe.Background image is the US Geological Survey 30 m digital elevation model.ecologically diverse TCSI landscape.This knowledge gap highlights the importance of developing both understanding and useful tools to guide management strategies across diverse landscapes.Despite the widespread practice of forest restoration, land managers have few decision support tools to guide management strategies to enhance hydrologic benefits of forest restoration at landscape scales (>100 km 2 ).For example, Dickerson-Lange et al. (2021) developed a decision system to help landowners determine where and how forests can be managed to optimize water storage.Specifically, they suggest different silvicultural approaches towards enhanced water storage based on the climatic regime: heavy thinning resulting in large canopy gaps for warmer regions with early snowmelt; light thinning resulting in small canopy gaps for colder climates; and denser forest cover for windy locations.The framework from Dickerson-Lange et al. (2021) provides an excellent starting point for determining where to selectively thin within the TCSI based on climate alone, without considering vegetation characteristics.
. Similarly, a stakeholder-driven process with broad representation across different interest groups in the central Sierra, in partnership with the US Forest Service (USFS) and the University of California, Berkeley, designed and implemented the Sagehen Forest Restoration Project.The project mechanically or hand thinned 1040 ha across the Sagehen Creek Experimental Forest (SCEF) watershed to reduce the risk of largescale, high severity fire and improve forest resilience, including habitat for the forest-dependent American marten (Martes Americana; USDA et al., 2013).In the Yuba River watershed, the North Yuba Forest Partnership proposed to selectively thin 6000 ha through the Yuba Forest Resilience Bond, in hopes of managing the forest sustainably and responsibly, despite a changing climate and increased wildfire.
models for each of the five subdomains to determine the most important metrics for predicting the change in RAIM from forest restoration and to inform land managers.For the site-specific RFs, we analysed the effect of elevation, slope, northness, eastness, vegetation height, air temperature, precipitation, wind speed, LAI, and change in LAI (ΔLAI) on the response variable ΔRAIM across the three water years and two thinning scenarios.

F
I G U R E 3 Example 10 (light) and 20 m (heavy) thinning scenarios in the Sagehen subdomain.(a) Histogram of vegetation height before and after 10 and 20 m virtual thinning.(b) Map showing areas where operability is compatible (shaded regions) with removing trees >10 m in height, within 300 m of established roads (black lines), and on slopes less than 35%.On the other hand, operability is compatible with removing trees <10 m in height across the entire region.F I G U R E 2 Calibration of SNOTEL sites for each subdomain, resulting in RMSE values at the West Shore of 106 mm, Sagehen of 209 mm, Yuba of 210 mm, and East Shore of 72.8 mm, while the WY15-17 snow depth RMSE for French Meadows is 133.8 mm.snow retention from Lundquist et al. (2013) and Safa et al. (2021), who found enhanced snow retention in lower density forests and forests with warmer winters.While most of the East Shore's response to thinning results in negative ΔRAIM, we do see some positive ΔRAIM values for the taller and denser vegetation.The differences in both vegetation height and density, as well as climate, between these five subdomains highlight the complex relationship between forest restoration and snowpack response.

F
I G U R E 4 Histograms of changes in rain + snowmelt (ΔRAIM) resulting from thinning for the light (left column, a) and heavy (right column, b) thinning scenarios across all subdomains and years.Note that the magnitude and sign of ΔRAIM can change among water years and between thinning scenarios.
Importance of predictor variables from the TCSI-scale random forests model on calculating ΔRAIM.Random forests predictions for the TCSI region encompassing all available lidar for the (b) light and (c) heavy thinning scenarios.wouldn't necessarily remove all vegetation shorter than 10 or 20 m in height.While our results are applicable across the central Sierra Nevada, our recommendations have limited transferability outside of the TCSI, particularly limited in different climate or forest ecosystems.However, we demonstrate how snow dynamics modelling can be strategically applied to subdomains and subsequently scaled up to inform landscape restoration being conducted across large regions.

3. 4 |
A landscape-scale decision support tool for TCSI At the landscape scale, important patterns emerge in both the light and heavy thinning scenarios that could inform larger scale forest restoration planning.Overall, we predict positive ΔRAIM due to forest restoration across 64.1% and 69.1% of the TCSI, with an average ΔRAIM of 50.9 ± 132.1 and 70.4 ± 132.6 mm a À1 for the light and T A B L E 1 Study site characteristics.
decision support tools developed by Dickerson-Lange et al. (2021) that uses local climate controls to predict snowpack response to thinning.The Dickerson-Lange et al. (2021) results indicate that light (for "Cold" classifications, $40% of the TCSI) to heavy forest restoration (for "Warm" classifications, $55% of the TCSI) should lead to increased snow retention.In other words, their model suggests that <5% of the TCSI would see negative snowpack response to thinning.However, using the SnowPALM-trained RF model, we find negative ΔRAIM across 35.9% and 30.9% of the TCSI region, for the light and heavy thinning scenarios, respectively.Dickerson-Lange et al. (2021) classify each region by the average December-February temperature and wind speed, as well as the maximum April-May temperature.For each of the Dickerson-Lange et al. (2021) classifications, we see regions where SnowPALM predicts both positive and negative ΔRAIM from forest restoration (Figure 9b).For example, the "Cold & Late" classification is bimodal, with both large positive and negative ΔRAIM values depending on where vegetation is removed.Additionally, for each of their classifications, we predict statistically significant positive and negative ΔRAIM across the TCSI.While Dickerson-Lange et al. (2021) patterns generally capture the climate variability that highlights the largest possible response along the Sierra crest (which they classify as "Cold and Late"), they also predict "Windy" classifications scattered throughout the TCSI, rather than the windy locations concentrated leeward of the crest in the eastern part of the TCSI (Figure 10).The difference in wind-prone classification could be because our study used NLDAS hourly wind speed, while Dickerson-Lange et al. (2021) used topographic convexity of the digital elevation model as a proxy for wind exposure.The differences from Dickerson-Lange et al. (2021) highlight the potential benefits from including vegetation thinning scenarios and forest structure in our SnowPALM-derived results, as well as the complexities of estimating a robust landscape-scale decision support tool.The RF model explicitly incorporates forest structure information, making it directly applicable to our study area.We demonstrate the importance of forest structure in controlling response at most of the subdomains (Figure4).The reliance on lidar datasets is also aF I G U R E 9 (a)Dickerson-Lange classification for the TCSI region with subdomains overlaid.(b) Histograms of SnowPALM RF results for the light (grey) and heavy (red) thinning scenarios for each of the Dickerson-Lange classifications.challenge because large lidar extents are expensive and timedemanding to collect.A second advantage of our decision supporttool is that it is trained and validated on smaller spatial extents, specifically our five subdomains.However, we have limited options to validate the thinning effects of our model (besides comparing with previous results likeVarhola et al., 2010).As large landscape forest restorations become more common, constructing and validating decision support tools will be critical to making the best and most costeffective management decisions.We believe our decision support tool and its development will meet the needs of the TCSI project and provide a possible template for other areas.4| CONCLUSIONSThe pace and scale of forest restoration is increasing across the western United States in an attempt to prepare the landscapes for changing climates and to avoid the worst impacts from droughtinduced mortality, wildfire, and loss of water yields.Our results using the high-resolution SnowPALM model show that interactions of climate and vegetation structure lead to both increases and decreases in snow and RAIM.That can depend on inter-annual climate.The resulting decision support tool highlights both larger patterns of forest removal increasing and decreasing RAIM, as well as the potential for finer scale information associated with local terrain and forest structure.In this part of the Sierra Nevada, managing forests to reduce fuels can contradict objectives to improve snowpack retention.We calculated that the most important variables for predicting increases in snowmelt following thinning were lower wind speed, higher ΔLAI, greater precipitation, and higher pre-thinning LAI.Our results lead to different management recommendations than other empirical models (i.e.,Varhola et al., 2010) or decision support tools based on gridded climate data (i.e.,Dickerson-Lange et al., 2021).Specifically, we show large negative ΔRAIM from thinning and more influence of local forest structure than these other available sources of information.Variability in forest treatments and other natural disturbances, like insects and wildfire, make our results hard to generalize but suggest the benefit of using high-resolution models to examine potential benefits of forest restoration.Moreover, the fate of RAIM after it reaches the ground surface, whether it infiltrates into the soil and flows down to a stream, is likely site specific and would require further hydrologic observations and modelling.Despite decades of research on the topic of post-forest thinning snowpack response, these compensating and cross-scale processes remain challenging to predict, validate, and communicate for improved forest restoration planning. .