LiDAR measurement of seasonal snow accumulation along an elevation gradient in the southern Sierra Nevada , California

We present results from snow-on and snow-off airborne-scanning LiDAR measurements over a 53 km 2 area in the southern Sierra Nevada. We found that snow depth as a function of elevation increased approximately 15 cm per 100 m, until reaching an elevation of 3300 m, where depth sharply decreased at a rate of 48 cm per 100 m. Departures from the 15 cm per 100 m trend, based on 1 m elevation-band means of regression residuals, showed slightly less steep increases below 2050 m; steeper increases between 2050 and 3300 m; and less steep increases above 3300 m. Although the study area is partly forested, only measurements in open areas were used. Below approximately 2050 m elevation, ablation and rainfall are the primary causes of departure from the orographic trend. From 2050 to 3300 m, greater snow depths than predicted were found on the steeper terrain of the northwest and the less steep northeast-facing slopes, suggesting that ablation, aspect, slope and wind redistribution all play a role in local snow-depth variability. At elevations above 3300 m, orographic processes mask the effect of wind deposition when averaging over large areas. Also, terrain in this basin becomes less steep above 3300 m. This suggests a reduction in precipitation from upslope lifting and/or the exhaustion of precipitable water from ascending air masses. Our results suggest a cumulative precipitation lapse rate for the 2100–3300 m range of about 6 cm per 100 m elevation for the accumulation period of 3 December 2009 to 23 March 2010. This is a higher gradient than the widely used PRISM (Parameter-elevation Relationships on Independent Slopes Model) precipitation products, but similar to that from reconstruction of snowmelt amounts from satellite snow-cover data. Our findings provide a unique characterization of the consistent, steep average increase in precipitation with elevation in snow-dominated terrain, using high-resolution, highly accurate data and highlighs the importance of solar radiation, wind redistribution and mid-winter melt with regard to snow istribution.


Introduction
In mountainous regions of the western United States, snowmelt is the dominant contributor to surface runoff, water use by vegetation, and groundwater recharge (Bales et al., 2006;Earman and Dettinger, 2011).Because of their importance and vulnerability of mountain snowpacks in a warmer climate, several researchers have recently developed scenarios for changes in annual and multiyear mountain water cycles, including trends in water storage and runoff, groundwater recharge, and feedbacks with vegetation (Peterson et al., 2000;Marks et al., 2001;Lundquist et al., 2005;Maxwell and Kollet, 2008;Barnett et al., 2008;Anderson and Goulden, 2011;Trujillo et al., 2012).
Given the challenges in measuring the spatial distribution of mountain precipitation, the processes controlling its distribution remain poorly understood.However, since the large majority of precipitation in the middle and upper elevations of the southern Sierra Nevada falls and accumulates as snow, with limited ablation through much of the winter, we can P. B. Kirchner et al.: LiDAR measurement of seasonal snow accumulation examine snow accumulation to assess processes governing the distribution of precipitation.
Snow accumulation across the mountains is primarily influenced by orographic processes, involving feedbacks between atmospheric circulation, terrain and the geomorphic processes of mountain uplift, erosion and glaciation on the Earth's surface (Roe, 2005;Roe and Baker, 2006;Pedersen et al., 2010;Kessler et al., 2006;Stolar et al., 2007;Galewsky, 2009;Mott et al., 2014).Orographic precipitation is well documented and central to determining the amount of snow water equivalent (SWE) in mountainous regions.The Sierra Nevada, a major barrier to land-falling storms from the Pacific, is ideally oriented to produce orographic precipitation, and exerts a strong influence on the upslope amplification of precipitation and the regional water budget (Pandey et al., 1999).Despite this well-developed conceptual understanding, our ability to apply this knowledge at spatial and temporal scales relevant to questions of regional climate and local water-supply forecasting are limited by a lack of accurate precipitation measurements across mountains (Viviroli et al., 2011).Additionally, long, narrow land-falling bands of extratropical Pacific water vapor, referred to as atmospheric rivers, frequently deposit large fluxes of orographic precipitation as they ascend over the Sierra Nevada (Neiman et al., 2008;Ralph and Dettinger, 2011).Atmospheric rivers deposit approximately 40 % of total winter snowfall in the Sierra Nevada, linking ocean-atmosphere interactions and the terrestrial water balance (Dettinger et al., 2011;Guan et al., 2012Guan et al., , 2013a)).
Current mountain-basin operational SWE estimates are made with a limited set of snow-course and continuous in situ point measurements from snow pillows.Measurements at these index sites are used to develop statistically based runoff estimates for the subsequent spring and summer.While this approach has provided operationally robust predictions in years near the long-term normal, snow accumulation both varies from year to year and changes in response to long-term climatic conditions, and, in recent decades, has trended outside the statistical normal (Milly et al., 2008).Hence our current methods are becoming less reliable, and accurate predictions require a more comprehensive approach to understanding the processes affecting precipitation and the probabilities of extremes (Rahmstorf and Coumou, 2011).
Accurate estimates of the amount and spatial distribution of both precipitation and SWE are essential given the shift toward spatially distributed models for forecasts of runoff, moisture stress and other water-cycle components (Rice et al., 2011;Meromy et al., 2012).Current operational measurements for precipitation and SWE are limited by scale and by the heterogeneity of snow-accumulation processes, and do not provide spatially representative values (Viviroli et al., 2011;Bales et al., 2006;Grünewald and Lehning, 2011).Uncertainty in watershed-scale SWE and precipitation estimates result in part from the lack of measurements at both the rain-snow transition and highest elevations, as well as the lack of representative measurements across different slopes, aspects and canopy conditions (Molotch and Margulis, 2008).
Remotely sensed snow properties from satellites and aircraft are used in research, and on a limited basis in forecasts.In both cases these measurements can be blended using statistical or spatially explicit models to produce discharge forecasts (Rice et al., 2011;Molotch et al., 2004;Fassnacht et al., 2003;Bales et al., 2008;Kerkez et al., 2012).A recent review highlighted the promise of aircraft LiDAR measurements for snow-depth mapping at high spatial resolution and vertical accuracy, using repeat snow-on and snow-off LiDAR flights (Deems et al., 2013).The emergence of quality research data sets for snow mapping offers opportunities to assess LiDAR accuracy and coverage in complex, forested terrain, and also its potential for providing a much-needed spatial "ground truth" for watershed-scale snow depth (Harpold et al., 2014).
Research reported in this paper was aimed at determining the influences of terrain and orographic precipitation on patterns of seasonal snow accumulation along a 1650 m elevation gradient in the southern Sierra Nevada.Three questions were posed in this research: (i) what is the magnitude of the average elevation lapse rate for snow accumulation, (ii) what is the variability in snow accumulation at each elevation along an elevation gradient, and (iii) to what extent do local terrain and wind redistribution influence this pattern?It was also our aim to evaluate consistency between LiDARestimated SWE and prior model-based estimates of accumulated SWE and total precipitation.

Methods
Our approach involved analysis of (i) LiDAR-based snowdepth estimates derived from two LiDAR acquisitions, one when the ground was snow-free and one near peak snow accumulation on 23 March; (ii) continuous ground-based measurements of snow depth, SWE, wind speed and air temperature, plus operational bright-band radar observations; and (iii) model estimates of SWE and precipitation.The LiDAR data were used to estimate snow depth across the study area at a 1 m 2 spatial resolution in open areas without canopy cover.The ground measurements were used in interpreting the spatial patterns and in estimating SWE, and the brightband radar in determining the rain-snow transition elevation for precipitation events, an important metric for interpreting snow depth and SWE along elevation gradients.

Location
Our study area is centered at approximately 36.5 • N, 118.7 • W and includes the 53.1 km 2 area covered by the two LiDAR flights in the southeastern part of the 135 km 2 Marble Fork of the Kaweah River watershed, located in Sequoia National Park in the southern Sierra Nevada, California (Fig. 1).Elevations of the LiDAR acquisition were 1850-3494 m, with aspects predominantly trending northwest, about orthogonal to the regions southwest prevailing storm tracks.The land features include glaciated lake basins, cirques and stepped plateaus at the highest elevations.Soils are characterized by moraine deposits and well-drained granitic soils at the lower elevations, and rock outcrops with pockets of course shallow soil at the higher elevations.The vegetation cover below 3000 m consists primarily of coniferous forests that transition with increasing elevation from a giant sequoia grove, to mixed-conifer forests, and to red fir forests.Above 3000 m there are increasing areas of bare rock with subalpine forests and alpine meadows in locations with soil (Fig. 2b).

LiDAR altimetry
Airborne-scanning LiDAR altimetry was collected by the National Center for Airborne Laser Mapping (NCALM) using an Optech Gemini ® ALTM 1233 airborne-scanning laser (Zhang and Cui, 2007).The two campaigns were conducted in the 2010 water year: 23 March for snow-covered, and 15 August for snow-free conditions (Harpold et al., 2014).The instrument settings used for acquisition generated an average point density greater than 10 points per square meter, and a fine-scale beam-sampling footprint of approximately 20 cm (Table 1).Ground points were classified by NCALM through iteratively building a triangulated surface model with discrete points classified as ground and non-ground (Shrestha et al., 2007;Slatton et al., 2007).The nominal horizontal and vertical accuracy for a single flight path are 0.5 and 0.11 m, respectively, but higher accuracy was likely achieved, particularly where flight paths overlapped.
A digital surface model (DSM) was created by using firstreturn points and discarding outliers > 100 m (tallest trees are approximately 85 m) and returns below −0.1 m, where values in the range of −0.1 to 0 m were classified as 0. A continuous-coverage bare-earth digital elevation model (DEM) was created through kriging of ground points using a linear variogram with a nugget of 15 cm, a sill of 10 m, a range of 100 m, and a search radius of 100 m, where the minimum number of points was five (Guo et al., 2010).We used a 1 m 2 gridded model for representing our data, as this is the smallest footprint that most closely matches the expected beam sampling footprint and uncertainty in horizontal accuracy.After interpolation, digital models of mean elevation and point-return density grids were georegistered to a common grid for snow-on and snow-off flights.The average point-return densities were 8 m −2 for the surface model and 3 m −2 for the bare-earth model.Grids with no point returns in either flight, primarily under forest canopy, were not used.
The accuracy of the LiDAR altimetry was evaluated by using 352 georegistered 2.5 m × 2.5 m grid samples of the point cloud along the paved highway in the western part of the domain; because the highway is plowed regularly, surface heights do not change with snow accumulation.These samples had a bias of +0.05 m and a standard deviation of 0.07 m, which is below the estimated combined two-flight instrumental elevation error of 0.11 m (Xiaoye, 2008;Zhang and Cui, 2007).A possible explanation of the 0.05 m bias for the snow-on flight is that some sections of the road had a small amount of snow remaining after plowing.
A 1 m 2 gridded DSM of the vegetation canopy > 2 m was created by subtracting the DSM from the DEM.In order to accurately determine snow depth, values were further classified into two groups, where snow depth was either greater or less than the coincident vegetation height.This allowed us to consider, for further analysis, only snow from open slopes or where it had accumulated in the gaps between trees.To reduce the amount of error, we eliminated locations with slopes greater than 55 • , warranted by the high number of erroneous values and known issues of vertical inaccuracies due to slope angle (Schaer et al., 2007;Deems et al., 2013).Additionally, we eliminated areas with rapid annual vegetation growth that had negative snow-depth values (e.g., areas within a wet meadow).Lastly, we filtered out areas with open water, buildings and parking lots where returns were not representative of local snow accumulation.Mean snow depth for each 1 m 2 elevation band with > 100 m 2 area was computed from the snow-depth grid.Additionally, a 5 m elevation model, aggregated from the 1 m 2 bare-earth model, was produced to remove scaling biases in the analysis of slope and aspect (Kienzle, 2004;Erskine et al., 2006).

Spatial analysis
To analyze possible correlations between terrain steepness and snow distribution, we calculated the first derivative of slope and snow depth, over distances of 5-100 m, using the 1 m 2 mean snow depths and the corresponding mean slope for each 1 m elevation band, computing the correlation at 5-100 m using 5 m steps.Using the derivatives identifies transition areas.
For quantifying the combined effect of slope and aspect on snow depth, we indexed aspect on a scale of 1 to −1 using methods adapted from Roberts (1986): where V A is the aspect value; A is the azimuth variable or direction for which the calculation is being indexed to; and FA is the focal aspect, e.g., FA = 0 • is north and FA = 45 • northeast.The aspect value was further scaled by the sine of the slope angle, yielding 0 in flat terrain and approaching 1 as the mean slope increases to 90 • : where I A is aspect intensity and S the slope angle.The method of scaling the cosine of aspect by sine of slope for A = 0 • is referred to as "northness" (Molotch et al., 2004).

Ground measurements
Meteorological data were obtained from six meteorological stations in the flight area for the period from the first significant snowfall on 3 December 2009 to the 23 March Li-DAR acquisition date, henceforth referred to as the snowaccumulation period.At these stations, temperature was measured using Vaisala HMP-35 and Campbell T-108 sensors, with wind speed and direction measured using RM Young 5103 sensors.All meteorological stations measure hourly average wind speed, and two stations - Wind sensors are between 4.2 and 6.5 m above ground level, and we scaled wind speeds to 10 m using a logarithmic profile to estimate saltation thresholds: where V 10 is wind velocity at 10 m, V z is measured velocity, z is instrument height, and k the site specific roughness length.
To identify periods with the greatest potential for wind redistribution of snow, we filtered for times when temperature was below 0 • C and wind velocity was above the minimum saltation threshold of 6.7 m s −1 established by Li and Pomeroy (1997a).
Snow depth was measured continuously by 26 ultrasonic snow-depth sensors (manufactured by Judd Communications, Salt Lake City, UT) placed on meteorological stations and over or near snow pillows.These snow-depth sensors have an effective beam width of 22 • , and were mounted up to 4.6 m above the ground on a steel arm extending from a vertical steel pipe anchored to a U-channel post.This arrangement provided a snow-depth observation area of up to 2.3 m 2 over flat, bare ground, with sampling area decreasing as snow depth increases.
The LiDAR measurements, plus ground-based snowdensity measurements, were used to develop estimates of SWE versus elevation.Paired snow-depth and snow-pillow SWE measurements were part of the California Cooperative Snow Survey network, and data were acquired from the California Department of Water Resources (data available at http://cdec.water.ca.gov/, 2014) for all 16 operable stations on the western slope of the southern Sierra Nevada within 100 km north and 50 km south of the study area (Fig. 1).One snow pillow (GNF) is located 2.5 km west-southwest of the LiDAR acquisition area.Daily snow densities were estimated by dividing the daily mean SWE from the snow pillows by snow depth from the collocated ultrasonic depth sensors.To minimize the error from intermittent noise associated with snow pillows, we used the daily average SWE and did not consider measurements under a 20 cm SWE threshold.This procedure was necessary because complete snow coverage of the snow pillow is unlikely for shallow snow and the combined uncertainties of depth sensors and snow pillows can yield significant error in density measurements (Johnson and Schaefer, 2002).In addition, accumulated precipitation measurements from Alter-shielded Belfort gauges at Giant Forest (GNF), Quaking Aspen (QUA) and Charlotte Lake (CRL) and manually measured daily precipitation from Lodgepole ranger station were compared with SWE measurements to estimate total precipitation (data available at http://cdec.water.ca.gov/, 2014).All instrumental data were formatted, calibrated and gap filled by interpolation or correlation with other sensors and aggregated to daily means prior to analysis (Moffat et al., 2007).Under 1 % of the meteorological data required filtering or gap filling; snow-pillow data required slightly more (< 5 %) and snow depth required up to 20 %.Stations with data gaps >2 days with no nearby station for interpolation were not used in our analysis.

Bright-band radar
The transition elevation where hydrometeors turn from frozen to liquid, or freezing level, was determined from analysis of hourly Doppler radar data from wind profilers located upstream of the LiDAR-acquisition area.Radar reflectance is greatest, or brightest, in the altitude range where precipitation changes from snow to rain, owing to a difference in the dielectric factor for water and ice and the aggregation of hydrometeors (White et al., 2009;Ryzhkov and Zrnic, 1998).We analyzed bright-band altitudes and thus identified freezing levels from observations collected over the 2010 water year snow-accumulation period (3 December-23 March) from the three nearest upwind locations, i.e., Punta Piedras Blancas, Lost Hills, and Chowchilla, California (data available at http://www.esrl.noaa.gov/psd/data/,2014) (Fig. 1).

Model reanalysis
We calculated spatial SWE from LiDAR snow-depth measurements using mean snow-density measurements from the 16 snow-pillow sites.These values were compared with two scales of the widely used PRISM (Parameter-elevation Relationships on Independent Slopes Model) precipitation estimates, plus SWE estimates from two different MODIS-based SWE reconstruction models (Daly et al., 2008(Daly et al., , 1994;;Rice and Bales, 2011;Rice et al., 2011).Using the available 4 km and 800 m PRISM model output, we summed precipitation for December-March at the spatial extent of the LiDAR acquisition.The 4 km data were monthly for the 2010 water year and the 800 m data were monthly 30-year mean climatology.We then calculated the cumulative precipitation for each 1 m elevation band across the elevation gradient of both data sets, and aggregated values to the resolution of the comparative data.
One reconstruction data set gives 2000-2009 accumulation-period means of the entire Kaweah River watershed, calculated at a 500 m resolution, based on 300 m elevation-bin averages of MODSCAG snow-cover data, local ground-based meteorological measurements, and a temperature-index snowmelt equation that was calibrated with snow-pillow data (Rice and Bales, 2013).Fractional snow-cover was adjusted for canopy using 2 standard deviations of the elevation-band mean.
The second reconstructed SWE data were developed using the algorithm developed by Molotch (2009) and applied to the Sierra Nevada as described in Guan et al. (2013b).Fractional snow cover was adjusted for canopy using vegetation data from the Global Forest Resource Assessment 2000.
The Guan et al. (2013b) values were a subset taken from a December-March Sierra Nevada-wide calculation.The primary difference between this method and the one developed by Rice and Bales (2013) is that the Guan et al. (2013b) method includes an explicit treatment of all radiative and turbulent fluxes, whereas the Rice and Bales (2013) method uses a degree-day melt-flux calculation.

LiDAR-measured snow depth
Of the 53.1 km 2 planer footprint of the LiDAR survey, 0.8 km 2 was over water or in areas that exceeded filter thresholds of the DSM.An additional 0.01 km 2 of area with slope > 55 • , roads and buildings, and rapidly growing meadow vegetation were also removed from the analysis.The total snow-covered area where both LiDAR and ground returns were available at a density > 1 m −2 was 40.2 km 2 , and of this area 5.0 km 2 was under canopy and also eliminated from this analysis.This left an area of 35.2 km 2 remaining for analysis, and of this < 0.2 km 2 , mostly below 2300 m, was snow-free.Mean snow depth, measured by LiDAR, increased with altitude from 1850 to 3300 m elevation, with depths decreasing above 3300 m (Fig. 2a).Up to 3300 m, snow depth shows a strong correlation with elevation (R 2 = 0.974, p < 0.001), increasing at 15 cm per 100 m elevation, with a steep increase in snow depth at 2000-2050 m.Above 3300 m, snow depth sharply decreased at a rate of −48 cm per 100 m (R 2 = 0.830, p < 0.001).The mean "open" snowcovered area in each 1 m elevation band was 1.7 ha, with a range of 0.1 to 7.3 ha.Overall, 67 % of the study area (excluding water or developed areas) was free of canopy, including most of the 5.6 km 2 area above 3300 m.The increase in snow depth with elevation up to 3300 m is accompanied by a decrease in canopy cover with elevation.Canopy cover, based on the canopy-height model, is greater than 40 % below 2600 m, and near zero above 3200 m (Fig. 2b).

Wind and topographic effects
Hourly average wind speed at the six meteorological stations showed that the highest potential for redistributing snow was from the westerly directions, with a few periods of strong winds from the northeast at Topaz (Fig. 3).Winds were highest at the three stations above 2800 m and, to a lesser extent, at one lower-elevation station, Giant Forest, which is located in an exposed area free of upwind vegetation.Only five instantaneous gusts over 6.7 m s −1 were recorded at Panther during the snow-accumulation period, and in one instance at Wolverton; no hourly averages at these sites were over 6.7 m s −1 .
Snow depths were lowest on the southwest-and southeastfacing slopes, and highest on the northwest-and northeastfacing slopes (Fig. 2c).This pattern was most pronounced at elevations above 2400 m, and depths were low especially in the southeast between 2300 and 2700 m, which is a small fraction of the area at this elevation (Fig. 2c).The aspect with the least overall area is northeast and the greatest areal aspect representation faces northwest.
The changes in mean snow depth and slope (Fig. 2a and d), over 5-100 m averaging lengths, show an (anti)correlation at −0.16 to −0.36, with the most negative correlation at 35 m (data not shown).The most-rapid changes in slope with elevation show the least increase in snow depth; this is most evident up to 3300 m, above which the terrain becomes flatter (Fig. 2e).
The combined effects of slope and aspect express the "aspect intensity" (I A ), where higher values represent more terrain at that aspect and/or greater slope angles (Fig. 4a).This analysis reveals the slope-and aspect-feature space of the study area, where the predominant sloped aspects of north, southwest, west and northwest have positive I A values.Conversely, south, northeast, east and southeast have negative values closer to zero and are therefore less represented in the study area.At elevations < 2000 m, moderate east-and southeast-facing slopes, indicated by the slightly negative I A values, quickly rise to steeper north, northwest and west slopes, as indicated by the higher and positive I A values (Fig. 4a).Near 2400 m, southwest aspects become more predominant than north, as indicated by the crossover in I A values, and at higher elevations aspect becomes equally represented by west, southwest and northwest, with some southerly aspects (negative north I A values) above 2800 m (Fig. 4a).
To evaluate the terrain effects secondary to elevation, we applied a regression to all snow depths as a function of elevation using the slope (0.15) and intercept (−169) from the snow depth measured by LiDAR at 1850-3300 m (Fig. 4b).The residuals from this regression were then correlated with each of the predominant I A values (Table 2, Fig. 4c and d).I A snow-depth anomalies for the lowest elevations (1850-2051 m) were negatively correlated with the southeast, at the mid-elevations (2051-3301 m) most positively correlated with the northwest, and at the highest elevations (3300-3494 m) most positively correlated with the southwest (Fig. 4c and d).

Bright-band radar
The radar-sounding data include 8287 hourly observations (353 missing) from the 3 sites.While individual observations of freezing levels ranged from 200 to 2700 m, the 95th percentile values were in the range of 950-2550 m (Fig. 5).The greatest variability and highest mean freezing level occurred at the coastal station of Punta Piedras, with the lowest values at the furthest-inland station of Chowchilla.This decline in mean freezing levels going inland from the coast suggests that the snow level drops as the air mass moves inland.The third quartile of the freezing level of the Chowchilla station is 2063 m; this closely tracks the break in the coefficient of variation and correlation between snow depth and elevation observed from LiDAR at 2050 m (Figs.2a and 5), as well as the steep increase in snow depth from 1950 to 2050 m elevation (Fig. 2d).

Ground measurements
Accumulated precipitation and SWE track each other closely at the two higher-elevation sites (CRL and QUA), but measured (GNF) SWE was up to 21 cm less than total precipitation at the lowest site, showing some melt prior to the LiDAR acquisition (Fig. 6).Total precipitation at the lowest total precipitation gauges was 75 cm at GNF (2027 m) and was 72 at Lodgepole (2053 m).The LiDAR flights were 18 days after mean peak depth and three weeks before mean peak SWE (Fig. 7a, b).The mean and standard deviation (shown in paranthesis) of snow depth during LiDAR acquisition, recorded by the 26 depth sensors in the Wolverton and Panther areas, plus the 16 operational sensors co-located with snow pillows, was 210 (38) cm.This was 19 % less than the mean peak depth of 266 (44) cm recorded on 4 March.However, the mean SWE recorded by the 16 snow pillows during LiDAR acquisition was 82 (16) cm, 2 % less than the mean peak SWE of 83 (20) on 14 April.Two snow pillows, the lowest (GNF) at 2027 m and the most southerly (QUA) at 2195 m, reached peak SWE one week before acquisition, on 15 March, and had ablated 9 and 7 % SWE, respectively, prior to the time of the LiDAR acquisition (Fig. 7b).All other snow pillows either gained SWE or ablated < 5 % in the period prior to the snow LiDAR acquisition.Snow depths measured at the snow-pillow sites on the days of the LiDAR flights failed to show the elevation patterns apparent in the LiDAR depths (Fig. 8).
Daily density values calculated for the 16 snow pillows for 1 February to 30 April indicate a general trend of increasing density and consistent intrasite patterns of accumulation and densification corresponding to stormy and clear conditions (Fig. 7).Over the three month period, density decreased with each accumulation event and increased through densification as the snowpack settled, metamorphosed and integrated free water from melt or rain.At the time of the LiDAR flights the mean density was 384 kg m −3 , with a range of 83 kg m −3 and standard deviation of 42 kg m −3 , across the 1036 m elevation range represented in these data.The combined measurement error of snow-pillow and depth-sensor instruments used in the density calculation can be greater than the range of variability reported here (Johnson and Schaefer, 2002).We found low spatial variability in density that showed no significant relationship with elevation at our sites.This observation concurs with other studies of mountain snowpacks finding spatial consistency in the density of mountain snowpacks (Jonas et al., 2009;Mizukami and Perica, 2008).

Model reanalysis
The 4 km resolution PRISM data were comprised of seven grid elements in the study domain, whereas the 800 m product had approximately 4225 grid elements.Both PRISM data sets show a small upward trend in precipitation and elevation up to ∼ 300 m and a reversal of this trend at the higher elevations.The 4 km and 800 m PRISM data demonstrate similar magnitudes of increase in precipitation with elevation, 2-3 cm per 100 m, respectively.Because of this small precipitation lapse rate, the PRISM estimates diverge from the LiDAR values below about 2800 m.Total precipitation measured at two locations near the lower extent of the Li-DAR footprint during the accumulation period was 72 cm at Lodgepole (2053 m) and 75 cm at Giant Forest (2027 m)  (data available at http://cdec.water.ca.gov/, 2014) (Fig. 1).When compared with the LiDAR SWE estimates in Fig. 9, both stations show slightly more precipitation.
Total precipitation was also measured at two additional snow-pillow sensors: CRL (3170 m) and QUA (2195 m) (Fig. 1).From Fig. 6 it is apparent that snow does not account for all of the precipitation at elevations below 2200 m, but does above this elevation where rain had little influence for the accumulation period prior to the LiDAR flight.Thus the LiDAR data should reflect total precipitation above 2200 m.
The LiDAR SWE and the two reconstructed-snowmelt calculations have similar slopes of about 6 cm per 100 m (Fig. 9).The calculations from Guan et al. (2013b) most closely match the LiDAR estimates up to 3300 m, where those from Rice and Bales (2013) are offset by 20-40 cm but show a slight decrease in depth at the highest elevations.In contrast, the two PRISM precipitation models deviate from the LiDAR SWE estimates at elevations below 2800 m and have markedly different slopes.

Discussion
The overall increase in precipitation with elevation observed with airborne LiDAR is consistent with the orographic effect of mountains on precipitation (Roe, 2005 2006).Variability of the snow accumulation along the elevation gradient and deviations from a regular increase with elevation can be attributed to the interactions of topography, wind and storm tracks.Deviations from a linear increase are apparent at the lower rain-snow-transition elevations, at higher elevations near the ridge, and at intermediate elevations that have a variety of aspect and steepness characteristics.

Variability of orographic trends
Elevations over 3300 m showed the greatest negative departure from the overall orographic trend, likely due to the southwest-to-northeast-trending terrain flattening out and no longer providing the necessary lift for the same rate of adiabatic cooling (Figs. 2 and 4).Above 3300 m, the reduced lift over flatter terrain, the exhaustion of precipitable water as storms rise less steeply, and the horizontal displacement of falling snow likely all contribute to declining precipitation at the higher elevations (Mott et al., 2014;Houze Jr., 2012).These processes have been approximated in the Sierra Nevada through simulations based on the convergence of the boundary layer and slope of the local terrain but, until now, have been difficult to observe (Alpert, 1986).
However, as other researchers have noted, it is also difficult to identify the effects of specific storms on snowpack ablation due to the variability of atmospheric conditions close to the Earth's surface (Lundquist et al., 2008).The extent to which high-altitude Sierra Nevada catchments receive more precipitation than adjacent low-altitude areas varies from storm to storm, and from year to year, from occasions during which nearly equal amounts of precipitation fall at high and low altitudes to occasions when 10 or more times as much precipitation falls at the higher altitudes (Dettinger et al., 2004).In the northern Sierra Nevada, the blocking and associated terrain-parallel southerly flow of air masses, referred to as the Sierra barrier jet, can enhance lower-elevation precipitation (Neiman et al., 2008).Conversely, in the central Sierra Nevada, it has been reported that, seasonally, the ratio between higher-versus lower-elevation annual winter-season total precipitation averages about 3, but in some years the ratio drops to as low as 1 (as in 1991) or rises to as much as 4 or 5 (Dettinger, et al., 2004).While the particular orographic patterns reported here could be unique to the 2010 water year, similar patterns have been observed in the mountains of Europe, and previous works have shown consistency in the interannual spatial patterns of snow accumulation (Grünewald et al., 2014;Sturm and Wagner, 2010a;Deems et al., 2008).

Wind redistribution and radiation effects
The high spatial resolution of LiDAR snow-depth measurements points to two possible controls on snow distribution.While wind affects snow accumulation during a storm, the combined effects of wind and radiation are apparent in post-depositional changes in snow depth.As the same topographic variables influence both wind and radiation, separating the effects based on an analysis of seasonal snow-depth is challenging.
While wind patterns from a single station may be a poor indicator of the wind fields influencing snow redistribution across the entire domain, we expect snow transport by wind to be coarsely defined by the consensus of the local station's wind direction when temperatures are below zero within 24 h of a snowstorm (Figs. 1 and 3).However, the Topaz Lake station, located in smooth terrain with limited upwind influence, may best represent the wind patterns of the free atmosphere and predominant southwest storm winds.We attribute the inconsistent wind direction of other stations to the terrain-induced turbulence of the free atmosphere upwind of the stations.The M3 and Emerald Lake sites have upwind obstacles, and the Wolverton and Panther stations have low wind speeds, reflecting the muting effect of tall forest cover on wind speed and consequently snow redistribution (Fig. 3).
Consistent with prevailing winds from the southwest, we observed more accumulation on the northeast slopes and less on the southwest; however, in our domain, northeast has the least total area of all aspect quadrants, and hence these areas may be underrepresented in the analysis (Fig. 2c).
The aspect intensity variable (I A ) combines the influences of slope and aspect, and serves as a proxy for several processes affecting snow depth, e.g., radiation, upslope orographic deposition, and potentially wind and gravitational redistribution.As a result, some local anomalies, such as deepsnow-patch development, are likely masked when considering topography and snow depth as elevation-band means.
Examining residuals from a linear orographic trend by I A suggests that the steeper, northwest-facing slopes at the midelevations and northerly slopes at the lowest elevations show the greatest snow depths, likely due to the combined effects of wind deposition and lower radiation influx (Fig. 4c and  d).Conversely, low-to mid-elevation slopes prone to the combined effects of ablation and wind erosion have the least snow.These findings suggest that departures from the overall orographic trend can be observed in the elevation profile using I A , but there are limitations to the approach as used here.
It is also possible that there is limited utility in extrapolating prevailing winds from meteorological stations to predict effects of wind on snow redistribution because of the turbulence from local terrain.Research into the relationship between slope, aspect and wind has revealed that small-scale slope breaks and surface roughness have the most-significant effects on where snow accumulates locally (Li and Pomeroy, 1997b;Winstral et al., 2002;Fang and Pomeroy, 2009;Pomeroy and Li, 2000).While not part of this analysis, classification of downwind terrain has also been effective for identifying snow-patch development and persistence of localized wind deposition, offering a deterministic explanation for the spatial stationarity of snow (Winstral et al., 2002;Schirmer et al., 2011).The I A variable may also be effective for classifying locations where these processes are likely to occur.

Sublimation
Wind-driven sublimation may also play a role in the departure from the linear increase in snow depth at the higher elevations, where the highest wind velocities and thus greatest suspension of snow occur (Fig. 3).
In dry intercontinental locations, sublimation rates can be in excess of 50 %, but are much lower in the maritime climate of the Sierra Nevada and lowest during the accumulation period (Ellis et al., 2010;Essery and Pomeroy, 2001).Studies conducted at 2800 and 3100 m in the Emerald Lake basin, located in the center of our measurement domain, found net losses due to evaporation and sublimation of < 10 % for the period between 1 December and 1 April (Marks and Dozier, 1992;Marks et al., 1992).Consequently, we consider the 2010 water year cumulative loss due to sublimation and snowmelt to be limited (< 10 %) prior to the 23 March Li-DAR acquisition at all elevation bands, with more melt occurring at the lowest elevations and on the southeast-facing slopes, as indicated by the loss of SWE measured at the lowelevation snow-pillow sites and reduced snow depths on the southeast mid-elevation slopes (Figs. 2c,6,and 7).

Rain-snow transition
At lower elevations, below 2050 m, a mix of rain and snow precipitation appears to influence the amount of seasonal snow accumulation.Local SWE measurements are only available at one location below 2050 m (GNF), and this station does show a very small loss of SWE in mid-February as a result of a rain-on-snow event (Fig. 6).Nevertheless, given the expected large storm-to-storm variation in freezing level, the relatively sharp transition in slope of LiDAR-measured snow accumulation at about 2050 m suggests that most precipitation above this elevation fell as snow in the winter of 2010.
In addition, seasonal snow at the lowest elevations and on south-facing slopes has greater positive net energy exchange (from radiation or condensation), and is most susceptible to melt during the accumulation period.LiDAR snow-depth results show lower depths on south-facing versus greater depths on north-facing slopes (Fig. 2c).

Snow density
Our 23 March, calculations of snow density based on snow depth and snow pillow measurements are uncorrelated with depth or elevation but varied < 11 % from the mean, within the combined uncertainty of the sensors used to calculate them (Figs. 1 and 8).Elder et al. (1998), Anderton et al. (2004) and Anderson et al. (2014) found the variability of spring snow density to be insignificantly correlated with ele-vation in their studies, while Zhong et al. (2014) found negative correlations with elevation in a meta-study of densities in the former USSR.A range of results has also been reported for the snow-density correlation with depth showing both positive and negative correlations depending on the age of the snow and season (Arons and Colbeck, 1995;McCreight and Small, 2014).
These seemingly contradictory findings can be explained by the seasonal and climatic effects on snow depth and the snowpack energy balance and their affect on snow density.Snow depth is often positively correlated with elevation and the energy balance of the snowpack often negatively correlated with elevation; the magnitude of these effects depends on season and climate (Jonas et al., 2009;Sturm et al., 2010b).For example, in winter, when there are low levels of solar influx on low-albedo snowpacks, snow depth, which is positively correlated with elevation, has a greater influence on density.Conversely, in springtime, or in a warmer climate, a warming snowpack may reverse any previous correlation or be uncorrelated with elevation.Thus our assumption of uniform density may not be accurate for early winter but provides a reasonable estimate for spring snowpack conditions when the LiDAR snow-depth measurement was made.

Other measures of orographic trends
Although orographic precipitation is a well-documented first-order process, in the Sierra Nevada it is not well described at the watershed to basin scale owing to the very limited availability of ground-based precipitation measurements.Each set of comparative measurements used in this study provides a different index of orographic response: (i) LiDAR is a one-time snapshot of snow depth; (ii) point SWE data are small samples from highly variable spatial values; (iii) reconstructed snowmelt, or retrospective gridded SWE, reflects precipitation minus evaporation and sublimation; and (iv) PRISM is a retrospective precipitation estimate, based largely on lower-elevation stations.Nevertheless these complementary data offer spatially relevant indices of seasonally accumulated precipitation.
As Fig. 8 shows, snow depths from snow-pillow sites fail to capture the elevation patterns apparent in the LiDAR data.This pattern is also apparent in the SWE values from the same sites (Fig. 7b).While the shallowest depth is registered at the lowest elevation site (GNF, 2027 m), where a greater percentage of precipitation falls as rain, the other sites do not show a consistent increase in depth with elevation.Thus current operational measurements in the Sierra Nevada are insufficient to capture orographic trends in snow depth and precipitation.
The less steep increase in precipitation with elevation seen in the two PRISM profiles versus the LiDAR results are thought to be primarily due to the limited number of mountain stations used to calculate the PRISM trends.SWE loss from ablation and rain versus snowfall are important P. B. Kirchner et al.: LiDAR measurement of seasonal snow accumulation components of the observed LiDAR lapse rates at lower elevations, particularly below 2050 m; these processes should have only a small influence above that elevation.Evidence for this can be seen in three locations of coincident SWE and cumulative precipitation measurements (Fig. 6).The accumulated SWE and total precipitation at the two higher-elevation stations, CRL and QUA, are in close agreement, and the lowest station, GNF, shows 21 cm more total precipitation and a slight loss of SWE on 18-23 March, prior to the date of LiDAR acquisition, demonstrating that measurable rain and melt occurred at the site.In addition, precipitation station in the Kaweah Basin near the LiDAR footprint (LDG, 2053 m) had an accumulation-period total of 72 cm, higher than the LiDAR SWE estimate and lower than both PRISM estimates for the same time period (Fig. 9).The difference in annual precipitation at these sites versus annual SWE accumulation reflects in part the contribution of both rain and snow, as well as mid-winter melt, at this elevation.Thus, divergence of the PRISM and reconstructed SWE at elevations below 2200 m is expected.Temperature records in the area suggest only a small amount of winter melt at 2100 m, with very little winter melt and precipitation as rain above 2400 m (Rice and Bales, 2013).
The general pattern of SWE reconstructed from snowmelt by Guan et al. (2013b) compares with the LiDAR data, being somewhat higher at the highest elevations, lower in the midelevations, and similar at the lower elevations.Even though the reconstruction was based on energy-balance modeling, the match is somewhat surprising given the coarseness of the reconstruction model relative to the complex topography of the basin.
The Rice and Bales (2013) reconstruction, in which snowmelt was indexed to amounts and rates at the snowpillow sites, has less SWE, particularly at the mid-to lower elevations.This offset may stem in part from the higher, 106 % of average 2010 seasonal precipitation versus the lower, 90 % of average, precipitation in the 2000-2009 snowmelt reconstruction period.Further, the reconstructed SWE estimates by Rice and Bales (2013) are based on a temperature-index calculation, versus a full energy-balance approach used by Guan et al. (2013b).Also, some offset in both reconstructed SWE estimates may reflect a bias in snowcovered-area estimates, which have a 500 m spatial resolution and are heavily influenced by canopy.In other words, the LiDAR data represent open areas, and the reconstructed SWE values represent the full domain, but are empirically corrected for canopy.

Conclusions
The current results show elevation as the primary determinant of snow depth near the time of peak accumulation over 1650 m on the west slope of the southern Sierra Nevada.Li-DAR data reveal patterns potentially associated with orographic processes, mean freezing level, slope, terrain orientation and wind redistribution.Snow depth increased approximately 15 cm per 100 m elevation from snow line to about 3300 m, equivalent to approximately 6 cm SWE per 100 m elevation.This lapse rate is nearly equal to the SWEreconstruction approach, but higher than the widely used PRISM precipitation data.Localized departures from this trend of +30 to −140 cm from the kilometer-scale pattern of linear increase with elevation are seen in an elevation profile of 1 m elevation bands.Interestingly, snow depth decreased by approximately 48 cm per 100 m elevation from 3300 to 3494 m elevation.Both PRISM and SWE reconstructions show a leveling-off or reductions in SWE at higher elevations as well.
The characterization of snow depth and SWE elevation lapse rates is unique given the high accuracy and high spatial resolution of these data.Moreover, the analysis of the residuals from this elevation trend reveals the role of aspect as a controlling factor, highlighting the importance of solar radiation and wind redistribution with regard to snow distribution.While previous works have come to similar conclusions, the use of LiDAR data reveals these signals in a spatially explicit manner.As LiDAR data become more available, the analyses performed here provide a framework for evaluating the sensitivity of snow-distribution patterns to variability in location and climate.

Figure 2 .
Figure 2. (a) Snow depth (blue) with regression lines and approximate 3 December-23 March bright-band radar freezing level noted, and snow depth percent coefficient of variation (dark red); (b) percent canopy cover; (c) 35 m running average of mean snow depth and stacked area by elevation for each 90 • quadrant of aspect; (d) mean slope of 1 m elevation band; and (e) first derivative of mean slope (green) and snow depth (blue) over 35 m running average.

Figure 3 .
Figure 3. Top: hourly average wind speed and direction for 3 December-23 March accumulation period.Bottom: periods with highest probability of snow redistribution.Radius scale in m s −1 , azimuth in degrees, and north at 0 •

Figure 4 .
Figure 4. (a) Aspect intensity of LiDAR domain by elevation; (b) residuals of mean 23 March snow depth using model regression of slope, 1850-3300 m, from Fig. 2a; (c) regression of residuals for 1850-2050 m (black, dashed line) and 3301-3494 m (gray, solid line) showing departures from elevation trend for NE and SW aspect intensity; and (d) 2051-3300 m (green, dotted line) in SE and NW.

Figure 5 .
Figure 5. 3 December-23 March accumulation period hourly bright-band freezing level recorded at three wind-profiler stations upwind of the study area; locations shown in Fig. 1.Circles are 5th and 95th percentile.

Figure 6 .
Figure 6.Accumulated gauge precipitation and snow-pillow SWE for the three sites with both measurements.Locations are shown in Fig. 1.

Figure 7 .
Figure 7.In situ measurements of (a) snow depth, (b) SWE and (c) density for all west-slope snow-pillow and depth sensors in sites located within 1 • latitude of study area.Upper halves of the panels show data for individual stations, with highest and lowest elevations plotted in bold.Lower halves of the panels show mean in black, with +1 standard deviation shaded in gray; vertical blue line indicates 23 March LiDAR acquisition dates.Figure 1 shows station names and locations.

Figure 8 .
Figure 8. Observed snow depth on the 23 March LiDAR acquisition date for all west-slope snow-pillow sites equipped with depth sensors, plotted over mean LiDAR snow depth (dark gray) and 1 standard deviation (light gray).Giant Forest (GNF), Farewell Gap (FRW) and Chagoopa Plateau (CHP) are within 21 km of the measurement domain.Chilkoot Meadow (CHM) and Poison Ridge (PSR) are the sites furthest to the northwest.Locations shown in Fig. 1.

Figure 9 .
Figure 9. Elevation trend of cumulative precipitation for the Kaweah River watershed for two scales of PRISM and two reconstructions of SWE from daily snowmelt estimates, for December-March, with 23 March 2010 LiDAR SWE estimate and two cumulative precipitation gauge measurements.

Table 1 .
Target parameters and attributes for LiDAR flights.
Case Mountain, Wolverton and Panther stations are in forest openings; Emerald Lake is an alpine cirque; and Topaz and M3 are in alpine fell fields.
Wolverton and Panther -recorded maximum wind gusts at 10 s scan intervals.The M3, Topaz, and Emerald Lake stations are managed by the University of California Santa Barbara; Giant Forest is operated by the California Air Resources Board (data available at http://mesowest.utah.edu/,2014); and Case Mountain is managed by the Bureau of Land Management (data available at http://www.raws.dri.edu/,2014).The Giant Forest station is located on an exposed shrub-covered slope; the Hydrol.Earth Syst.Sci., 18, 4261-4275, 2014 www.hydrol-earth-syst-sci.net/18/4261/2014/

Table 2 .
Regression of snow-depth residuals with aspect intensity (I A ).

B. Kirchner et al.: LiDAR measurement of seasonal snow accumulation
; Roe and Baker, P.