Meteorological data come from SNOTEL stations and NLDAS Phase-2 (Xia et al., 2012; NLDAS project, 2021) downscaled spatially and to account for local topography and forest structure using built-in routines in SnowPALM (Broxton et al., 2014). Sites 3 and 4 are collocated with SNOTEL stations 308 and 617, respectively, which provide local air temperature and precipitation. For sites 1 and 2, air temperature and precipitation forcing data are developed using nearby SNOTEL measurements (Table 1) adjusted using regression relationships with locally-measured values of precipitation and air temperature made over 1–2 years (Geonor T200B wind-fenced weighing precipitation gage, radiation-shielded Type T thermocouple). At all sites, relative humidity, barometric pressure, wind speed, wind direction, downward shortwave and longwave radiation, and potential evapotranspiration come from the downscaled NLDAS data. While site-level weather forcing data are generally preferred, NLDAS data downscaled within SnowPALM agreed well with the available ground measurements of air temperature, humidity, and precipitation. Furthermore, model calibration to Snowtography-observed snow depths showed that NLDAS-data can be useful as model forcing to simulate snowpack dynamics at the high spatial resolution required to diagnose forest structural effects (see results below).

At each site, a network of automated trail cameras is maintained (Moultrie model M-999i, Campark model T180, and Wasoda Trail Camera model G300, all with time-lapse programmability) collecting photographs of snow depth measurement stakes (Snowtography) to generate daily time series of snow depth (Broxton et al., 2019, 2020; Payton et al., 2021; Dwivedi et al., 2022). Snowtography stakes (n = 16–29 per site) are spaced 5 m apart along transects that span forest stand and gaps in cardinal directions to capture spatiotemporal differences of snow accumulation and ablation due to differences in shading (along north–south transects) and wind (along east–west transects). Snow depth at each stake is measured using repeat photography that is interpreted manually. Snow stakes are marked with alternating red and white bandsone inch tall (~2.5 cm), and the estimated uncertainty of the measurement is ±1 cm. During 1–2 days after snowstorms, small snowdrifts (1–3 cm) can occur against the snow stakes, while during prolonged dry periods, small depressions (1–3 cm) can form around each stake, presumably due to altered radiation balance. In such cases, technicians are directed to interpret the snow depth of the neighborhood (~1 m radius, congruent with the resolution of SnowPALM modeling). Snowtography measurements began between 2017 and 2020 and run through winter 2022 for this study, providing 3–6 years of daily snow depths per site (Table 1).

Airborne lidar data were mostly collected using low-altitude airplanes. Snow-off lidar data were collected by Quantum Spatial, Inc. between 2013 and 2020 (Quantum Spatial Inc., 2013, 2014, 2020). Snow-on lidar collections were conducted on February 1, 2017, March 7, 2017, and March 4, 2019 at sites 2, 3, and 4 (van Leeuwen and Quantum Spatial Inc., 2017; Quantum Spatial Inc., 2019). At Site 1, snow-on lidar was collected on February 10, 2023 using a DJI Matrice 600 Pro drone with a DJI Zenmuse L1 Lidar+RGB sensor. Following our previous studies (Broxton et al., 2019; Dwivedi et al., 2022), maps of canopy height, canopy cover, and ground elevation were prepared at 1-m spatial resolution using Fusion software (McGaughey, 2020). Snow depth maps were created by differencing the ground returns from snow-on and snow-off lidar collections. Lidar-based spatial snow depth maps were previously bias corrected using manually sampled snow depths during field campaigns (Broxton et al., 2019). On the lidar flight dates, RMSE between lidar and Snowtography depths ranged from 7 to 13 cm (Supplementary Figure S2).

Previous works suggest several metrics for classifying snowpack stability (Sturm et al., 1995; Fierz et al., 2009; Hammond et al., 2018; Petersky and Harpold, 2018; Musselman et al., 2021). Their application to the present sites was challenging due to lack of data such as snow grain size, shape, microstructure, hardness, liquid water content, and high-resolution satellite data. Here a simple snow persistence metric (SPM) was defined based on daily SWE that is the ratio of days with minimal snowpack ablation (<5 mm/day) between the first and last days with SWE > 2.5 mm (Uecker et al., 2020). High SPM indicates persistent, seasonal snowpack while low SPM indicates ephemeral snowpack. Cumulative and mean radiation values for each model pixel were calculated using days with snow cover, which varies by model pixel, to quantify the energy balance over snowpack. To assess net hydrologic impacts of forest structure, mass fluxes were computed by SnowPALM daily and then aggregated at each pixel over a winter duration computed for each site and year from the first day of snow cover ≥99% of site pixels until the last day of snow cover ≥1%, where snow cover is defined by SWE > 2.5 mm. In this paper, the traditional hydrologic year definition is applied, meaning that references to the day of year are counted from October 1.

Following Mazzotti et al. (2020a), each 1-m pixel at each study site was classified into one of five snow environments created by the structure of surrounding forest: Under canopy (UC): canopy cover >10% for canopy height > 2 m; Cool Edge (CE): uncovered but shaded by nearby canopy to the south; Warm Edge (WE): uncovered, exposed to solar radiation from the south, and near canopy to the north; Overlapping Edge (OE): small forest gaps with nearby canopy to both north and south; and Open (O): uncovered and far from canopy. This approach was refined for the present study sites using cool edge and warm edge lengths of 10 and 5 m, respectively, based on observed daily ablation rates following several significant snowstorms within our Snowtography datasets (Supplementary Figure S2). Therefore, pixels that were both within 10 m to the north of canopy and within 5 m to the south of canopy were classified as Overlapping Edge.

The SnowPALM model (Broxton et al., 2014) was used to represent 3-D mass and energy balance of the snowpack and canopy environment at 1 m spatially and a daily time step. SnowPALM recently assessed the impact of forest thinning and wildfire on snowpack in conifer forests (Harpold et al., 2020; Krogh et al., 2020; Moeser et al., 2020). Important features of SnowPALM include evaporation and sublimation from canopy, snowpack, and wet soils, geometrically explicit interactions of solar radiation with canopy height and density (i.e., shadowing), canopy temperature and resulting longwave radiation balance, and influences of local topography and canopy on wind distribution of snowfall, using an adaptation of the landscape wind index approach (Winstral and Marks, 2002). SnowPALM considers trees as having uniform cross-section for their full height, extending vertically from the canopy height model to the ground, with constant density (i.e., from 1 to 100% closure) determined from the canopy density map. In this sense, SnowPALM parameterizes trees based on the information available from airborne lidar, and it does not explicitly consider the 3-D canopy morphology over the full tree height, which is generally only feasible for smaller spatial extents (Russell et al., 2021; Hojatimalekshah et al., 2023). The model was calibrated using snow throughfall, snow depth, and canopy snow storage (for details see Section S1 in Supplementary material). Prior work comparing ground measurements georeferenced using differential GPS systems with airborne lidar snow depth maps suggested that the uncertainty in lateral positioning was <1 m lidar (Broxton et al., 2019, Figure 6), providing good confidence for matching lidar and snowtography data to model pixels. During ground-based field surveys at site 4, replicate samples of snow depth (n = 5) collected within a 1-m radius have an average standard deviation of 6–7 cm [among sets of five replicates indicated as black dots in Broxton et al. (2019), Figure 6]. This is considered an upper estimate of the uncertainty involved in representing a 1-m square model pixel with a single snowtography stake, because 1-m radius is larger than 1-m square and spatial variability in snow depth tends to be greatest at site 4 with the most persistent snowpack. Due to the complexity of the SnowPALM model, calibration was performed manually using the modified Kling-Gupta Efficiency (KGE’; Gupta et al., 2009) and the spatial efficiency metric (SPAEF; Koch et al., 2018) to assess model fit compared with time series measurements from Snowtography and lidar snow depth maps, respectively. Both KGE’ and SPAEF range from −∞ (for the worst fitting model) to 1 (for the best fitting model). The final trained model has a unique parameter set for each site that was held constant across all study years.

The trained models were run for winters 2017 through 2022 at all sites. For sites 1 and 2, where forest thinning occurred during the model simulation period but before the Snowtography or lidar measurements, the post-thinning forest structure from lidar was used (section 3.3). The effects of snow environment with respect to canopy were evaluated with one-way ANOVA and multiple comparison of group means (MATLAB functions anova1 and multcompare).

Daily liquid water input (LWI) was computed as the sum of rainfall that does not contribute to snowpack (Dwivedi et al., 2022) and snowmelt computed by the model. LWI and other mass fluxes (e.g., canopy, snowpack, and total sublimation) were only computed for winter, defined using snow cover (see section 4.1). At sites 3 and 4, LWI was previously estimated using 1-D point runs of the Snow-17 model at each Snowtography stake using 1-D canopy cover for 2018–2020 (Dwivedi et al., 2022). Here, the results of this 1-D Snow-17 LWI modeling are extrapolated across each site using canopy cover at each 1-m snow pixel. The site mean and spatial variability of LWI were then compared for the 1-D canopy-cover based estimates vs. the spatially explicit 3-D SnowPALM model for the warm/dry year 2018 and cool/wet year 2019.

The relative abundance of each of the five snow environments varies due to forest structure (Figure 2). Site 3 is the only site with no known disturbance for at least several decades and accordingly has the greatest fraction under canopy (78%). Much of the remaining area at site 3 is small gaps classified as overlapping edges (12%). Mechanical thinning occurred within the last 5–10 years throughout site 1a, the southwest portion of site 1b, and the northern half of site 1, which has resulted in a “clump/gap” forest structure at this site (Belmonte et al., 2019). Approximately half of sites 1 and 2 are covered by canopy, with the remaining area distributed among open area and cool, warm and overlapping edges. The southern portion of Site 4 was moderately affected by 2011 wildfire and has a similar mix of snow environments as sites 1 and 2. In the following analyses, results are summarized by snow environments for each site.

Snow physics and lidar mapping simulates the observed daily snowpack dynamics at the Snowtography measurement locations in each snow environment. The observed and simulated snow depth patterns for an intermediate winter (Figure 1), 2020, are shown (Figure 3, see Supplementary Figures S4, S5 for other years). Snow disappearance date is generally accurate within ±2 days. For several 2020 events at sites 3 and 4, modeled snow depth increases exceed observed increases, likely due to model overestimation of snow/rain fraction and/or underestimation of snowfall density. Fortunately, the resulting small data-model residuals tend to be consistent for all snow environments across a site, and therefore the model results remain useful for quantifying the impacts of 3-D forest structure. Ablation rates are generally well-predicted.

The trained SnowPALM model reproduces the spatial patterns of lidar snow depth snapshots (Figure 4) with SPAEF ranging from 0.43 to 0.53 at all sites except site 2, where spatial patterns in the modeled map were generally similar to the lidar, but relative accuracy (SPAEF = 0.25) is reduced in part by low snow depths (360 mm at site 2 as compared to 560–944 mm at other sites) (Figure 4; Supplementary Figure S6; Supplementary Table S1).

Among sites, incoming shortwave radiation at the snow surface is highest at site 1 and the northern half of site 2 (Table 2), where recent mechanical thinning has reduced canopy cover (Figure 2). Site average incoming longwave radiation at the snowpack surface varies with site air temperature (Figure 1; Table 1) and is therefore greatest at site 3 and least at site 4 (Table 2). Average net radiation generally declines from a maximum of 18 Wm⁻² to a minimum of 1 Wm⁻² from site 1 to site 4, although the warmer average air temperature at site 3 gives it a greater net radiation average than site 2 (11 vs. 7 Wm⁻², respectively).

Total site-average sublimation increases from 13 ±6% to 21 ± 6% across sites 1–4, suggesting greater total sublimation with increasingly persistent snowpack (Table 2). Canopy sublimation dominates at the lowest/most ephemeral site, while snowpack sublimation increasingly dominates the total loss with increasing site snow persistence. However, direct comparisons of sublimation at the site scale are complicated by differences in canopy cover (Table 1). When considering only pixels with at least 10% canopy cover, canopy sublimation consumes 12 ± 5% of winter precipitation at site 1a declining monotonically to 4 ± 3% at site 4 (Table 2). Meanwhile, snowpack sublimation increases from 9 ± 1% to 19 ± 5% across sites 1–4.

Across sites, average peak SWE increases from 22 to 26% of winter precipitation at the lower/warmer sites up to 48% at the highest/coldest site 4, consistent with increasing snow persistence from 0.32 to 0.66 (Table 2). The site-averaged snowpack persistence pattern is consistent with the SNOTEL-based pattern (Figure 1 and Section 5.1), but all peak SWE:P and SSM are lower, which is consistent with SNOTEL locations as places with most/longest snow (Molotch and Bales, 2005; Broxton et al., 2019; Dwivedi et al., 2022). However, LWI as a percentage of winter precipitation decreases across this same gradient from an average of 87–79%, contrasting with the idea that peak SWE is a good proxy for liquid water inputs in snowmelt-driven watersheds. Instead, LWI is effectively equal to winter precipitation minus total sublimation (i.e., LWI ~ Winter P – Total Sublimation, Table 2). The portion of LWI derived from snowmelt increases from 86 to 90% from site 1 to site 4 (the remaining LWI is winter rainfall). Average snow cover duration increases monotonically across the site gradient of snow persistence, nearly doubling from an average of 57 days at site 1 to 105 days at site 4.

This section describes spatial patterns observed in site-wide maps. Statistical comparisons of these variables by snow environment will be presented in the following section.

Spatial maps of winter-integrated radiation balance and sublimation from the trained SnowPALM model (Figure 5) are consistent with our snow environment classification (Figure 2). Due to tree shading, shortwave radiation is low under canopy, especially for dense, contiguous canopy (Figure 5, e.g., most of site 3, undisturbed portions of sites 2 and 4). Within gaps, shortwave radiation is partially attenuated by tree shading at the cool edges. Incoming longwave radiation is greatest under canopy and somewhat enhanced at canopy/gap edges, consistent with thermal emission to the snowpack from canopy warmed by the sun and inhibition of nighttime thermal emission from snowpack under canopy.

Spatial variability in net radiation shows the combined effects of shortwave shadows and longwave emission near canopy (Figure 5). At each site, the highest net radiation occurs at warm edges, where unshaded shortwave is additive with longwave from nearby canopy. Under-canopy net radiation is second-highest, followed by overlapping edges and cool edges.

To illustrate finer-scale details indicative of forest structural effects on snow processes, the end-members of snow persistence, sites 1a and 4, are the focus of Figures 6–9. Within a given site, 3-D forest structure induces large ranges of spatial variability in the number of snow-covered days, with up to 30–80 days difference across a site (Figure 6; Supplementary Figure S7). Temporally-averaged peak SWE ranges from ~10 to 45% of P at the least-persistent site (1) and from 10 to 60% at the most persistent site (4) across each site’s mosaic of snow environments. Liquid water inputs vary between ~70 and 90% of total winter precipitation. At site 1a, spatial patterns in LWI are strongly controlled by under canopy/uncovered, with relatively small effects of various uncovered environments (Figure 6), consistent with the dominant role of canopy sublimation and lesser role of snowpack sublimation at the ephemeral end of the site gradient (Table 2). At site 4, however, there are larger LWI gradients across open areas between nearby canopy (e.g., gap shown in Figure 6, bottom row, southwest corner of site 4), suggesting that spatial variability in snowpack sublimation increases with winter duration. Statistical comparisons are presented in the following section.

The spatial variability in snow dynamics and LWI created by forest structure is illustrated for sites 1a and 4 with cross sections traversing gaps from cool edges through open areas and ending at warm edges (see Figure 6, left column, for location of each cross-section). These cross sections are deliberately placed across gaps to illustrate the range of spatial variability related to nearby forest structure and are not meant to represent site-average conditions, which are summarized in a subsequent section. Surprisingly, LWI amount increases when moving northward from cool to warm edges, even though net radiation increases and peak SWE decreases (Figure 7). The magnitude of variability in LWI/P is larger at the most persistent snowpack site, increasing from 77% at the cool edge to 86% at the warm edge, as compared to a LWI/P increase of only about 2% along the transect at the ephemeral site. Snowpack sublimation, a dominant control of LWI spatial variability in gaps, increases with snow covered days (i.e., time exposure of the snowpack) but decreases with increasing net radiation, at both sites (Figure 7).

The time courses of net radiation and snow dynamics (Figure 8) illustrate mechanisms underlying the spatial variability across gaps (Figure 7). As expected, net radiation is highest at warm edges throughout the winter, due to direct shortwave exposure combining with thermal radiation from the warm, south sides of nearby trees. Net radiation differences between open areas and warm edges begin to emerge more strongly in the second half of the winter, likely because solar angles steepen, shortening tree shadows and exposing open areas to increasing direct shortwave radiation. Most of the cumulative winter snow sublimation differences between sites and among locations at site 4 occur during the ~30-day difference in snow-covered duration at the end of winter. LWI occurs frequently throughout the winter at site 1, consistent with ephemeral snowpack. At the high/cold site 4, open areas and cool edges have, on average, persistent seasonal snowpack, while warm edges tend to have snow dynamics similar to a low/warm site with ephemeral snowpack. The ephemeral snowpack characteristics of the warm edge starts around day 60 and lasts until around day 180 (Figure 8) for an average LWI rate of ~1.95 mm/day. The open and cool edge environments accumulate snow throughout winter and melt in one pulse during spring, with sustained high intensity input after day 150 and lasting to around day 200, with an average LWI rate of 3.52 mm/day during that period.

While Figures 5–8 highlight the range of spatial variability in snowpack, LWI, and underlying mass and energy fluxes, this section provides summaries of forest structural impacts aggregated across each site and indicates their statistical ranking. These differences are illustrated at Sites 1 and 4, the end-member sites with least and most persistent snowpack (Figure 9, see Supplementary Figure S9 for other sites).

The patterns of average peak SWE across snow environments are mostly consistent with H1 (Figure 9). At site 1, the ranking of average peak SWE was cool edge>open>overlapping edge ~ warm edge>under canopy. At site 4, the ranking was the same except that overlapping edge ~ open.

Variability of liquid water inputs across snow environments only partially supported H1. Like peak SWE, LWI was lowest in the under-canopy environment at both sites 1 and 4, due to greatest total sublimation (Figure 9; Supplementary Figures S8, S9), consistent with H1. At both sites, under-canopy locations consistently partitioned an additional 8–12% of winter P to sublimation (reducing LWI) as compared to uncovered environments. Among the uncovered environments, cool and overlapping edges showed the greatest sublimation and least LWI, opposite of the pattern in peak SWE and the order of LWI proposed in H1 (Figure 9). These differences were greater at the most persistent snowpack site (4) than the ephemeral site (1), consistent with greater spatial variability of snowpack sublimation across uncovered gaps at site 4 than site 1 (Figures 5, 7, 8).

At the ephemeral snowpack site 1, snow covered duration for the under-canopy environment was 10–24 days less than any of the four uncovered environments (Figure 9; Supplementary Table S2; Supplementary Figures S8, S9), with cool edges having the greatest snow-covered duration. At the most persistent snowpack site 4, snow-covered duration tended to be greatest for cool edges, least for warm edges, and overall similar to the total sublimation patterns for all four non-covered environments, consistent with the idea that longer integration time is an important control of spatial patterns in LWI (H2). Under-canopy snowpack at site 4 tended to persist for a similar or slightly longer duration than the open area, in contrast to reduced snow-covered duration under canopy at site 1.

This section compares model predictions for the spatial distribution and site-average LWI by 1-D and 3-D approaches for sites 3 and 4, where 1-D modeling was previously conducted (Dwivedi et al., 2022). 3-D forest modeling with SnowPALM consistently predicts ~30 mm less site-averaged LWI than the 1-D forest cover approach (Figure 10), likely because the Snow17 model used in the 1-D approach does not account for snowpack sublimation or other effects such as snow redistribution due to wind or longwave enhancement due to vegetation. Fractionally, this difference is more important in a warm/dry year than a cool/wet year. Aside from this bias, histograms of LWI across the study domains show roughly similar distributions during the warm/dry year at each site. During a cool/wet year, however, the 3-D model produces a more widely spread and bimodal distribution of LWI as compared to the unimodal, densely concentrated 1-D prediction, consistent with the prediction of H2 that LWI spatial variability regulated by canopy structure would be more important in the most persistent seasonal snowpacks and cold/wet winters.

3-D forest structure strongly regulated both canopy sublimation and net radiation to snowpack (Figures 5, 9; Supplementary Table S2), which only partially supports Hypothesis 1. As expected, peak SWE was greatest at cool edges (Figures 6, 7, 9), which have no canopy interception and lowest net radiation (Figures 5, 7), followed by open areas or overlapping edges, further supporting H1 (Figures 6, 7, 9). Critically, magnitudes of peak SWE were low compared to winter-integrated LWI at all but the most persistent snowpack site (Table 2; Figures 6, 7, 8), illustrating that winter rainfall and midwinter snowmelt act to disconnect peak seasonal SWE from LWI at ephemeral snowpack sites. This decoupling between peak SWE and LWI undermines a key assumption of our first hypothesis that sublimation would vary positively with net radiation. Instead, we found that LWI was often maximized at warm edges (Figures 6, 7, 9), where canopy sublimation ~0 (Figure 5) and net radiation are greatest (Figures 5, 7, 8). We infer that net radiation tended to drive earlier snowmelt (shorten winter duration), decreasing snowpack sublimation through decreased exposure time of snowpack to the atmosphere (Figure 8) thereby increasing LWI (Figure 7). For most environments at ephemeral snowpack sites, and for warm edges even at a persistent snowpack site, frequent snowmelt during winter both increased total LWI and advanced its seasonal timing (Figure 8).

Our second hypothesis that 3-D forest structure would have greater impacts on LWI at the most persistent snowpack sites and/or during cold/wet winters was supported across sites and years (Figures 7–10). Greater fractions of LWI deriving from snowmelt and less from rain at the higher/colder sites (Table 2) increases the importance of canopy structure, since snowfall is intercepted at greater rates than rain, and snowpack integrates canopy structure effects through modified energy environments. When the mosaic of environments created by 3-D forest structure was simplified into five classes, the relative importance of 3-D structure for LWI (as a fraction of winter P) was indeed more important at sites and in years when snow covered duration was longer (Figures 7, 10; Supplementary Table S2), although it was LWI at warm edges, not cool edges, that was enhanced. Accordingly, the additional complexity of 3-D modeling, when compared to the simple 1-D canopy cover approach to scaling snow processes, was more warranted in a cool/wet year with longer snow-covered duration (Figure 10).

Across the elevation gradient of snow persistence represented by these Arizona sites, canopy interception is an important control of partitioning winter precipitation between sublimation losses and liquid water inputs (Table 2; Figure 9). Therefore, removing canopy is expected to increase both peak SWE and LWI. However, uncovered areas have different snow dynamics and hydrologic partitioning depending upon the three-dimensional structure of nearby forest, and this matters more when winter is longer, especially during the final days and weeks of winter, when energy fluxes are highest (Figure 7; Supplementary Table S2). While we expected the greater net radiation at warm forest edges to drive greater sublimation (H1), we found instead more rapid snowmelt, decreasing exposure to sublimation (Figure 7).

This work has important implications for forest management in the interior United States Southwest, which is typified by thin, ephemeral snowpacks and high radiation inputs and where nearly all snowmelt occurs below tree line (Broxton et al., 2021). Our results suggest that forest management reducing canopy would benefit LWI amount by 5–11% of winter P when contrasting under-canopy with the mean of all other environments (Figure 9). As a simple thought experiment, consider the effects of reducing canopy cover from 66 to 33% in a watershed with 400 mm winter precipitation. Assuming an average 9% increase in LWI for the portion converted from canopy to any uncovered environment, LWI across the watershed is predicted to increase by 3% of winter P, or 12 mm. However, our results show that forest structure matters, at least when snow is relatively persistent, and that further LWI amount enhancements might be possible with post-thinning structure prioritizing an abundance of warm edges (Figures 8, 9). It should be noted, however, that LWI inputs on warm edges occur episodically in small pulses throughout winter, contrasting with the more concentrated and higher intensity inputs at cool and open edges at site 4 during spring. Therefore, it is unclear if the increased LWI at warm edges and warm sites enhances surface or subsurface runoff or if it is evaporated and transpired during the more frequent warm, snow-free periods. Furthermore, although the sites in this study do not have an elevational precipitation gradient, higher/colder sites with longer winters and more persistent snowpack usually have greater precipitation inputs, making them overall more important for streamflow generation.

The results of this study emphasize the potential of forest management to significantly affect snowmelt timing (Figure 9). Earlier snowmelt is expected to generate streamflow more efficiently because it occurs at a time of low evaporative demand and plant uptake (Knighton et al., 2020; Robles et al., 2020; Biederman et al., 2022; Gordon et al., 2022). Here, we found that reducing tree cover tended to delay snowmelt at warm sites with ephemeral snowpack but advance melt timing when snowpack is persistent (Figure 9). While this pattern of canopy having opposing effects on snow cover duration at warm vs. cool sites is consistent with prior frameworks, our finding that canopy prolongs snow cover at site 4, with mean DJF temperature of −0.1°C (Table 1), suggests a higher threshold for the switching behavior than previously reported (-1°C, Lundquist et al., 2013 or − 4.5°C, Safa et al., 2021). The spatial structure of the remaining forest has the potential to influence snowmelt timing, with snow cover at cool edges lasting longer than at warm edges by an average of 15 days (Figure 9) and up to 60 or more days (Figure 7). In large-scale efforts to restore forests to historically lower density and improved fire resilience through mechanical thinning and managed fire, such as the Four Forest Restoration Initiative,¹ forest cover reduction goals are relatively fixed and generally driven by fire-risk reduction (Robles et al., 2014). However, the present results suggest there is an opportunity to consider how the spatial arrangement of thinning influences the amount and timing of liquid water inputs. In watersheds similar to our study sites, thinning that prioritizes cool edges is likely to delay snowmelt, enhancing the “water tower” benefits of storing water resources at high elevation, shortening the duration of summer drought stress for vegetation (Belmonte et al., 2022; Sankey and Tatum, 2022; Dwivedi et al., 2023), but exposing snowpack to additional sublimation and reducing water input volumes. Thinning that instead prioritizes warm edge creation is suggested to advance snowmelt timing, decreasing sublimation losses and potentially enhancing streamflow through increasing of LWI and decoupling of water and energy availability, although this would depend in part on the rate of snowmelt, with higher snowmelt rates (e.g., mm/day) producing greater streamflow per unit of SWE (Barnhart et al., 2016). Even if warm edges enhance LWI to the root zone, earlier conclusion of snowmelt may increase vegetation water stress (Dwivedi et al., 2023). Across large landscapes, treated areas that result in warm versus cool edge environments may be of similar magnitudes, suggesting that these impacts may counteract one another. However, we are pursuing further research to evaluate if variations in forest treatment attributes such as orientation of openings or perimeter to area ratios can result in meaningful shifts in hydrological outcomes at the landscape scale.

This snow-focused work has important implications for other ecohydrological outcomes. Coupling of these snow analyses with soil moisture measurements is important to evaluate the accuracy of the LWI time series and quantify the impacts of these moisture inputs on soil moisture. We recently conducted such research with a limited number of soil profiles and a one-dimensional soil moisture model to estimate drought stress and deep percolation (Dwivedi et al., 2023), but modeling work across hillslope at catchment scales is a vital step toward understanding hydrologic feedbacks at the hillslope scale, including plant-water relations (Tai et al., 2017, 2018; Thomas et al., 2024) and larger-scale streamflow prediction (Safeeq et al., 2013; Watts et al., 2016). Second, forest managers urgently seek assessment of how forest thinning affects successional trajectories, fuels moisture, and seedling success, all of which could be accomplished through coupled studies of snow, soil moisture, and forest ecology (Jackson Leonard, USFS Rocky Mountain Research Station, personal communication Dec 1, 2023).

There are several important avenues for advancement of this work. First, while snow depth is generally more variable than density, density becomes an increasingly important source of uncertainty when snow depth is shallow (Raleigh and Small, 2017), and field measurements of snow density can be used with machine learning to map density (Broxton et al., 2019). Second, simplistic snow environment definitions (e.g., ≤ 10 m north of any pixel with canopy density ≥ 10% is considered a cool edge) ignore the important complexities observed on the ground. Because radiation penetrates some distance northward from a gap into the forest, the “warm edge” is most likely a transitional environment rather than being truncated at the forest edge (Supplementary Figure S3). In more windy locations, wind effects on spatial distribution of snow have been shown to occur over length scales of up to 55 and 98 m on the windward and leeward sides of forest stands (Webb et al., 2020). Snowtography stations installed since 2023 include elements to characterize the length scale of warm edges in various forest types including across north–south transitions, intended to quantify radiation length scales, and along east–west transitions, to capture wind effects along the prevailing wind direction. Third, there is a need to expand this work to colder, snowier, and windier sites at higher elevations and latitudes (we established Snowtography sites in Colorado and Wyoming, United States, in 2021–2023). It is possible that in more energy-limited environments, snowpack will be relatively persistent across the forest structural mosaic, and that instead of driving earlier melt, the increased radiation at warm edges will indeed drive greater sublimation, consistent with our Hypothesis 1. Finally, mostly level sites were selected for this study to isolate canopy effects, whereas slope and aspect are strong controls on snow processes (Broxton et al., 2014, 2021).