Independent Evaluation of the SNODAS Snow Depth Product Independent Evaluation of the SNODAS Snow Depth Product Using Regional Scale Lidar-Derived Measurements Using Regional Scale Lidar-Derived Measurements

Repeated Light Detection and Ranging (LiDAR) surveys are quickly becoming the de facto method for measuring spatial variability of montane snowpacks at high resolution. This study examines the potential of a 750 km 2 LiDAR-derived dataset of snow depths, collected during the 2007 northern Colorado Cold Lands Processes Experi-5 ment (CLPX-2), as a validation source for an operational hydrologic snow model. The SNOw Data Assimilation System (SNODAS) model framework, operated by the US National Weather Service, combines a physically-based energy-and-mass-balance snow model with satellite, airborne and automated ground-based observations to provide daily estimates of snowpack properties at nominally 1 km resolution over the coter-10 minous United States. Independent validation data is scarce due to the assimilating nature of SNODAS, compelling the need for an independent validation dataset with substantial geographic coverage. Within twelve distinctive 500m × 500 m study areas located throughout the survey swath, ground crews performed approximately 600 manual snow depth measurements 15 during each of the CLPX-2 LiDAR acquisitions. This supplied a dataset for constraining the uncertainty of upscaled LiDAR estimates of snow depth at the 1 km SNODAS resolution, resulting in a root-mean-square di ﬀ erence of 13 cm. Upscaled LiDAR snow depths were then compared to the SNODAS-estimates over the entire study area for the dates of the LiDAR ﬂights. The remotely-sensed snow depths provided a more 20 spatially continuous comparison dataset and agreed more closely to the model estimates than that of the in situ measurements alone. Finally, the results revealed three distinct areas where the di ﬀ erences between LiDAR observations and SNODAS estimates were most drastic, suggesting natural processes speciﬁc to these regions as causal inﬂuences on model uncertainty. 25


Introduction
Meltwater from mountain snowpacks is an important ecological component of Earth's water cycle. However, quantifying the amount of water stored in a snowpack from year to year remains difficult. Millions of people in the western United States rely on water that descends from the Rocky Mountains, where over 70 % of the annual water supply 5 is delivered from melting snow (Carroll et al., 2006). Hydrologic models provide important information that assist water managers in mitigating flood disasters and water shortages each year.
The primary goal of most hydrologic snow models is to provide estimates of snow water equivalent, or SWE, over large mountain regions, but in addition most models 10 include routines to estimate secondary snow properties. The methods used to estimate snowpack characteristics such as depth and density vary between models; some use empirical methods from available historical data, while others are more physics-based. Even so, SWE is but a function of depth and density, and if validation is achieved for either of these so-called secondary model components, then higher confidence can 15 be placed into corresponding SWE estimates. Since snow depth varies considerably more than bulk density over space  and is also inherently easier to measure, this study purports to examine the snow depth prediction component of a gridded, spatially-distributed snow model.
Due to the resolution capability and gridded nature of most distributed snow models, 20 many small-scale features that affect spatial variability are either averaged or not considered, influencing the bulk estimates of total SWE and overall depth in each grid cell (Marchand and Killingtveit, 2005). Nevertheless, sub-grid spatial properties have been shown to have a significant effect on the accuracy of spatially distributed snow models (Luce et al., 1999;Liston, 2004;Skaugen and Randen, 2013), but the datasets required 25 for parameter estimation and optimization are small and spatially sparse in mountain, tundra and shrubland environments (Elder et al., 1991;Sturm et al., 2001a, b;Hiemstra et al., 2002;Liston and Sturm, 2002;. Considerable

3143
This document was originally published by Copernicus Publications in Cryosphere Discussions. This work is provided under a Creative Commons Attribution License. Details regarding the use of this work can be found at: http://creativecommons.org/ licenses/by/3.0/. doi: 10.5194/tcd-8-3141-2014 variability in the spatial snow distribution can be introduced through the interaction between wind and snow with terrain and vegetation (Elder et al., 1991;Blöschl, 1999;Liston et al., 2007). In fact, wind has been shown to be the dominant influence on spatial variability of snow in complex terrain (Pomeroy et al., 1993;Winstral et al., 2002;Sturm and Wagner, 2010). Without prior knowledge of the spatial snow distribution in 5 a given area, arbitrary manual snow measurements will not provide accurate estimates of snow depth in alpine regions (Elder et al., 1991;Anderton et al., 2004;Erickson et al., 2005).
Various studies have shown that LiDAR (Light Detection and Ranging) surveys can provide spatial information on mountain snow depths at high-resolution over large areal 10 extents that comprise various physiographic regimes (Hopkinson et al., 2004;Deems et al., 2006;McCreight et al., 2012). The first Cold Lands Processes Experiment (CLPX-1) of 2002-2003 in the Colorado Rocky Mountains, was the first large-scale coordinated study to use LiDAR acquisitions for the assessment of snow properties over a range of areas . Since then, numerous campaigns have used 15 LiDAR to quantify spatial variability of snow depths in mountain terrain. Deems et al. (2006) used fractal analysis of the CLPX-1 LiDAR snow depths to determine scalebreaks, while Trujillo et al. (2007) found that spatial distributions of snow depth are strongly controlled by both wind redistribution and vegetation interception of snow over uneven surface topography in five of the CLPX-1 intensive study areas. More recently, 20 LiDAR has been used with simple statistical models to determine scale invariance due to vegetation and wind direction (Trujillo et al., 2009) as well as to verify high-resolution dynamical snow models (Mott et al., 2011).
Even though depths can vary greatly over space in a snow pack, the overall distribution of snow has been found to exhibit spatial similarities from year to year (Hiem- 25 stra et al., 2006;Sturm and Wagner, 2010;Winstral and Marks, 2014). Repeated LiDAR surveys throughout single seasons  and over multiple seasons (Deems et al., 2008) have found similar results through fractal analyses of the snow depth distributions. By comparing findings 3144 This document was originally published by Copernicus Publications in Cryosphere Discussions. This work is provided under a Creative Commons Attribution License. Details regarding the use of this work can be found at: http://creativecommons.org/ licenses/by/3.0/. doi: 10.5194/tcd-8-3141-2014 from large-scale LiDAR snow depth surveys to operational hydrologic models, we can pinpoint causes of any shortcomings and subsequently refine model results.
Developed by the National Operational Hydrologic Remote Sensing Center (NOHRSC) and first operationally implemented in 2004, the Snow Data Assimilation System (SNODAS) estimates various snow properties by merging satellite, airborne, 5 and ground-based snow data with modeled approximations of snow cover (Barrett, 2003). Historical model output from SNODAS is stored and archived at the National Snow and Ice Data Center (NSIDC) in Boulder, Colorado for every day that the model has been executed since its inception. These eight snow properties are the primary estimates that are made available to the public: A large portion of model fidelity is directed towards SWE prediction rather than any of the other model outputs because the amount of total water storage within snow-20 packs is far more important for water managers. The physically-based energy-and mass-balance NOHRSC Snow Model (NSM), described by Carroll et al. (2006), is the primary component of SNODAS, while an assimilation step gives analysts the ability 3145 This document was originally published by Copernicus Publications in Cryosphere Discussions. This work is provided under a Creative Commons Attribution License. Details regarding the use of this work can be found at: http://creativecommons.org/ licenses/by/3.0/. doi: 10.5194/tcd-8-3141-2014 to decide every day whether to augment the model estimates with any available remote sensing or ground-based measurements. Ultimately, the final model product has a spatial resolution of approximately 1 km 2 over the coterminous United States.
Independent validation data for SNODAS is scarce as a consequence of the framework's data assimilating nature which ensures that all available data at the model scale 5 (i.e. 1 km 2 ) is used to adjust estimates of the NSM (Barrett, 2003). An alternative validation method has been to perform comparisons of SNODAS with other hydrologic models and satellite remote sensing products. Rutter et al. (2008) compared various NSM properties with two energy-balance snow models, but found difficulty in constraining model uncertainties due primarily to the high sub-grid spatial variability exhibited in 10 mountain snowpacks. Other studies have used SNODAS as the validation source for large-scale hydrologic models such as the Noah land surface model (Barlage et al., 2010), and SWE retrieval using satellite-based microwave radar remote-sensing platforms (Azar et al., 2008).
To our knowledge, only two validation studies of SNODAS' performance have been 15 conducted using independent datasets and each of those studies relied on extensive in situ measurement campaigns. Clow et al. (2012) performed snow surveys of snow depth within 45 SNODAS pixels over a three month period in 2007. The results revealed that SNODAS performed satisfactorily for predicting snow depth in forested areas, but depth estimates in alpine areas were poor in comparison to manual measurements 20 chiefly due to sub-grid scale variability from wind redistribution of snow. This discrepancy was addressed by applying a correction factor to account for wind redistribution of snow in the wind-affected alpine areas. In another study, Anderson (2011) intensively sampled three SNODAS pixels in the mountains just north of Boise, Idaho over the course of two winter seasons and found that SNODAS slightly under predicted 25 snow depths in heavily-forested areas but maintained reasonable estimates of SWE overall. Each of the studies required an enormous amount of manpower and time to obtain the independent datasets for proper comparison, but came to somewhat differ-

3146
This document was originally published by Copernicus Publications in Cryosphere Discussions. This work is provided under a Creative Commons Attribution License. Details regarding the use of this work can be found at: http://creativecommons.org/ licenses/by/3.0/. doi: 10.5194/tcd-8-3141-2014 ent conclusions about the model performance most likely due to the individual locations of the collected data (Idaho and Colorado, USA).
As the technology has become more widespread over the last decade, LiDAR for snow research has become increasingly popular. The advantages of LiDAR for spatially characterizing snow depths over large remote areas are finally being used to assess 5 lower resolution hydro-meteorologic snow models. Melvold and Skaugen (2013) used six parallel 500 m × 80 km LiDAR surveys, each separated by 10 km, to investigate the Norwegian operational temperature index snow model, seNorge. After upscaling the LiDAR-derived 2 m resolution snow depths to the spatial resolution of the 1 km 2 gridded model output, the modeled results were found to accurately represent the remote 10 sensing estimates despite the lack of sub-grid spatial information within the model structure.
To our knowledge, this study for the first time examines the performance of a physically-based hydro-meteorologic snow model in complex terrain using airborne LiDAR-derived snow depths that have been subjected to ground validation error analy- 15 sis.

Study area
The second Cold Lands Processes Experiment campaign (CLPX-2, 2006-2008) was a multi-faceted mission designed to cover a much larger coincident extent than the previous campaign (CLPX-1, 2002(CLPX-1, -2003 three years prior. The primary objective 20 of CLPX-2 was to acquire snow volume backscatter measurements from NASA's POLSCAT (POLarimetric SCATterometer) airborne Ku-band radar system and the necessary ground truth measurements ) for validation of the proposed NASA Snow and Cold Land Proceses (SCLP) and ESA Cold Regions Hydrology Highresolution Observatory (CoreH 2 O) satellite missions (Rott et al., 2010). The airborne 25 LiDAR portion of the campaign was intended to be an ancillary validation dataset for the radar measurements.

3147
This document was originally published by Copernicus Publications in Cryosphere Discussions. This work is provided under a Creative Commons Attribution License. Details regarding the use of this work can be found at: http://creativecommons.org/ licenses/by/3.0/. doi: 10.5194/tcd-8-3141-2014 Three Intensive Observation Periods (IOPs) were organizeded over a 9 km × 84 km rectangular swath to the south and east of the town of Steamboat Springs in northern Colorado, USA (Fig. 1). During both IOP-1 (early December 2006) and IOP-3 (late February 2007), airborne LiDAR surveys were performed to provide high-resolution surface elevation change datasets to aid in the POLSCAT validation process. Covering  surfaces owing to the minimal penetration depth on the order of only 1 cm (Deems et al., 2006). The pulse rate of 32 500 Hz, combined with the aircraft's speed, altitude, and scan rate, resulted in raw point clouds with nominal point spacings of 2.0-2.5 m, depending on surface roughness, canopy coverage and scan angle relative to the aircraft. 5 The LiDAR vendor filtered vegetation returns from ground returns using a minimum block mean algorithm and proprietary software to create vegetation-filtered point clouds for each flight with updated nominal point spacings of 2.5-3.0 m, again depending on the terrain, canopy cover and scan angle. Various alternative filtering algorithms were explored during the course of this study, but the decision was ultimately made to utilize 10 the vendor-filtered data in order to maintain consistency over the large variety of landscapes. Next, we applied the open-source Points2Grid interpolation tool, employing an inverse distance weighting scheme, to produce a 5 m Digital Surface Model (DSM) for both of the vegetation-filtered point clouds. Because the CLPX-2 LiDAR scans were never acquired over an absolutely snow-free surface, as many of the higher elevations 15 had already received snow by 3 December, the interpolated surfaces were differenced to provide a raster of the estimated change in total snow height between 3 December and 22 February at 5 m resolution (Fig. 1). This 5 m gridded product of LiDAR-estimated changes in snow depth will hereafter be referred to as L Est . Though less dense than the original CLPX-1 acquisitions, the CLPX-2 LiDAR footprint covered many more ter-20 rain, vegetation, and snowpack classes, thereby providing a useful comparison tool for hydrologic snow models over large spatial extents.

In situ measurements
All remote sensing methods are subject to an appreciable amount of measurement uncertainty which should be quantified by ground validation if possible. Originally in-25 tended as the primary ground truth dataset for each POLSCAT radar acquisition over the 2006-2007 winter season , the CLPX-2 intensive manual measurement campaign was arranged and completed by a team of 12-15 researchers 3149 during each IOP. Twelve 500 m × 500 m "hourglass" transects, henceforth referred to as HG sites with each comprised of 47-50 evenly spaced snow depth measurements, were manually sampled during IOP-1 and IOP-3 within a day of each of the CLPX-2 LiDAR acquisitions (Fig. 1). The HG sites were chosen to represent physiographically distinctive regions of the CLPX-2 survey swath. Ground crews made measurements at 5 preprogrammed waypoints loaded onto mapping grade handheld GPS units in order to maintain location consistency for each survey. We estimate the resulting relative pointto-point horizontal uncertainty between the HG surveys to be less than 2 m while the HG transect locations themselves can be approximated to 7 m accuracy in absolute space. The repeated HG depth measurements were differenced to provide a similar 10 comparison metric of snow depth change to the L Est dataset.

SNODAS snow depths
SNODAS estimates of snow depth were downloaded from the NSIDC for the two dates of the CLPX-2 LiDAR acquisitions (3 December 2006 and 22 February 2007), then spatially referenced to the UTM coordinate projection. The two rasters of snow depth 15 were differenced to provide 1 km gridded model estimates of snow depth change, hereafter referred to as M Est . Figure 2 depicts M Est over the area surrounding the LiDAR swath, along with the locations of all nearby SNOTEL stations that can be used for model assimilation. In view of the fact that snow was already present over much of the study area on 3 December, it is paramount that no snow melt occurred between the 20 survey dates in order to properly assess the snow depth component of SNODAS using LiDAR estimates alone. Therefore, we examined the SNODAS estimates of snow melt due to incoming solar radiation and sublimation (Fig. 3) and the six nearby SNOTEL stations to determine whether the LiDAR estimates of snow depth change are a viable comparison tool for SNODAS. As can be seen in Fig. 3

SNODAS/in situ measurement comparison
To provide a link to the previous ground-based SNODAS validation studies, we examined the ability of manual measurements from the twelve HG sites to represent model estimates of snow height change. Mentioned previously, Clow et al. (2012) averaged depth measurements from snow surveys carried out within 45 individual SNODAS pix-5 els to perform model validation. We employed the same method to assess SNODASpredicted snow depth changes using the CLPX-2 in situ HG transects. The mean change in snow depth at each HG site was found with an associated interquartile range, similar to the method used by Clow et al. (2012). Then, a coincident 1 km 2 SNODAS estimate of depth change was constructed over each HG transect site from 10 the areal coverage fraction of the overlapping model pixels, creating an area-weighted average of M Est centered over each HG measurement site. This spatial averaging was performed because the CLPX-2 campaign was not designed during the planning phase to be a validation source for SNODAS and the HG transects were therefore not aligned within individual model pixels.

Characterizing LiDAR uncertainty
The 750 km 2 CLPX-2 LiDAR dataset (L Est ) overlaps 980 individual SNODAS pixels completely 1 , supplying a statistically robust dataset for determining contributing factors to model uncertainty. However, LiDAR measurements are fundamentally estimates themselves and require uncertainty assessments, which was available from the HG 20 in situ measurements. To account for the horizontal position uncertainty in both the HG and L Est datasets, the 5 m gridded L Est snow depth changes were averaged in a 10 m 1 Though previously stated as nominally 1 km 2 , the actual resolution of SNODAS is 30 arcseconds because the model is implemented in the geographic coordinate system (Barrett, 2003). At the latitude of CLPX-2, that is equivalent to approximately 0.83 km 2 .

3151
This document was originally published by Copernicus Publications in Cryosphere Discussions. This work is provided under a Creative Commons Attribution License. Details regarding the use of this work can be found at: http://creativecommons.org/ licenses/by/3.0/. doi: 10.5194/tcd-8-3141-2014 radius around each reported in situ measurement location and treated as a separate point measurement for comparison purposes. To perform the model comparison, the L Est 5 m snow depth changes were binned into the spatial extents of the 980 overlaid SNODAS M Est grid cells. Statistics were calculated within each 1 km pixel, resulting in a mean, standard deviation, and interquar-5 tile range of LiDAR-estimated snow depth change over the CLPX-2 study area at the SNODAS model resolution. The mean L Est snow depth changes are portrayed in Fig. 4.

Results and discussion
We examined how well the twelve averaged HG measurement transects were able to substantiate the SNODAS estimates of snow depth change. Figure 5a shows the 10 relationship between M Est and the mean manually-measured changes in snow depth at each of the HG sites, while Fig. 5b shows the difference between the HG in situ measurements and M Est as a function of the mean HG measurements themselves. The trend of this limited dataset suggests that as the mean snow depth within a model pixel increases above approximately 40 cm, the ability of SNODAS to estimate the 15 amount of total snow depth change decreases substantially. Also, it is not clear from the small sample size of the in situ data what physical factors could be influencing such discrepancies and a much more spatially extensive dataset, such as the CLPX-2 LiDAR, is required for determining the underlying causes of model error.
The exhaustive CLPX-2 in situ HG measurement campaign provided an ideal dataset 20 for limiting uncertainty in the large-scale LiDAR surveys of 3 December 2006 and 22 February 2007. The changes in snow depth as measured by the standard probing method and interpolated from the LiDAR surveys were compared throughout all twelve HG sites individually, and then averaged on a site by site basis (Fig. 6), again similar to the averaging scheme used by Clow et al. (2012). In this instance, the RMS 25 difference of 12.9 cm between mean HG and mean L Est snow depth changes for all twelve CLPX-2 HG site locations is well within the bounds of conventional airborne 3152 This document was originally published by Copernicus Publications in Cryosphere Discussions. This work is provided under a Creative Commons Attribution License. Details regarding the use of this work can be found at: http://creativecommons.org/ licenses/by/3.0/. doi: 10.5194/tcd-8-3141-2014 LiDAR uncertainty estimates (Baltsavias, 1999;Hodgson and Bresnahan, 2004). Also, the averaged HG measurements and L Est observations are highly correlated with a linear r-squared coefficient of 0.942. It can be seen in Fig. 6 that the mean L Est snow depth change consistently underestimates the mean HG measurements by 5-30 cm, with deeper snowpacks tending to be more sensitive to higher spatial variability and 5 in situ under sampling. A variety of factors likely attribute to the underestimates by the LiDAR-derived snow depth change when compared to the HG point measurements. One viable explanation would be the difference in measurement support between LiDAR and in situ measurements. While the manual depth measurement support was less than one centimeter 10 (the size of the depth probe tip), the LiDAR-derived snow depth change was averaged over all the 5 m pixels within a 10 m radius of the reported in situ measurement location. This smoothing of the snow height change could result in lower reported values than that of the point measurement when scale breaks in snow depth are present that are lower than the LiDAR-averaged areal extent. Even so, the correlation between 15 the ground-based and airborne data shows that we should be able to trust the 1 km spatially-averaged, LiDAR-derived estimates of snow depth change within a range of approximately ±13 cm.
The comparison between M Est and mean L Est snow depth change within the model pixels (Fig. 7) determined an r 2 = 0.72, signifying a reasonably strong correlation be-20 tween the two estimate datasets. Since snow melt between the LiDAR flights was found to be an insignificant portion of the snowpack evolution (Fig. 3), the actual changes in snow depth over the study area were primarily influenced by accumulation, densification, and redistribution factors (sublimation effects are beyond the scope of this paper).
To investigate the primary cause of disagreement between M Est and L Est , seven po-25 tential physiographic parameters were culled from the LiDAR data to perform a regression analysis. In addition to the LiDAR-derived snow depth changes and the vegetationfiltered elevations, vegetation height and canopy coverage across the survey swath was calculated at 5 m resolution using both the raw and vegetation-filtered December

3153
This document was originally published by Copernicus Publications in Cryosphere Discussions. This work is provided under a Creative Commons Attribution License. Details regarding the use of this work can be found at: http://creativecommons.org/ licenses/by/3.0/. doi: 10.5194/tcd-8-3141-2014 LiDAR point returns. Vegetation heights and elevations were each upscaled to the 1 km SNODAS resolution in a similar fashion to the LiDAR snow depth change, while the vegetation density was calculated by finding the number of 5 m pixels within each 1 km SNODAS grid cell that contained LiDAR first returns greater than 50 cm above the filtered surface. Lastly, the interquartile range of the 1 km averaged (upscaled) variables 5 was determined to result in a group of 7 individual predictor variables for regression analysis. The upscaled snow depth changes were overwhelmingly found to best predict the 15 discrepancy between M Est and L Est , indicating that none of the other six factors singularly influenced SNODAS performance over the study area. Figure 8 shows the plot of mean L Est snow depth changes against the M Est − L Est differences within each model pixel. Within that plot, the pink vertical and blue horizontal stripes between −13 and +13 cm on each axis represent the minimum attainable resolution of the LiDAR-derived 20 snow depth values determined from the HG ground-based measurement data. Three regions have been circled in the figure, each corresponding to portions of the difference dataset that were found to be outside the uncertainty levels of the LiDAR-derived changes in snow depth.

3154
This document was originally published by Copernicus Publications in Cryosphere Discussions. This work is provided under a Creative Commons Attribution License. Details regarding the use of this work can be found at: http://creativecommons.org/ licenses/by/3.0/. doi: 10.5194/tcd-8-3141-2014 Contrasting the images of M Est − L Est (Fig. 9) and mean L Est (Fig. 4) reveals the geographic locations of the three regions within the survey swath containing the greatest SNODAS and LiDAR disagreements. Within these three regions specific physiographic factors are likely the cause of greater relative discrepancy.

5
The region with the lowest annual snow totals is comprised of pixels that SNODAS estimated to have a larger positive change in snow depth between LiDAR acquisitions. However, the LiDAR snow depth changes within these pixels are well below the trusted LiDAR uncertainty level (the pink vertical stripe). These pixels are located in the North Park region of the survey area, where the flat landscape is densely populated by low 10 sagebrush (≈ less than 30 cm) and high winds frequently scour the snow above and near the height of the sage throughout the winter. The snow that remains is subsequently packed between the low vegetation and the snow height changes very little throughout the year once it has reached a height similar to the sagebrush. SNODAS does incorporate a sublimation factor due to wind into the accumulation model, but the 15 nearest SNOTEL station used in the assimilation step is located nearly 15 kilometers to the southwest in a much different landscape so wind speeds and directions are not well-represented in the area. The locations of the region #1 pixels are roughly delineated in Fig. 9. 20 Pixels that comprise region #2 in Fig. 8 are where snow depths are similarly estimated by SNODAS to have changed more than observed by the LiDAR. However, the geographic location of the pixels are in a region with higher snow accumulation totals, which are above the lower LiDAR uncertainty level of 13 cm. Again delineated in Fig. 9, these pixels are nestled directly to the east of Rabbit Ears Pass where the Columbine 25 SNOTEL station provides assimilation data for SNODAS. Since the relative error of 3155 This document was originally published by Copernicus Publications in Cryosphere Discussions. This work is provided under a Creative Commons Attribution License. Details regarding the use of this work can be found at: http://creativecommons.org/ licenses/by/3.0/. doi: 10.5194/tcd-8-3141-2014 the LiDAR observations is small and a large altitudinal gradient can be seen in L Est changes in snow depth (Fig. 4), this discrepancy can possibly be attributed to SNODAS over-distributing the SNOTEL information to areas of lower elevations and vegetation types. Future study on regions such as this would be important for determining optimal precipitation forcings by the SNODAS data assimilation process.

Region #3: Rabbit Ears Pass
Finally, the region #3 pixels represent an area where the upscaled LiDAR changes in snow depth are significantly larger than the SNODAS estimates. These pixels occur primarily in topographically complex areas with exceptionally high snow totals and dense coniferous forests, once again outlined in Fig. 9. The probable controlling factor of un-10 derestimation by SNODAS in this region is the sub-kilometer scale heterogeneity of snow distribution caused by both vegetation and topography. SNODAS has been found to underestimate snow depths in similar forested alpine terrain (Anderson, 2011), so this result is not unexpected.

15
Over the past decade, high resolution snow depth information has become a highly sought-after data product by snow researchers and many scientific questions have been addressed using the spatial continuity and extent provided by LiDAR surveys. This study first examined the ability of ground-based measurements to constrain remote sensing uncertainty, and in turn compared the remote sensing estimates to an 20 operational hydrologic model for validation purposes. In this case, the CLPX-2 ground truth campaign was vitally important for quantifying uncertainty in the LiDAR snow depth estimates, revealing the necessity of similar campaigns to complement future LiDAR remote sensing missions.

3156
From the comparison study, three distinct regions were extracted that exhibited greater disagreement than could be explained by LiDAR estimate uncertainty alone. Distinct physiographic characteristics within these three regions ultimately affected the accuracy of the SNODAS predictions of snow height change between the two LiDAR acquisitions. 5 To further investigate model performance, more studies are needed from subsequent LiDAR-derived snow depth datasets that focus on the performance of SNODAS as a function of the distance from SNOTEL stations. Additionally, micro-scale wind redistribution effects could be applied within the model structure to assist in areas where blowing snow transport is a major cause of spatial variability. Finally, large-scale coin-10 cident density surveys would allow model validation with LiDAR-derived snow depths as well as in situ estimates of SWE, for which SNODAS is likely to be more accurate than compared with depth alone.

Author contributions
Adam Winstral supplied modeling expertise, Kelly Elder produced the archived LiDAR