Analysis of topographic controls on depletion curves derived from airborne lidar snow depth data

The spatial distribution of snow water equivalent (SWE) is affected by numerous processes at multiple scales including, but not limited to, orography, wind and avalanche redistribution, and ablation dynamics driven by the land-surface energy balance. Inherently, these factors are strongly influenced by topography and vegetation and therefore result in annually repeating spatial patterns (Kirnbauer & Bloschl 1994; König & Sturm 1998; Deems et al. 2008; Sturm & Wagner 2010). Our ability to characterize these persistent patterns of snow cover and their relationships with snow accumulation have direct relevance for forecasting water supplies (Potts 1937; Hannaford et al. 1979; Kolberg & Gottschalk 2010) and estimating land–atmosphere energy exchange and atmospheric processes (Flanner & Zender 2006). In addition, patterns of snow distribution have relevance for understanding the distribution of soil moisture (Harpold & Molotch 2015), ecosystem response to water availability (Trujillo et al. 2012), and wintertime animal movement and habitat (Parker et al. 1984; Bhattacharyya et al. 2014). Depletion curves that link the snow-covered area (SCA) with SWE are a frequently used method for accurately estimating the contributing area to snowmelt runoff (Anderson 1973; Luce et al. 1999) and for characterizing snow-albedo feedbacks within land-surface and general circulation models (Niu & Yang 2007; Swenson & Lawrence 2012).

The existing literature regarding depletion curves typically relates SWE and SCA using empirical relationships or probability distribution functions for modeling applications. These cover a variety of spatial scales from point-based derivations over hillslopes (Anderson 1973; Luce et al. 1999; Luce & Tarboton 2004), physically based gridded models applied to the watershed scale (Clark et al. 2011), and regional to global-scale land-surface models (Niu & Yang 2007; Swenson & Lawrence 2012; Driscoll et al. 2017). However, the scale at which snow heterogeneity is consistent is ambiguous, prompting statistical parameterizations of depletion curves based on the coefficient of variation (CV) of SWE (Liston 2004). Estimates of the CV of SWE and snow depth for different environments have been provided in the literature based on field studies (Liston 2004; Clark et al. 2011); however, these are broad classifications that may not hold true at horizontal scales larger than 100 m, in mixed landscapes, or for elevation differences larger than 200 m (Essery & Pomeroy 2004; Clark et al. 2011). Importantly, satellite-based methods of SWE estimation have limited utility resolving the CV of SWE and snow depth, given that the accuracy of SWE estimates is highly sensitive to SWE spatial variability (Vander Jagt et al. 2013).

Multiple studies have investigated sub-grid depletion curves of SWE based on terrain characteristics (Donald et al. 1995; Niu & Yang 2007; Swenson & Lawrence 2012; Helbig et al. 2015; Fassnacht et al. 2016; Driscoll et al. 2017). In particular, Niu & Yang (2007) and Swenson & Lawrence (2012) observed that snow cover depletion patterns across the USA varied substantially between cells with differing sub-grid variances in elevation. Accordingly, they developed scale-variant SWE depletion curve parameterizations based on sub-grid elevation variability. Topographic controls on SWE–SCA relationships exhibit significant interannual consistency (Fassnacht et al. 2016; Driscoll et al. 2017), suggesting that the shape of depletion curves is largely controlled by the terrain that affects both snow accumulation and ablation.

Few studies have explicitly investigated the relationship between snow depth and SCA, perhaps because snow depth has been perceived as less relevant than SWE for estimating water supply. Helbig et al. (2015) parameterized a snow depth depletion curve with the sub-grid standard deviation of snow depth at peak accumulation. Similarly, Egli et al. (2012) indicated that estimates of maximum snow depth distribution were critical to model the evolution of basin mean snow depth and SCA throughout the melt season. López-Moreno et al. (2014) observed a consistent spatial pattern of the CV of snow depth across multiple lidar acquisitions in time, and Egli & Jonas (2009) found that variation in year-to-year snow depth depletion was largely controlled by the standard deviation of snow depth at peak snow depth. Other works have also related the standard deviation of snow depth to the variability of the underlying terrain (Marchand & Killingtveit 2005; López-Moreno et al. 2014; Helbig et al. 2015). These studies suggest that topographic variability that affects the mean and standard deviation of snow depth will influence the relationship between snow depth and SCA.

These previous studies were largely based on theory, modeling results, or relatively sparse measurements of the snowpack in time or space because regular, spatially extensive measurements such as those of the U.S. National Aeronautics and Space Administration (NASA) Airborne Snow Observatory (ASO) to assess depletion curve characteristics have not been available until the last 7 years. Hence, we still lack a process-level understanding of the spatial variability in depletion curve characteristics and how individual terrain elements at scales less than 100 m affect this behavior.

The NASA ASO dataset offers a new opportunity to observe and investigate spatial differences in depletion curves that result from topographic influences on snow accumulation and ablation processes (Painter et al. 2016). The ASO dataset provides high spatial resolution time series of lidar-derived snow depth and hyperspectral measurements of snow properties over multiple years from which relationships between SCA and snow depth can be determined. ASO also produces a SWE dataset that combines the lidar-derived snow depth with modeled density, but we focus on snow depth depletion curves in this study to keep the analysis anchored by direct measurements. However, improvements in characterization of snow depth depletion curve variability are relevant to the estimation of basin SWE since snow depth varies an order of magnitude more than snow density (Mizukami & Perica 2008).

Our objective is to gain insight into the processes by which snow and terrain interact to produce differences in SCA. We do not attempt to parameterize a universal depletion curve but rather to assess the influence of terrain characteristics on depletion curve shape. Specifically, we ask: what are the primary physiographic controls on the characteristics of snow depth depletion curves?

This study was conducted on data collected by the ASO in the Tuolumne River basin in the Sierra Nevada Mountains in CA, USA (Figure 1). The basin area is 1,175 km² with land cover that is 48% vegetation, 50% rock/talus, 2% open water, and << 1% permanent snow/ice. The elevation range of 1,127–3,997 m encompasses four distinct ecological zones: the lower montane forest; upper montane forest; a sub-alpine mixed-meadow/coniferous forest zone; and an alpine zone with large granitic features, talus slopes, and limited herbaceous vegetation (NPS 2016).

The climate is characterized as Mediterranean with relatively mild winter temperatures and dry, warm summers. The snow season typically lasts from November to March with more than 70% of the yearly precipitation falling as snow. A long-term record of SWE since 1967 (snow depth is only available since 2004) from 17 snow pillows within 20 km of the Tuolumne Basin recorded a mean peak SWE of 0.8 m with a range of 0.3–1.5 m. The four study years were characterized by below average snowpack with 2014 and 2015 characterized as a severe dry snow drought (Harpold et al. 2017), experiencing only 36 and 35% of the climatological mean peak SWE, respectively. The years 2013 and 2016 also experienced below average snowpack conditions, but less severely, with the data from the snow pillows reporting 62 and 71% of climatological peak SWE (data available from http://cdec.water.ca.gov/). Runoff from the watershed flows into the Hetch Hetchy Reservoir, the water supply for the City of San Francisco and other Bay Area municipalities.

We use the ASO dataset which is derived from a paired scanning lidar and imaging spectrometer (Painter et al. 2016). We use the ASO 3 m resolution snow-free digital elevation model (DEM), 3 m vegetation height, and 3 m resolution snow depth datasets. The ASO flew approximately weekly from near peak SWE through the snowmelt season, with six flights in 2013, ten in 2014, eight in 2015, and eight in 2016. Painter et al. (2016) report a mean absolute vertical error of 8 cm and a bias of <1 cm at 3 m spatial resolution when compared with manually measured snow depths. Further details about ASO and data processing can be found in Painter et al. (2016).

The 3 m resolution gridded snow depth data are aggregated to two new 495 m (nominally 500 m) resolution datasets, containing the mean and standard deviation of snow depth for each 500 m grid cell, respectively. Similarly, fractional SCA (fSCA) maps at 500-m resolution are derived from the 3 m gridded snow depth data based on a binary snow-free/snow-on schema. Pixels at the edge of the domain were not included if the 500 m grid cell extended past the boundary of the basin and pixels containing water bodies were also removed. We chose 500 m resolution for the analyses herein because it is the scale at which daily satellite observations of fSCA exist from the MODerate-resolution Imaging Spectroradiometer (MODIS); hence, our results are theoretically transferable to applications using MODIS data. We consider the ASO snow-free classification to be robust, because the original 3 m gridded snow depth data were masked for snow-free areas based on spectral data from the spectrometer (Painter et al. 2016). ASO SCA has also compared well against SCA derived from 0.3 m resolution World View imagery, with Recall, Precision, Accuracy, and F statistics all >0.99 (Bair et al. 2016).

We derived depletion curves for each 500 m grid cell through time using the 33 flights available from the ASO. For each grid cell, we analyzed 33 snow depth data points versus 33 fSCA data points and characterized the depletion curve with a bilinear regression between snow depth (dependent variable) and fSCA (independent variable). Theoretical considerations would motivate a curvilinear fit because fSCA can, by definition, only increase to a value of one, whereas maximum snow depth is unbounded (Liston 2004). However, we observed distinctly linear relationships between snow depth and fSCA for low fSCA and snow depth values in each grid cell, which is consistent with previous studies (Swenson & Lawrence 2012). As stated above, the objective of this work is not to parameterize a universal depletion curve but rather to evaluate the relationships between topographic characteristics and snow depletion dynamics. Hence, we quantify the depletion curve for each 500 m pixel based on a segmented bilinear regression, resulting in two slopes and a breakpoint value between the two linear segments.

We tested the statistical relationship for a significant breakpoint (p < 0.1) in the bilinear regression based on a Davies' Test as implemented in the R package ‘segmented’ (Muggeo 2003, 2020). The coordinate of the breakpoint represents the fSCA and snow depth at which the fitted slope of the depletion curve changes (Figure 2(a)–2(c)). Grid cells without a significant breakpoint were fit with an ordinary least squares linear regression and the slope tested for significance. Grid cells with no breakpoint but with significant slopes less than 1.6 were included in our analyses under the assumption that these slopes reflect a similar process as the first slope of the bilinear regression but did not accumulate snow beyond the theoretical breakpoint (Figure 2(d)). Visual inspection of the data indicated that significant linear slopes greater than 1.6 typically represented grid cells that did not fully deplete within the ASO observation period and the x-intercept for these linear models were often much greater than zero (Supplementary Figure S1).

Consequently, we examined the topographic controls on the first slope of the bilinear regression and slope values of the ordinary linear regression below values of 1.6, and we denote these as the ‘depletion slope’ – directly defined by the breakpoint coordinates and a y-intercept of zero. Significant y-intercept differences from zero (p < 0.05) are rare, and we consider this to be a function of vertical uncertainty in the measurements. The upper slope of the bilinear regression was found to be uncorrelated with topography. This is consistent with the concept that the upper portions of the curves describe dynamics of snow depths greater than the sub-grid topographic variability. Thus, they are not explicitly addressed herein.

We analyze 2,596 grid cells (500 m) for relationships between topographic variables and the depletion slope. Focusing on the depletion slope is particularly revealing with regard to the compounded result from various snow distribution processes. Fundamentally, steeper depletion slopes describe a relatively heterogeneous snow depth distribution where snowpacks covering rough terrain result in areas of relatively deep snow accumulations. As such, for an equivalent decrease in snow depth, the decline in fSCA for grid cells with these deep pockets of snow is relatively small compared with grid cells with shallower, more uniform snow distributions. Intuitively, these different pockets of snow are influenced by the underlying and surrounding topography.

We use regression tree analysis (Breiman 1984) to identify potential physiographic controls (independent variables) on the depletion slope (dependent variable). We consider a range of physiographic variables that have demonstrated physical or statistical relationships with snow depth previously under the assumption that the snow depth–fSCA relationship is determined both by the magnitude and variability of snow depth. We investigate three different scales of control including grid cell scale (500 m) physiography, inter-grid cell effects, and sub-grid variability within each grid cell (Table 1).

First, the pixel-scale physiographic variables represent potential physical process controls on the magnitude of snow depth: elevation (elev), slope (zness), aspect (northness; eastness), and vegetation height (vegheight). Second, the inter-pixel topographic variables, topographic position index (tpi) and relief, quantify the topographic variability within a 3 × 3 pixel window around each 500 m grid cell and also represent potential controls on the magnitude of snow depth. Third, variables representing sub-grid cell scale variability describe terrain roughness based on elevation, slope, and curvature within each 500 m grid cell and represent potential controls on the variability of snow depth. Further information on the pixel-scale and inter-pixel-scale variables are provided in Table 1 and the references therein. Further detail on the sub-grid cell variables is included below as these variables are less prevalent in the literature.

At the coarser scale, we quantify the standard deviation of elevation (stdelev), the standard deviation of slope (stdslope), and the standard deviation of maximum curvature (stdmaxcurv). Previous studies have shown the length scale at which snow distribution and depletion processes change to be between 6 and 40 m (Deems et al. 2006; Trujillo et al. 2007; Schirmer et al. 2011), so we use a 45 m resolution DEM in order to explicitly treat variation native to processes operative at, or greater than, this scale and the 3 m DEM to capture the processes at finer scales.

Regression trees have commonly been used to estimate the spatial distribution of snow (e.g. Elder et al. 1998; Balk & Elder 2000; Molotch et al. 2005). A regression tree can model nonlinear interactions between the topography and depletion slope. We grew regression trees until additional splits no longer improved the cross-validated r² by 0.001. We subsequently pruned the regression trees to an optimal complexity by selecting the number of splits corresponding to within 1 standard error of the minimum mean squared error (Breiman 1984; R Core Team 2015; Therneau et al. 2015).

Relationships between snow depth and fSCA are remarkably consistent in all 4 years despite significant interannual variability in weather patterns, i.e. points from different years plot along the same depletion curve in Figure 2. On average, the model for each grid cell explained 81% of the variance in snow depth (Figure 2(e)). Moreover, 84% of depletion slopes fall within the 95% confidence interval of the mean cross-validated depletion slope (where four depletion slopes were calculated with the data for each year iteratively removed). This suggests that the snow depth–fSCA relationships presented here are robust even in a year with an otherwise anomalous snowpack.

The range of the 2,596 depletion slopes analyzed with respect to topography is 0.20–3.06 with a mean of 0.73 (Figure 3(a)). The spatial pattern of grid cells that had statistically significant depletion slopes corresponds to higher elevation grid cells (mean 2,822 m with a range of 1,276–3,607 m) of the basin (Figure 3(a) and 3(b); see Figure 1 for elevation contours). At lower elevations, shallow snow depths and a lack of clear depletion dynamics – likely driven by the drought conditions during the study period – rendered the regressions statistically insignificant (white, Figure 3(a)). The highest depletion slopes (darkest colors, Figure 3(a)) roughly correspond with deeper snow and occur along the southern and northern basin divides. The lowest values (lightest colors, Figure 3(a)) are clustered in the southeast.

The temporally averaged standard deviations of snow depth (Figure 3(c)) exhibit a similar spatial pattern as the depletion slope (Figure 3(a)). A positive relationship between depletion slope and snow depth standard deviation exists across the basin with r² = 0.58; p < 0.01 (Figure 3(d)). This supports our interpretation of the depletion slope that grid cells with steeper depletion slopes have more heterogeneous snowpacks.

None of the topographic variables has a strong individual relationship with the depletion slope (Table 2).

The highest magnitude correlation with the depletion slope is displayed by the standard deviation of maximum terrain curvature (r = 0.27) followed by the Δh/Δw parameter (r = 0.26) and the standard deviation of terrain slope (r = 0.25). The prevalence of three different sub-pixel topographic variables indicates that sub-pixel terrain variability exerts substantial influence on the depletion slope. This finding agrees well with previous work, showing snow depth variability is well correlated with topographic variability, in part due to wind redistribution (Winstral et al. 2002; Erickson et al. 2005). All correlations are statistically significant except for eastness, possibly due to gradients in precipitation shadowing from north to south in the Tuolumne Basin.

The regression tree analysis explains 31% (95% CI: ±0.2%) of the variance in the depletion slope (Figure 4). Vegetation height explains the most variance with the depletion slope, with taller vegetation corresponding with smaller depletion slopes and shorter to no vegetation corresponding with larger depletion slopes. This implies that taller vegetation leads to more homogeneous snowpacks and shorter vegetation leads to more heterogeneous snowpacks, which agree with the previous work showing that areas with forests have more homogenous snow distribution due to reduced wind redistribution (Hiemstra et al. 2002; López-Moreno & Latron 2008). We interpret this split in the regression tree as indicative of differences in accumulation and depletion dynamics above and below the timberline. We note that northness, elev, and multiple sub-grid cell variables are used for subsequent splits both above and below the tree line. A key difference between the regression tree results and the single variable correlation presented in Table 2 is the nonlinearity inherent in the regression tree results. This is apparent with the appearance of elev, stdmaxcurv, and northness multiple times at different levels of the tree. They explain a higher variance than other variables for individual subsets of the data.

We observe that increased terrain roughness at sub-grid scales – represented by the stdmaxcurv, stdelev, and Δh/Δw variables – correlates with larger depletion slopes at multiple nodes in the regression tree. Lower sub-grid terrain variability is associated with lower depletion slope values, consistent with expectations of relatively homogenous snow distribution. This finding is intuitive, given that these areas have shallower terrain depressions in which to trap snow. In the latter stages of depletion, decreases in snow depth in these areas will yield relatively large decreases in SCA because large portions of bare ground will be exposed nearly simultaneously across the grid cell. We also tested sub-grid standard deviation of vegetation height as another measure of sub-grid variability that could potentially be important near the tree line, but found mean vegetation height to explain more variance.

Northness occurs twice in the regression tree, and in each case, larger northness values are associated with increased depletion slopes, i.e. north-facing grid cells exhibit more heterogeneous snow distribution which is congruent with Lehning et al. (2011). This may be a result of preferential snow accumulation on northeast-facing slopes due to preferential deposition of precipitation and snow drift (Lehning et al. 2008; Dadic et al. 2010). Moreover, flow separation across ridge-tops can produce cornices which can exhibit particularly heterogeneous snow depths (Anderton et al. 2004). We also considered increased small-scale terrain variability, e.g. boulder fields, as a confounding factor but found no correlation between Δh/Δw and northness.

Higher elevations are associated with steeper depletion slopes and therefore a more heterogeneous snowpack. Further work is warranted to explore the physical processes linking the depletion slope and elevation, but several possibilities exist. Increased snow heterogeneity at higher elevations can be produced by more prevalent wind redistribution effects, i.e. at higher elevations, increased snowfall and wind speeds and lower-density snow result in greater wind redistribution relative to lower elevations. In addition, ridgelines that exhibit cornices are typically at the higher elevations. Conversely, shallow depletion slopes would result from a less-distinct accumulation season, e.g. at the lowest elevations, where snowfall results in broad, shallow coverage but quickly melts away. In addition, we might also expect grid cells with high elevation and high terrain curvature to have steep depletion slopes as a consequence of avalanche redistribution which leaves narrow areas of high snow depth.

The relationship between snow depth and fSCA appears to be consistent across 4 years during the depletion period, suggesting that remotely sensed fSCA contains a large amount of information regarding snow distribution. Moreover, we show, for the first time, that repeat lidar scans of snow depth distribution are a novel data source for developing fSCA–snow depth relationships over large areas (i.e. >1,000 km²). The establishment of these relationships from explicit observations provides a framework for improving depletion curve parameterizations in hydrologic, land-surface, and climate models with the expansion of the ASO program to other mountain ranges and other efforts to explicitly map snow depth and fSCA (e.g. Nolan et al. 2015). Incorporating physically representative empirical factors that capture accumulation and ablation dynamics over multiple scales could result in improved spatial representation of snow-free and SCAs in models thus improving energy balance calculations.

The findings presented in this paper are conceptually consistent with a global classification system of SWE variability used for parameterizing probability distribution depletion curves (Sturm et al. 1995; Liston 2004). It is encouraging that depletion curves developed from high-resolution spatio-temporal data similarly describe more heterogeneous snowpacks in areas with greater topographic variability and greater wind speeds.

We see a significantly larger range in CV of snow depth than previously reported with the mean observed sub-grid cell CV (3 m spatial resolution snow depth within each 500 m grid cell) ranging from 2.4 to 34.0 across all dates. This is compared with field values of CV of snow depth for mountainous regions of 0.1–3 (Clark et al. 2011). We consider this inconsistency to be a function of scale. The 3 m spatial resolution snow depth dataset from ASO samples the complete range of terrain and snow depths compared with the studies reviewed in Clark et al. (2011), and thus captures snow depth extremes.

The results of our regression tree analyses are specific to the Tuolumne basin for two very dry years and two dry to average years. Future analyses of repeat lidar measurements across time and space are warranted to determine the spatio-temporal transferability of our results. In this regard, the regression tree analysis misses the extremes of the range of depletion slope values, which is not surprising given an r² value of 0.31. There is considerable variability in the observed depletion slopes (Figure 3(a)) that may have resulted from local processes that our predictor variables do not capture. This suggests that there were either too few grid cells with very high depletion slopes for a meaningful statistical relationship at the extremes of the observed parameters, or that a key variable is missing from our analysis. Furthermore, grid cells without a breakpoint but with linear slopes greater than 1.6 were excluded from the topographic analysis (see Methods and Supplementary Information for details); ASO flights did not extend late enough into the year to observe the lowest values of fSCA and depth for these grid cells.

Herein, we show that snow depletion curve shape is inherently linked to snow depth distribution. Hence, further analysis of remotely sensed fSCA depletion may improve understanding of the processes controlling snow distribution. In particular, the fSCA and snow depth components of the depletion slope may be analogous to the scale break or fractal dimension reported in previous studies that investigated fractal scaling in snow depth distributions (Deems et al. 2006, 2008; Trujillo et al. 2009; Schirmer & Lehning 2011). Further work aimed at identifying relationships between topography and the snow depth value corresponding to the depletion curve breakpoint could provide a means to parameterize depletion curve shape within hydrologic models. This could provide useful information not only about the depletion slope during ablation, but also about the capacity of the terrain to trap snow in concavities.

We also performed the same analysis with SWE from ASO to link to the previous literature but did not find a consistent bilinear relationship. We might expect differences in the shapes of depletion curves derived from SWE rather than snow depth because density strongly increases early in the ablation season, thus offsetting decreases in snow depth. Consistent density during melt as fSCA is decreasing may then contribute to a more dynamic relationship between SWE and fSCA. ASO SWE estimates use field-measurement constrained, model-based snow density information with density uncertainties of 12–30 kg m⁻³ (Painter et al. 2016). Future work evaluating SWE depletion with respect to topography is warranted. Given the extensive use of SWE depletion curves in hydrologic modeling, improved linkage between snow depth depletion curves and SWE depletion curves would also prove valuable.

We derived individual depletion curves for 500 m grid cells using a lidar-derived, high-resolution, spatio-temporal dataset of observed snow depth and analyzed the effect of physiography on these depletion curves. We found that relationships between snow depth and fSCA (i.e. depletion slope) were robust over the 4 years of study with a mean r² value of 0.81. Additionally, 84% of presented depletion slopes fall within the 95% confidence interval of the mean cross-validated depletion slope. These depletion slopes exhibited significant variability across the watershed. In this context, we show that a positive relationship with r² = 0.58; p < 0.01 between depletion slopes and snow depth variability exists across the basin. We also show that sub-pixel and pixel-scale terrain variables explain 31% of the spatial variability in the depletion slope. In particular, increased vegetation height and decreased sub-pixel terrain variability were associated with more homogeneous snowpacks and lower depletion slopes. These results illustrate that repeat, distributed snow depth measurements such as those from the NASA ASO can provide insights to the influence of topography on the evolution of snow distribution. Such understanding has important implications for developing parameterizations of snow cover depletion curves across physiographic gradients. Given that parameterizations of snow cover depletion underpin sub-grid representation of energy and water fluxes across a range of earth system models, the results and approach presented herein has potentially broad applicability.

Data are available from https://doi.org/10.17605/OSF.IO/A96BP. D.S. was supported by a NASA Earth and Space Science Fellowship (NNX14AL27H). N.P.M. was supported by NASA Grants NNX17AF50G and 80NSSC17K0071, and USDA grant 2012-67003-19802. T.H.P. acknowledges the NASA Terrestrial Hydrology Program, NASA Applied Sciences, and the California Department of Water Resources. Part of this work was performed at the Jet Propulsion Laboratory, California Institute of Technology under a contract with NASA. Special thanks to contributors to R, the tidyverse and raster packages, and GDAL (Bivand et al. 2015; GDAL/OGR contributors 2015; Hijmans 2015; R Core Team 2015; Wickham et al. 2019).

All relevant data are available from an online repository or repositories (https://doi.org/10.17605/OSF.IO/A96BP).