What drives basin scale spatial variability of snow water equivalent during two extreme years ?

Introduction Conclusions References


Introduction
A majority of earth's moving freshwater originates in snow dominated mountainous areas (Viviroli et al., 2003), with 60 to 75 % of annual streamflow in the Rocky Mountain region of the western United States originating from snowmelt (Doesken and Judson, 1996).A comprehensive understanding of the distribution of the seasonal mountain Introduction

Conclusions References
Tables Figures

Back Close
Full snowpack and estimation of its snow water equivalent (SWE) is essential to improve hydrologic models used for forecasting water availability.Additionally, the recent shift towards earlier snowmelt in regions of the western US (e.g.Stewart, 2009;Clow, 2010) necessitates a more accurate accounting for future water resources planning.Mountainous landscapes have complex topography as well as strong and highly variable climatic gradients yielding spatial and temporal (seasonal and interannual) variability in snowpack properties.Determining the meteorology and related feedbacks that drive hydrologic processes in these areas is challenging in such complex terrain and requires spatial scaling (Bales et al., 2006).Often the resolution of available SWE measurements is much larger than the scale needed to characterize the correlation length of its spatial variability (Blöschl, 1999).Across the western US, the Natural Resource Conservation Service (NRCS) SNOwpack TELemetry (SNOTEL) and snow course network provide operational snowpack measurements of snow depth and SWE and thus calculated average density at a daily and monthly time step, respectively.NRCS operational stations were established to measure the snowpack for water supply forecasts, yet, they have been shown to represent SWE only as a point location rather than surrounding area (Molotch and Bales, 2005).Nonetheless, SNOTEL and snow course sites are the most widely available and utilized ground based measurements of SWE and relied upon heavily for estimating basin scale snow distribution.
Research on the spatial distribution of snow has emphasized the statistical relation between snow properties and terrain characteristics, the latter as a surrogate for the driving meteorology.These studies have used SNOTEL data to interpolate SWE over large basins (e.g.Fassnacht et al., 2003), as well as snowpack field measurements over small catchments (e.g.Elder et al., 1991).However, few studies have analysed snow's spatial variability at the basin scale using both operational and field measurements.Operational measurements can provide regional knowledge on the spatial distribution of snow (e.g.Fassnacht et al., 2003), yet cannot accurately characterize the spatial variability of the snowpack at the basin scale (Bales et al., 2006).It has been Introduction

Conclusions References
Tables Figures

Back Close
Full recommended that future research should focus on more accurate estimation of SWE at the basin (100s to 1000s km 2 ) and regional (10 000s to 100 000s km 2 ) scale to effectively assess and manage mountain water resources (Viviroli et al., 2011), thus there is need to supplement operational data with additional field-based snowpack measurements at this scale of interest to evaluate the spatial variability of SWE and provide additional ground truth measurements within the scale extent of remote sensing observations.At the basin scale, an approach to reducing the sampling effort needed for more measurements is to use snow depth as a surrogate for SWE by developing a model for snow density, as manual snow density measurements require more time and effort than snow depth measurements.Recent studies have attempted to characterize the spatiotemporal characteristics of snow density (e.g.Mizukami and Perica, 2008;Fassnacht et al., 2010), or to develop reliable methods for modeling snow density and thus estimating SWE from snow depth measurements (e.g.Jonas et al., 2009;Sturm et al., 2010).The objectives of this research were: (1) to develop a basin scale snow density model that can be used to estimate SWE from snow depth measurements; this is a different domain and scale than used in previous studies, and (2) to combine operational SNO-TEL and snow course measurements, as suggested by Dressler et al. (2006), with supporting field-based snowpack measurements to evaluate what is driving variability of the snowpack at the basin scale.

Study area and datasets
This study was conducted in the Cache la Poudre basin located in northern Colorado and a small portion of southeastern Wyoming (Fig. 1).The upper portion of the basin that contributes to the canyon mouth (gaged by Colorado Division of Water Resources gaging station Cache la Poudre River at Canyon Mouth) has an area of 2729 km 2 and ranges in elevation from 1591 to 4125 m.We focus on this portion of the basin because it is responsible for the majority of hydrologic input to the river system.Spruce-fir (Picea Introduction

Conclusions References
Tables Figures

Back Close
Full engelmannii and Abies lasiocarpa), lodgepole pine (Pinus contorta), and ponderosa pine (Pinus ponderosa) forests cover a majority of this area, with the alpine community located at the highest elevations and the mountain shrub community located at the lowest elevations.Snow is the dominant form of precipitation within the basin, and the hydrograph peak is driven by snowmelt generally occurring in late May to June.The majority of winter moisture moves into this region by pacific frontal storm tracks from the west, southwest, or northwest, however, systems moving north from the Gulf of Mexico can also bring substantial snowfall to the Front Range of Colorado (Barry, 2008).
The NRCS operational stations located within the study area and in a 15 km buffer around the basin were analysed (Fig. 1), yielding a total of 10 SNOTEL stations and 17 snow courses.Deadman Hill and Joe Wright, the two long-term SNOTEL stations located within the Cache la Poudre basin, have a mean (1980 to 2012) peak SWE of 538 mm and 690 mm, respectively (Fig. 2).The lowest snow year recorded was 2002 at Deadman Hill and 2012 at Joe Wright, while the maximum snow year was 2011 at both SNOTEL stations.Despite the similar elevation of the two stations, historically Joe Wright has a greater accumulation of snow than Deadman Hill.
Field snow surveys were conducted on and about 1 April 2011 and 2012 within the study area.At each field sampling location, snow density (ρ s ) and/or snow depth (d s ) measurements were taken and Universal Transverse Mercator (UTM) geographic coordinates were recorded using a Global Positioning System (GPS).Eleven measurement points of snow depth (using a snow depth probe to the near cm of depth) were averaged across a one-meter interval in one of the four cardinal directions to account for the small scale spatial variability at a location (e.g.López-Moreno et al., 2011).Snow density is a conservative variable that varies less spatially than depth (Logan, 1973;Fassnacht et al., 2010), thus, fewer snowpack density measurements were made across the study area than depth.Three methods of measuring snow density were used at each site.A cylindrical metal can with a diameter of 15.3 cm was used to measure snow density if the snowpack was less than 50 cm.A cylindrical plastic snow sampling tube with a diameter of 6.6 cm was used to measure snow density for snowpacks greater than 50 cm Figures

Back Close
Full and less than 150 cm.Additionally, a Federal Sampler (diameter of 3.77 cm) was used to measure the snow density for snowpacks greater than 150 cm, but it was also used at some locations shallower than 150 cm.Each of the field-based surveys included multiple transects, with each transect consisting of a number of snowpack sampling locations with a systematic sampling spacing of approximately 500 m.A total of 42 field sampling locations were monitored on and about 1 April 2011 and 121 field sampling locations on and about 1 April 2012.The location of snow survey transects were selected based on accessibility as well as representation of snow producing regions within the study area.The high elevation areas located around the Colorado State University Pingree Park Campus, Cameron Pass, and Deadman Hill were the focus within the Cache la Poudre basin (Fig. 1).The 2011 field-based snow survey was completed over the span of three days (31 March through 2 April 2011), while the 2012 survey was completed over four days (29 March through 1 April 2012).Small amounts of precipitation was recorded at SNO-TEL stations within the study area during the 2011 and 2012 survey time period, however the majority of change to the snowpack during these periods were due to melt, compaction, and/or metamorphism.Changes in snow depth were accounted for using daily SNOTEL snow depth measurements to standardize the field-based snow depth measurements to a single date for each survey.The average change in snow depth among SNOTEL stations was added to our field-based snow depth measurements outside of the standardized date to adjust for the change in snow depth over that period.Snow depth measurements from the 2011 survey were standardized to 2 April, while 2012 measurements were standardized to 31 March.Introduction

Conclusions References
Tables Figures

Back Close
Full

Snow density model
Since the range of variability of snow density is more conservative than snow depth and SWE (e.g.Logan, 1973;Fassnacht et al., 2010) 1) were evaluated.These snow courses range in elevation from 2408 m to 3261 m and are (generally) measured on or about the first of the month from January through June each year.For the analysis, snow density values greater than 600 kg m −3 and less than 50 kg m −3 were omitted.Additionally, due to the limited precision and possibly the lack of accuracy for snow density measurements in shallow snowpacks, data for snow depth less than 0.13 m and/or SWE less than 50 mm were also omitted.This selection of data resulted in 3262 data records of snow depth, snow density, and SWE, with 10.3 % of the original data being removed.
The seasonal variability of snow density is largely dictated by time of year, while inter-annual variability is minimal (Mizukami and Perica, 2008).Snow density tends to increase gradually throughout the snow season due to crystal metamorphism, settling, and compaction.Therefore, snow density tends to increase with the day of year (Mizukami and Perica, 2008) and with increasing snow depth (Pomeroy and Gray, 1995).Topography and tree canopy also impact snow densification, as they can be surrogates for solar radiation.However, snow courses are often located in flat open Introduction

Conclusions References
Tables Figures

Back Close
Full A multiple linear regression model was developed to predict snow density considering snow depth (d s ), Julian day (DOY), elevation (z), UTM Easting (UTM e ), and UTM Northing (UTM n ) as independent variables.The statistical software R (R Development Core Team, 2012) was used for all statistical analyses.The final independent variables included in the multiple linear regression model were selected based on an all-subsets regression procedure (Berk, 1978), which assesses a criterion statistic for every possible combination of independent variables.Mallows' C p (Mallows, 1973), which assesses the fit of a regression model and increases a penalty term as the number of predictor variables increases, was used as a criterion for the all-subsets regression.Candidate models that showed the best Mallows' C p values were then evaluated using the Akaike information criterion (AIC) statistic (Akaike, 1974), which is also a measure of the relative goodness of fit of the statistical model that introduces a penalty for increasing the number of model parameters.The variance inflation factor (VIF) was used to quantify the severity of multicollinearity between independent variables.A VIF score greater than 4 may suggest multicollinearity between variables (Kutner et al., 2005).Model diagnostics were evaluated using residual plots to check the model assumptions of normality, linearity, and homoscedasticity and were used to determine if variable transformations were necessary.Final model selection was based on the results of criterion statistics and model diagnostics (Kutner et al., 2005).
The multiple regression model provides an estimate of snow density for each snow depth measurement and their product yields an estimate of SWE.To assess the accuracy of the snow density model, several methods of model evaluation were performed.Calibration was performed by comparing modeled snow density as well as calculated SWE with observed values from the original dataset; explained variance was computed.Model validation with two sets of independent data was also completed to test model transferability to predict independent data.The two independent datasets included field-based measurements from the 2011 and 2012 snow seasons (n = 84), as Introduction

Conclusions References
Tables Figures

Back Close
Full well as historic first of the month SNOTEL measurements (n = 121) at sites that are not co-located with a snow course.Additionally, a 10-fold cross validation procedure, which runs 10 iterations of removing a random selection of the dataset and fitting the regression to the remainder of the data, was used to compare modeled values to the observed values removed for each iteration.Performance of the final snow density model was determined from the residuals of both observed snow density as well as calculated SWE through the calculation of the Nash-Sutcliffe Coefficient of Efficiency (NSCE) and Root Mean Squared Error (RMSE) performance statistics.

Basin scale SWE variability
Topographic variables that are thought to potentially drive the spatial distribution of snow at the scale of interest were derived from a 30 m resolution digital elevation model (DEM) of the study area.The DEM was downloaded from the USGS National Elevation Dataset (NED) (http://ned.usgs.gov).Additionally, canopy density was obtained from the National Land Cover Database (NLCD 2001) (http://www.mrlc.gov).A value of the derived terrain and canopy variables (spatial data grids) was extracted for each sampling location based on its corresponding 30 m DEM pixel.A description of the derivation and importance of each of the spatial data grids is provided below.
Location within the study area is represented by UTM Zone 13N Easting and Northing coordinates for each field-based and operational sampling location.A 30 m resolution spatial data grid of UTM Easting and Northing was created for the study area in ArcGIS (ESRI, 2011) by assigning the mean UTM value to each pixel.Spatially continuous coordinates of UTM Easting and Northing can be correlated with the distribution of snow in various ways that depend on site location and scale.Previous studies have used distance to a mountain barrier and distance to ocean or source of moisture (e.g.Fassnacht et al., 2003;López-Moreno and Nogués-Bravo, 2006), which can also be represented by UTM Easting for the study site due to its geographic orientation.Furthermore, given the scale of the study, UTM Easting and Northing represent different Introduction

Conclusions References
Tables Figures

Back Close
Full regions within the study area that are thought to display different patterns of snow accumulation and ablation due to the variability of meteorology and storm tracks.Elevation was extracted for each sampling location directly from the 30 m DEM.Snow accumulation has long been shown to be a function of elevation (e.g.Washichak and McAndrew, 1967;Dingman, 1981) due to orographic precipitation patterns and the effect of air temperature (Doesken and Judson, 1996).
Slope was derived from the 30 m DEM using the Spatial Analyst tools within ArcGIS to provide an output spatial data grid with a value of slope (in degrees) for each pixel.The degree of slope impacts the stability of the snowpack (influencing snow accumulation and redistribution) and input of solar radiation (influencing melt) (Anderton et al., 2004).Previous studies have successfully used slope angle as an explanatory variable within statistical models describing the distribution of snow (e.g.Erxleben et al., 2002;Winstral et al., 2002).
Aspect (in degrees) was also derived from the 30 m DEM using the Spatial Analyst tools within ArcGIS.Aspect can be problematic as an independent variable due to its continuous range of 0 to 360 • , thus normalizing this variable is necessary.Degrees of northness and eastness were calculated to normalize the aspect variable (Fassnacht et al., 2001(Fassnacht et al., , 2012)).Degree of northness is the product of the cosine of aspect and the sine of slope (Molotch et al., 2005), while degree of eastness is the product of the sine of aspect and the sine of slope.Exposure of slope aspect controls solar radiation input, which influences snowpack stability, densification, and ablation (McClung and Schaerer, 2006).
Solar radiation was derived using the Area Solar Radiation tool in ArcGIS, which calculates incoming solar radiation across a DEM surface for a specified time interval.Given the latitude of the study area, the cumulative clear sky solar radiation (in WH m −2 ) from 15 November through 30 March was calculated for each pixel.Cumulative incoming solar radiation is calculated based on solar zenith angle and terrain shading, and does not consider the influence of forest canopy.Previous studies have successfully used solar radiation spatial data grids derived by similar methods within Introduction

Conclusions References
Tables Figures
Terrain curvature was derived from the 30 m DEM using the Spatial Analyst tools within ArcGIS to provide an output spatial data grid with a value of curvature for each pixel.Curvature is defined as the second derivative of the surface (Kimerling et al., 2011).This variable represents the local relief of terrain (i.e.concavity or convexity) in the direction of maximum slope, which, in terms of snow accumulation, primarily accounts for wind drifting from high exposure areas with steep slopes to low lying gullies (Blöschl et al., 1991;Lapen and Martz, 1996).
Maximum upwind slope (Winstral et al., 2002) is a terrain-based variable that has been shown to account for redistribution of snow by wind, which is especially important in alpine areas.However, this variable requires the knowledge of predominant wind direction to account for upwind terrain features, which is not measured across a basin scale, requiring a modeling approach (e.g.Liston and Sturm, 1998), thus it was not used in this study.
Canopy density is derived from Landsat Enhanced Thematic Mapper+ (ETM+) circa 2001 satellite data and DEM derivatives (Homer et al., 2007).The canopy density spatial data grid provides an estimated percentage of canopy cover for each pixel at a 30 m resolution.Canopy density can influence how snow is distributed across space as it is directly related to the amount of snow that is intercepted in the tree canopy.Snow sublimation from snow intercepted within the forest canopy is a major component of the overall water balance (Pomeroy and Gray, 1995).
Multiple linear regression was used to model 2 April 2011 and 31 March 2012 SWE based on its relation with independent physiographic variables identified above.A detailed description of the multiple linear regression methodology is provided above.Multiple linear regression models were developed using both field and operational snowpack measurements and also operational measurements only.At this scale of interest, operational data are commonly the only snowpack data available, thus it was useful to compare the results from operational data only to those results obtained from us-Introduction

Conclusions References
Tables Figures

Back Close
Full ing operational data and additional field-based measurements.The following notation will be used in this study: model O+F will refer to the multiple regression model using both field and operational snow measurements and model O will refer to the multiple regression model using only operational snow measurements.A total of four regression models were developed: model O+F11 (field and operational data from 2011), model O11 (operational data from 2011), model O+F12 (field and operational data from 2012), and model O12 (operational data from 2012).

Snow density model
The pairwise relations between snow depth, snow density, and SWE from the historic snow course records are presented in Fig. 3.A strong correlation exists between snow depth and SWE, which is best fit as a power function (Fig. 3a).There is considerable scatter about the linear fit for snow density versus snow depth (Fig. 3b), which suggests that additional variables should be included to describe the variability of snow density.Snowpack relations shown here are similar to those found in previous studies (e.g.Jonas et al., 2009;Sturm et al., 2010).The mean snow density from the snow course dataset is 287 kg m −3 with a standard deviation of 64.8 kg m −3 .SWE and snow depth have a greater standard deviation (178 mm, 0.46 m, respectively) compared to their mean (275 mm, 0.92 m, respectively) than snow density.Snow density is most highly correlated with Julian day, and also shows a strong positive correlation with snow depth and negative correlation with UTM Easting (Table 1).
The final snow density model takes the following form where ρ s is snow density, DOY is Julian day, d s is snow depth, z is elevation, and UTM e is UTM Easting.The variance inflation factor (VIF) is less than 2.2 for each variable within the final model, suggesting that multicollinearity between independent variables is not observed.The residuals of the regression model are normally distributed and do not violate the underlying assumptions of the regression (normality, linearity, homoscedasticity), thus no data transformations were necessary.
The calibrated model underestimated more dense snowpacks and overestimated less dense snowpacks, while calculated SWE showed generally unbiased residuals that tended to slightly increase with increasing observed SWE (Fig. 4).Performance statistics calculated from the residuals of calibration with the original dataset showed that predicted snow density explained 51 % of the total variance in the data with a RMSE of 45.3 kg m −3 , yet, calculated SWE was able to explain 94 % of the variance in the data and had a RMSE of 44 mm (Table 1).
Various methods of model evaluation were performed to test the utility of the regression model, including 10-fold cross validation, that all showed similar trends (Fig. 4) and comparable error estimates (Table 1) to model calibration.As expected, a minor increase in error estimation was observed for model validation with independent data, yet the minimal increase in error shows that the regression model should be transferable to independent data within the bounds of the original dataset.Thus, we used the snow density model to calculate SWE for WY 2011 and WY 2012 field-based snow depth measurements.

Basin scale SWE variability
A total of 51 and 127 snowpack measurements (both field-based and operational) were analysed from the 2 April 2011 (WY 2011) and 31 March 2012 (WY 2012) snow surveys, respectively (Fig. 1).The mean SWE and snow depth from WY 2011 were greater than WY 2012, yet the mean snow density and standard deviation of snow density was shown to be consistent among both years (Table 2).The WY 2011 was the maximum snow year on record within the study area, while WY 2012 was one of the least snowy Introduction

Conclusions References
Tables Figures

Back Close
Full years on record (Fig. 2); thus WY 2011 snowpack measurements were shown to have a higher mean SWE and snow depth, but also had a greater range of variability than that of WY 2012 (Table 2).From the average SWE among SNOTEL stations within the study area, the 1 April snow survey occurred before peak SWE in 2011, however, occurred after peak SWE in 2012.Analysis of the 1 April snowpack from these two water years allows for the comparison between the two extreme snow years (maximum and minimum) as well as between two different stages of the niveograph (during accumulation and melt).
Terrain and canopy variables derived within GIS at each of the snowpack measurement locations have similar averages when compared to the 50 % Snow Cover Index (SCI) (Richer et al., 2013) (Fig. 1), for both 2011 and 2012.Histograms of relative frequency (Fig. 5) show that the distribution of terrain and canopy variables sampled in 2011 and 2012 is similar to the 50 % SCI area distribution of these variables, suggesting that the snowpack measurement locations sampled during WY 2011 and WY 2012 are representative of the variability of physiography among the entire study area.
The range of values of terrain and canopy variables observed at operational stations tended to be smaller than the field-based station ranges (Fig. 5), which also suggests, the combination of field-based and operational measurements are more representative of the basin than the operational measurements alone.A formal Kolmogorov-Smirnov test (K-S test) for equality of distributions between a random sample (n = 244) of the continuous terrain and canopy variables within the 50 % SCI area of the basin versus the variables associated with the WY 2011 and WY 2012 measurement locations was completed.The K-S test shows that during both years the difference between the two samples for curvature, eastness, and canopy density is not significant enough (95 % significant level) to say they have a different distribution.However, a significant difference between the distributions of elevation, slope, northness, and solar radiation was observed for both years.The difference in elevation is obvious since field data are located more at higher elevations than the entire domain (Fig. 5a), and the operational data tend to be located in a small elevation zone (Fassnacht et al., 2012).Northness Introduction

Conclusions References
Tables Figures

Back Close
Full is highly correlated to solar radiation, and both are related to slope so the significance difference for each of these variables is partly based on their correlation.For avalanche safety purposes, manual measurements are usually on slopes less than 35 • , so steeper slopes can be underrepresented.Snowpack variables were shown to have a strong correlation with each other, with SWE and snow depth showing the strongest relation (consistent with the historic snow course dataset), while also showing to be highly correlated with elevation (Table 3).Bivariate relations showed SWE increased with increasing elevation, with the steepness of this slope being greater in WY 2011 than 2012.The strength of the correlation between SWE and elevation for WY 2011 (r = 0.75) and WY 2012 (r = 0.67) suggests that elevation is the most important physiographic variable for driving the distribution of SWE across the study domain, which is consistent with previous findings from studies evaluating SWE at the basin scale (e.g.Fassnacht et al., 2003;Jost et al., 2007;Harshburger et al., 2010).As UTM Northing increases, SWE decreases in WY 2011, suggesting northern regions of the study area receive less snow than southern regions (as suggested by J. Meiman, personal communication, 2010), yet this trend was not apparent in the low snow year of 2012.The apparent greater accumulation of snow in southern regions of the study area could be related to an upwind elevation gradient, with high peaks of Rocky Mountain National Park located in the southern portion of the study area, or due to the possibility of a dominant storm track that preferentially precipitates in southern regions before moving northward.SWE also decreased with increasing UTM Easting, which corresponds with both the effect of orographic precipitation within the study area (the continental divide is located on the western border of the study area), and also lower elevation regions receiving less snow than higher elevation regions.The other physiographic variables that are known to influence snow accumulation (e.g.forest canopy, aspect, and slope) did not exhibit a strong bivariate correlations with SWE; however, they may still be important in explaining variability of the datasets once the trends of elevation and UTM coordinates have been removed.

TCD Introduction Conclusions References
Tables Figures

Back Close
Full Multiple linear regression was used to model SWE for 2 April 2011 and 31 March 2012 with the field-based and operational snowpack dataset (model O+F ) and the operational snowpack dataset only (model O ) (Fig. 6).The final independent variables used within each model, beta coefficient values, and a summary of model calibration statistics is provided in Table 4 and Fig. 6.To satisfy the assumptions of the regression model and improve overall model performance, a square root transformation was made to SWE (dependent variable) and slope (independent variable) for model O+F11 , which explains 84 % of the total variance with an RMSE of 95.9 mm.Model O+F12 (no data transformations) explained 50 % of the total variance and showed an RMSE of 76.2 mm.The WY 2011 operational model (model O11 ) explains 90 % of the total variance of the data with an RMSE of 72.7 mm and includes a square root transformation of SWE.Lastly, model O12 explains 83 % of the total variance with a RMSE of 46.1 mm and includes a natural log transformation of slope.The VIF is less than 1.5 for each variable within all four of the multiple regression models, suggesting that multicollinearity between independent variables is not observed.Also, the residuals of each regression model do not violate the underlying assumptions of the regression (normality, linearity, homoscedasticity).
A comparison of the standardized error estimation between WY 2011 and WY 2012 models shows that the model O+F11 has a lower standardized typical magnitude of error (standardized RMSE) than model O+F12 , and describes more of the variance in the data (R 2 ) (Table 4).Similarly, model O11 has a lower standardized RMSE and greater R 2 value than model O12 , but the difference between these two models is less (Table 4).The difference among these performance statistics can partially be explained by the nature of each snow year (WY 2011 was the maximum snow year and WY 2012 was amongst the lowest) and sampling scheme.).The reduced model included UTM Easting, UTM Northing, elevation, eastness, and northness as independent variables and explained 67 % of the total variance with a standardized RMSE of 32 % (Fig. 6).These results show an improvement from the full model (model O+F12 ), suggesting that fewer data points may be increasing the model's ability to describe the variance of the data.Also, the reduced model showed a lower R 2 value than model O+F11 which suggests that the model performs better for the 2011 snow year due to the greater range of observed variability in the data.

Discussion
The snow density model developed across the study area performed relatively well in modeling SWE from independent snow depth measurements.Predicted SWE RMSE ranged from 44 mm (calibration data) to 66 mm (independent field validation data).Only 0.26 % (n = 11) of the validation data showed a residual value outside of one standard deviation of SWE from the original dataset (178 mm).Additionally, 80 % of all residual values (n = 2768) fell within ±50 mm.The variance of the model residuals were on average within 12.8 % of the observed values.Within site variability of SWE has been conservatively estimated to be 15 to 25 % (Jonas et al., 2009), which suggests that the error observed from the model is within the natural range of SWE variability at a site (Fassnacht et al., 2008).The small range of error suggests that estimating SWE from model developed here is its ability to estimate SWE from the most easily measured variable snow depth.Across basin scales, efforts are being made to estimate snow depth using satellite data, such as ICESat (Jasinski et al., 2012).This can improve snowpack estimates across varying domains of interest.This method is especially useful for fieldbased snow surveys at the basin scale, in which many snowpack measurements are required, and the assumption of a constant snow density across the study area is not valid (Lopéz-Moreno et al., 2013).
The snow density model is simple to develop and implement and an effective tool for obtaining estimates of SWE from snow depth measurements across basin scale domains.The model is however constricted to its spatial domain, range of physiographic inputs, as well as temporal coverage, thus, may not be applicable to areas outside of the study area, for elevations that are lower than 2408 m or higher than 3261 m, or for snow depths shallower than 0.20 m or deeper than 2.52 m.The snow course data were collected on or about the 1st of the month from January through June, and thus the model may be less suitable for mid-month days, and may not be useful before 1 January or after 1 June.
Similar snow density models have been developed from historic data for different domains.Jonas et al. (2009) developed a set of regression equations to model snow density using snow depth, day of year, elevation, and region for the Swiss Alps, while Sturm et al. ( 2010) employed a statistical method based on Bayesian analysis for the United States, Canada, and Switzerland using snow depth, day of year, and climate class.These previous studies and our research show that snow density is a conservative variable that varies spatially much less than snow depth and SWE.The previous studies used spatial domains that are orders of magnitude larger than what has been presented here, with the current data being at a finer resolution.While there are differences in the modeled scale, favorable results have been observed in each approach, suggesting this method is applicable for basin wide, regional, and global scales.
Despite WY 2011 being a maximum snow year and WY 2012 being a minimum snow year, the variables driving each SWE regression were similar and included elevation, Introduction

Conclusions References
Tables Figures

Back Close
Full location within the basin (UTM Easting and/or UTM Northing), and a variable related to slope, aspect, and/or curvature.The inclusion of elevation and geographic location within each regression as well as the strong bivariate correlations of these variables with SWE indicates that they may be consistent drivers of the spatial variability of SWE at the basin scale, representing how orographic precipitation and storm track patterns play a strong role in basin scale SWE distribution.Additionally, the inclusion of a terrain variable (slope, aspect, and/or curvature) within each model suggests that terrain characteristics are also important, but to a lesser degree, at this scale.However, given that various studies (e.g.Erickson et al., 2005;Fassnacht et al., 2012) have shown the spatial variability of snow accumulation to be described by different physiographic variables from year to year, additional years of data collection at the basin scale are needed for more complete evaluation.
Comparison of the error between model O+F and model O for WY 2011 and WY 2012 shows that model O has superior performance statistics for both years (Table 4).Model O11 and model O12 showed a similar strong performance to previous research using operational data at a comparable scale (e.g.Harshburger et al., 2010).This strong performance of the operational regression model, however, may not be representing the study area, as SNOTEL measurements have been shown to represent point locations rather than surrounding areas (Molotch and Bales, 2005) often having more snow (Daly et al., 2000), and tend to be located in areas with similar physiographic features (flat and open canopy areas located near tree line).
The spatial dataset of field-based snowpack measurements in this study is at a scale similar to remote sensing observations and modeling applications; these data and the approaches of empirical modeling (e.g.multiple linear regression) for characterizing the distribution of SWE at the basin scale can be used in those contexts for valida- physically based snow evolution modeling can provide insight for basin scale SWE distribution estimations.

Conclusions
We have used a combination of field-based and operational snow measurements to evaluate snowpack properties across the basin scale.This research was motivated by the need for additional ground truth snowpack observations at a scale that coincides with that of remote sensing observations and is especially pertinent to water resources forecasting.
A method for modeling snow density across the Cache la Poudre basin from historical snow course measurements was employed for estimating SWE from snow depth.
The independent variables of snow depth, day of year, elevation, and UTM Easting were used in a multiple linear regression model to estimate snow density.Statistics showed strong performance of SWE calculated from snow depth observations using the snow density model, and model validation suggests the model is transferable to independent data within the bounds of the original dataset.The methods here provide a pathway for estimating SWE from snow depth measurements, which is especially useful when evaluating snowpack properties at the basin scale, where time consuming field-based measurements of SWE are often not feasible.
The spatial variability of SWE within the Cache la Poudre basin was analysed using field-based and operational snowpack measurements.Bivariate relations of SWE and snow depth with terrain and canopy variables show that elevation, UTM Easting, and UTM Northing are most strongly correlated with SWE, and thus drive spatial variability at this scale.Multiple linear regression models were developed for WY 2011 and WY 2012 using both a combined dataset of field-based and operational measurements (model O+F ) and a dataset of operational measurements only (model O ).The continuity of field-based snowpack measurements, as provided within this study, is essential given the assumption of non-stationarity from hydroclimatic change (Milly et al., Introduction Conclusions References Tables Figures

Back Close
Full  , 2007).This examination of two very different snow years may represent the bounds of extremes and possibly the limitations due to non-stationarity.Continued field measurements of the snowpack will aid advancement of remote sensing and modeling applications, but more importantly continue to provide "ground-truth" observations for evaluating the complexities and uncertainties of the changing earth system.Introduction

Conclusions References
Tables Figures

Back Close
Full Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | areas, limiting the ability of the dataset to represent the variability explained by those variables.
Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | 841 + 1.05DOY + 17.2d s + 3.53 × 10 −2 z − 1.72 × 10 −3 UTM e (Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | WY 2011 showed much more variation in snow amounts than WY 2012, which could explain the difference in the RMSE.Additionally, the greater number of measurement locations (n = 127) in WY 2012 compared to WY 2011 (n = 51) could further explain the difference in R 2 between model O+F11 and model O+F12 .Given this difference in field-based sampling locations, Discussion Paper | Discussion Paper | Discussion Paper | a reduced model O+F for WY 2012 was developed including only WY 2012 field-based measurement locations that were co-located with WY 2011 measurement locations (n = 42 snow depth measurements though a snow density model works due to the conservative nature of snow density; 52 % of snow density data values ranged from 250 to 350 kg m −3 .Using historical operational measurements for development of a regional based snow density model has implications for future field-based basin scale sampling campaigns, suggesting a sampling scheme dominated by snow depth measurements may be successful for evaluating basin scale SWE variability.The strength and utility of the Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | tion.For instance, the observed patterns of SWE variability within this study, showing to be largely driven by elevation and geographic location, could be compared to the patterns of variability observed within a physically based snow evolution model.The comparisons of the statistical relation of the snowpack with terrain based variables and Figures Back Close Full Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | 2008) and indications of more extreme conditions (IPCC

Fig. 1 .Fig. 3 .
Fig. 1.The Cache la Poudre basin above the Cache la Poudre River at Canyon Mouth Colorado Division of Water Resources gaging station.The locations of field-based and operational snow measurements are shown.The 50 % SCI (Richer et al., 2013) is indicated by the transparent light gray color within the basin.