## Abstract

The annual consistency of spatial patterns of snow accumulation and melt suggests that the evolution of these patterns, known as depletion curves, is useful for estimating basin water content and runoff prediction. Theoretical snow cover depletion curves are used in models to parameterize fractional snow-covered area (fSCA) based on modeled estimates of snow accumulation and snowmelt. Directly measuring the spatio-temporal snow distribution, characterization of depletion curves, and understanding how they vary across mountainous landscapes was not possible until the recent U.S. National Aeronautics and Space Administration (NASA) Airborne Snow Observatory (ASO). Herein, for the first time, high-resolution spatio-temporal snow depth information from the ASO is used to derive observation-based snow cover depletion curves across physiographic gradients by estimating the slope of the fSCA–snow depth relationship (i.e. depletion slope). The depletion slope reveals important insights into snow processes as it is strongly related to snow depth variability (*r*^{2} = 0.58). Regression tree analysis between observed depletion slopes and physiography, particularly vegetation height and terrain roughness, displays clear nonlinear dynamics and explains 31% of the variance in depletion slope. This unique observation-based analysis of snow cover depletion curves has implications for energy and water flux calculations across many earth system models.

## HIGHLIGHTS

Relationships between snow depth and fSCA (i.e. depletion slope) were robust over the 4 years of study.

Significant spatial variability in depletion slope is well correlated with snow depth variability.

Increased vegetation height and decreased terrain roughness were associated with more homogeneous snowpacks and lower depletion slopes.

## INTRODUCTION

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?

## METHODS

### Site description

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^{2} 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.

### Data sources

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).

### Deriving depletion curves

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.

### Relating depletion slope to topography

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.

Variable name . | Short name . | Resolution . | Definition/Units . | Source . |
---|---|---|---|---|

Grid cell scale physiography | ||||

Elevation | elev | 495 m | Elevation from the summer flight; meters | Elder et al. (1991, 1998) |

Slope | zness | 495 m | Sin(slope); ranges 0–1; dimensionless | Balk & Elder (2000); Erxleben et al. (2002); Fassnacht et al. (2003) |

Northness | northness | 495 m | Cos(aspect); ranges 0–1; dimensionless | Balk & Elder (2000); Erxleben et al. (2002); Fassnacht et al. (2003) |

Eastness | eastness | 495 m | Sin(aspect); ranges 0–1; dimensionless | Balk & Elder (2000); Erxleben et al. (2002); Fassnacht et al. (2003) |

Vegetation height | vegheight | 495 m | Mean vegetation height measured by ASO from the summer flight; meters | Deems et al. (2006); Trujillo et al. (2007); Painter et al. (2016) |

Inter-grid cell variability | ||||

Topographic position index | tpi | 495 m from a 3 × 3 pixel window around each pixel | Elevation difference of a pixel from the mean of the surrounding pixels; meters | Revuelto et al. (2014); GDAL (2015) |

Relief | relief | 495 m from a 3 × 3 pixel window around each pixel | Maximum elevation difference between any two of the surrounding pixels; meters | Bookhagen & Burbank (2010); GDAL (2015) |

Sub-grid variability | ||||

Standard deviation of elevation | stdelev | 495 m from 45 m | Standard deviation of elevation from 45 m DEM; meters | Marchand & Killingveit (2005) |

Standard deviation of slope | stdslope | 495 m from 45 m | Standard deviation of slope from 45 m DEM; shown to detect changes in slope at multiple scales; radians | Grohmann et al. (2011); Fassnacht et al. (2016) |

Standard deviation of maximum curvature | stdmaxcurv | 495 m from 45 m | Standard deviation of maximum curvature from 45 m DEM; maximum curvature from SAGA GIS; dimensionless | Marchand & Killingveit (2005); López-Moreno et al. (2014); Conrad et al. (2015) |

Terrain variability | Δh/Δw | 495 m from 3 m | Terrain variability metric defined in Equation (1) in Helbig et al. (2015); computed from the 3 m DEM; dimensionless | Helbig et al. (2015) |

Variable name . | Short name . | Resolution . | Definition/Units . | Source . |
---|---|---|---|---|

Grid cell scale physiography | ||||

Elevation | elev | 495 m | Elevation from the summer flight; meters | Elder et al. (1991, 1998) |

Slope | zness | 495 m | Sin(slope); ranges 0–1; dimensionless | Balk & Elder (2000); Erxleben et al. (2002); Fassnacht et al. (2003) |

Northness | northness | 495 m | Cos(aspect); ranges 0–1; dimensionless | Balk & Elder (2000); Erxleben et al. (2002); Fassnacht et al. (2003) |

Eastness | eastness | 495 m | Sin(aspect); ranges 0–1; dimensionless | Balk & Elder (2000); Erxleben et al. (2002); Fassnacht et al. (2003) |

Vegetation height | vegheight | 495 m | Mean vegetation height measured by ASO from the summer flight; meters | Deems et al. (2006); Trujillo et al. (2007); Painter et al. (2016) |

Inter-grid cell variability | ||||

Topographic position index | tpi | 495 m from a 3 × 3 pixel window around each pixel | Elevation difference of a pixel from the mean of the surrounding pixels; meters | Revuelto et al. (2014); GDAL (2015) |

Relief | relief | 495 m from a 3 × 3 pixel window around each pixel | Maximum elevation difference between any two of the surrounding pixels; meters | Bookhagen & Burbank (2010); GDAL (2015) |

Sub-grid variability | ||||

Standard deviation of elevation | stdelev | 495 m from 45 m | Standard deviation of elevation from 45 m DEM; meters | Marchand & Killingveit (2005) |

Standard deviation of slope | stdslope | 495 m from 45 m | Standard deviation of slope from 45 m DEM; shown to detect changes in slope at multiple scales; radians | Grohmann et al. (2011); Fassnacht et al. (2016) |

Standard deviation of maximum curvature | stdmaxcurv | 495 m from 45 m | Standard deviation of maximum curvature from 45 m DEM; maximum curvature from SAGA GIS; dimensionless | Marchand & Killingveit (2005); López-Moreno et al. (2014); Conrad et al. (2015) |

Terrain variability | Δh/Δw | 495 m from 3 m | Terrain variability metric defined in Equation (1) in Helbig et al. (2015); computed from the 3 m DEM; dimensionless | Helbig et al. (2015) |

Variables were chosen based on their existing use in the literature; sources listed are by no means exhaustive.

*h*/Δ

*w*’ parameter, which represents the ratio of the typical height of a terrain feature to the typical width of a terrain feature. Larger Δ

*h*/Δ

*w*values indicate a greater relative protrusion of the terrain. We calculated Δ

*h*/Δ

*w*for each 500 m grid cell from the 3 m ASO DEM after the following equation (Helbig

*et al.*2015):where ∂

*and ∂*

_{x}z*are the first partial derivatives of elevation, i.e. the orthogonal slope components (Neteler*

_{y}z*et al.*2012).

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*^{2} 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).

## RESULTS AND DISCUSSION

### Depletion curves

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 slope*s 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 slope*s 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*^{2} = 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.

### Topographic controls on depletion curves

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

Topographic variable . | Kendall rank correlation . | p-value
. |
---|---|---|

stdmaxcurv | 0.27 | < 0.01 |

Δh/Δw | 0.26 | < 0.01 |

stdslope | 0.25 | < 0.01 |

vegheight | − 0.23 | < 0.01 |

stdelev | 0.21 | < 0.01 |

elev | 0.20 | < 0.01 |

northness | 0.19 | < 0.01 |

relief | 0.13 | < 0.01 |

tpi | 0.13 | < 0.01 |

zness | 0.09 | < 0.01 |

eastness | 0.01 | 0.41 |

Topographic variable . | Kendall rank correlation . | p-value
. |
---|---|---|

stdmaxcurv | 0.27 | < 0.01 |

Δh/Δw | 0.26 | < 0.01 |

stdslope | 0.25 | < 0.01 |

vegheight | − 0.23 | < 0.01 |

stdelev | 0.21 | < 0.01 |

elev | 0.20 | < 0.01 |

northness | 0.19 | < 0.01 |

relief | 0.13 | < 0.01 |

tpi | 0.13 | < 0.01 |

zness | 0.09 | < 0.01 |

eastness | 0.01 | 0.41 |

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.

### Implications

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^{2}). 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*^{2} 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^{−3} (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.

## CONCLUSIONS

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*^{2} 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*^{2} = 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.

## ACKNOWLEDGEMENTS

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).

## DATA AVAILABILITY STATEMENT

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

## REFERENCES

*segmented: Regression Models with Break-Points/Change-Points Estimation*. http://cran.r-project.org/package=segmented