Winter Snow Level Rise in the Northern Sierra Nevada from 2008 to 2017

The partitioning of precipitation into frozen and liquid components influences snow-derived water resources and flood hazards in mountain environments. We used a 915-MHz Doppler radar wind profiler upstream of the northern Sierra Nevada to estimate the hourly elevation where snow melts to rain, or the snow level, during winter (December–February) precipitation events spanning water years (WY) 2008–2017. During this ten-year period, a Mann-Kendall test indicated a significant (p < 0.001) positive trend in snow level with a Thiel-Sen slope of 72 m year−1. We estimated total precipitation falling as snow (snow fraction) between WY1951 and 2017 using nine daily mid-elevation (1200–2000 m) climate stations and two hourly stations spanning WY2008–2017. The climate-station-based snow fraction estimates agreed well with snow-level radar values (R2 = 0.95, p < 0.01), indicating that snow fractions represent a reasonable method to estimate changes in frozen precipitation. Snow fraction significantly (p < 0.001) declined during WY2008–2017 at a rate of 0.035 (3.5%) year−1. Single-point correlations between detrended snow fraction and sea-surface temperatures (SST) suggested that positive SST anomalies along the California coast favor liquid phase precipitation during winter. Reanalysis-derived integrated moisture transported upstream of the northern Sierra Nevada was negatively correlated with snow fraction (R2 = 0.90, p < 0.01), with atmospheric rivers representing the likely circulation mechanism producing low-snow-fraction storms.


Introduction
As the climate warms, the partitioning of snowfall to rainfall in snow-dominated mountain watersheds is likely to change [1][2][3][4][5]. Changes in precipitation phase alter seasonal snowpack dynamics, ecological processes, peak streamflow timing, and winter flood hazards [5][6][7]. These changes present different challenges for water supply forecasting and reservoir operations [8], particularly in states that depend upon snow-derived water resources and are at risk for winter flooding [9,10]. Continued climate change, increased drought risk [11], and more intense extreme precipitation events [12] will necessitate adaptive water management strategies.
During the past seven water years (WY; 1 October-30 September), California and Nevada experienced the full range of hydroclimate extremes. Precipitation reductions during a multiyear drought from WY2012 to WY2015 were found to be consistent with average conditions that were estimated during the extreme and persistent Medieval droughts [13], with anomalous positive temperature anomalies enhancing drought severity [14]. Conversely, WY2017 was the wettest year in the past century and followed a near-average WY2016. These wide-ranging hydroclimatic conditions provide incentive to study the physical mechanisms that may be influencing hydrologic variability in the Sierra Nevada, as this knowledge can be used to inform and prioritize adaptation strategies for sustainable water management under a changing climate [15].
This study focuses on the northern Sierra Nevada, a 150 km wide north-south trending mountain range with crest elevations ranging from 2000 to 3000 m ( Figure 1). Approximately 50% of annual precipitation occurs during the winter months (December-February; hereafter winter) [16]. The orientation of the northern Sierra Nevada orthogonal to prevailing westerly winds creates significant orographic precipitation effects that influence precipitation magnitude and spatial distribution (e.g., [17] and references therein). Water stored as snow that accumulates during the cool season (November-March) subsequently melts and runs off (typically during April-July) and provides the primary source of streamflow in the American, Yuba, Truckee, and Feather River basins ( Figure 1). All of these basins are managed for multiple purposes, including flood control, and to meet ecological and human consumptive demands. The relatively low elevation of these basins makes them susceptible to changes in precipitation phase due to warming [18], which has significant implications for winter flooding [9,19] and warm season water availability [15]. Because of the importance of the northern Sierra Nevada for water resources and flood prevention in California and Nevada, a novel network of hydrometeorological measurements exist upstream of the Sierra Nevada [20] and offer a unique means to explore hydroclimatic change in this region.  Snow level, or the elevation where snow melts to rain, is a primary control on the sensitivity of mountain snowpack accumulation to climate warming during precipitation events [21]. The brightband elevation, or the altitude of maximum radar reflectivity in the melting layer, provides an adequate estimate of the snow-level height and thus the phase of precipitation for given elevations [22,23]. We employed measurements of brightband elevations using a network of snow-level sensing radars established for real-time hydrometeorological monitoring in California to evaluate observed snow-level changes during the past 10 winters in the northern Sierra Nevada (Figure 1) in order to address the following objectives: (1) Test the hypothesis that snow levels have been rising during the past decade; (2) See if consistency exists between the changes in inferred precipitation phase deduced from snow-level radar observations and those estimated by an empirical relationship at weather stations developed by Dai [24]; and, (3) Explore plausible physical mechanisms controlling precipitation phase using reanalysis products and observational data, with the hypothesis that atmospheric rivers will strongly control snow fraction and snow level.

Snow Level Radar
We estimated snow levels on the basis of brightband heights measured at hourly intervals from December 2007 through February 2017 from a 915-MHz Doppler radar wind profiler [25], located in Chico, California ( Figure 1). These measurements were complemented by brightband elevations that were estimated by frequency-modulated, continuous-wave, snow-level radars [26] at Oroville and Colfax, California (Figure 1), spanning the meteorological winters from December 2010 (December 2011 for Oroville) to February 2017. The data for all three radars were obtained from the Earth Systems Research Laboratory [27]. The radar data and algorithm used to convert Doppler velocity and radar reflectivity into melting elevation are described in White et al. [28,29].
The brightband height data from the snow-level sensing radars as well as the precipitation and temperature observations from the weather stations (described in Section 2.2) were examined to identify possible outliers. These included snow levels above 4000 m and temperatures above 30 • C or below −15 • C, while hourly precipitation measurements greater than 50 mm were filtered to eliminate erroneous measurements. We applied a diurnal correction to account for over-and under-estimation to radar-derived snow levels [22]. We shifted snow-level values during daylight hours (9:00 A.M.-5:00 P.M. LST) upward 100 m for Chico and 20 m for Oroville and Colfax. We shifted snow-level values lower during nighttime hours (200 m for Chico and 150 m for Oroville and Colfax). We also applied a latitudinal correction [22] to Colfax consisting of a 20 m increase to facilitate comparisons with Oroville and Chico. Our trend analysis results did not change significantly with or without applying these corrections because these changes were systematically applied. The snow-level sensing radars are located upstream of the Sierra crest (Figure 1), therefore the brightband heights should be interpreted as the maximum values given potential mesoscale lowering of the snowline [21,30].

Station-Based Observations
We obtained hourly temperature and precipitation values from alter-shielded weighing gauges operated by the California Department of Water Resources [31] from two mid-elevation stations: Mineral, California (1511 m), and Blue Canyon (1609 m) ( Figure 1). These data span the winters of 2008 through 2017 and were included to compare snow-fraction estimates made at high temporal resolution against snow-level radar observations. To examine the relationship between lower temporal resolution observations and snow-level radar observations as well as extend the snow fraction estimates back further in time, we also used daily minimum and maximum values of temperature and precipitation from nine National Weather Service Cooperative Observer (COOP) [32] stations with >85% complete records [33]. Data from these stations, that are in the northern Sierra Nevada snow belt (>1200 m; Figure 1), spanned December 1950 through February 2017. In order to develop a relationship between observed precipitable water and snow levels, hourly GPS-observed precipitable water data at Petaluma, California was acquired from SuomiNet [34]. We compared precipitable water data only during hours (n = 1249) when snow-level observations were present for six hours during a calendar day. The reported correlation of determination (R 2 = 0.6) was maximized with a three-hour lag.
The fraction of precipitation falling as snow (hereafter, snow fraction) was estimated based upon the approach developed by [24] and employed by [35]. We summed hourly precipitation totals at 0.5 • C intervals from −15 • C through 25 • C with respective observations from Mineral and Blue Canyon. We calculated conditional frequencies of snow at each 0.5 • C interval using the hyperbolic tangent function developed by [24] but with parameters estimated for the Sierra Nevada ecoregion by [35]. The conditional probabilities can be interpreted as likelihoods of snow at given temperatures [22]. We multiplied the probability of snow by the total precipitation at each temperature interval. We then divided the sum of this value by the total precipitation to calculate the snow fraction for each water year. This calculation was repeated for the nine COOP stations at daily time steps for the WY1951-2017 period but with maximum temperatures [35]. We report the median snow fraction calculated as an average of median values from the nine stations.

Tests for Trend
To test the hypothesis that non-zero trends in snow level and snow fraction existed during the past decade (between WY2008-2017), we used the Mann-Kendall [36,37] test that was modified to account for serial correlation [38] and reported the Thiel-Sen slope [39] for this period.

Large-Scale Controls on Snow Level and Snow Fraction
We used 1 • gridded monthly sea surface temperatures (SST) from the COBE-SST reanalysis [40] and 2.5 • gridded output of daily precipitable water and 700 hPa geopotential heights from the National Center for Environmental Prediction (NCEP) reanalysis [41] from December 1950 through February 2017. We used integrated vapor (moisture) transport at a grid point that was located upstream of the northern Sierra Nevada from the NASA Modern Era Retrospective Analysis (MERRA) [42] and a catalog of atmospheric rivers developed using the atmospheric river detection method of Rutz et al. [43] for the winters spanning WY1981-2017.
To explore a plausible, large-scale, physical mechanism favoring precipitation-phase changes between snow and rain, we performed single-point correlations [16] to calculate Pearson's correlation coefficient (r) between each SST grid point (averaged over winter) and the time series of the nine-station COOP median snow fractions from WY1951 to 2017. We removed linear trends from each time series prior to performing correlations [44]. We only noted grid points as statistically significant when the p-value of the correlations satisfied the false discovery rate test [45] at the alpha = 0.05 level.
To further investigate an atmospheric circulation mechanism driving observed high and low snow fraction events, we used the Tahoe City COOP (Figure 1) to identify precipitation days exceeding the 90th percentile of non-zero precipitation days. These days were subset into warm (>3 • C) and cold (<−2 • C) events. We selected Tahoe City because it is near many important river basins in our study area and because the observed snow-fraction variability from this station reasonably approximated the nine-station median (R 2 = 0.67, p < 0.01). Under the two sets of temperatures, the snow-fraction calculation [24] implies that the cold days will be dominated by frozen precipitation, whereas warm days will result in the melting of the existing snowpack [22]. For the respective wet warm and cold combinations, we composited NCEP precipitable water and 700 hPa geopotential height fields over the eastern North Pacific domain. Next, to provide a continuous examination of the role of moisture transport on snow fraction, we correlated MERRA-derived integrated vapor transport with daily snow fraction bins from 0 to >0.95, with a bin size equal to 0.05 (or 5%). The number of atmospheric rivers identified by the method of [43] was counted for each bin.

Snow Levels and Snow Fractions
Snow levels at Chico were variable across years during the observed record ( Figure 2). The lowest median value was observed during WY2009, and the highest values were observed during the peak drought years of WY2014-2015 ( Figure 2). A significant non-zero trend (p < 0.001) during WY2008-2017 was identified with a Thiel-Sen slope of 72 m year −1 ( Figure 2). This trend was broadly consistent with the WY2011-2017 slopes at Oroville and Colfax (approximately 50 m year −1 ). Observed snow levels at Chico and Colfax were positively correlated (R 2 = 0.84, p < 0.01); Figure 3). Snow fractions at Blue Canyon and Mineral were moderately well correlated (R 2 =0.57, p < 0.05) and varied from as high as 0.55 to as low as 0.19. The correlation between Chico snow level and Mineral snow fraction was Snow levels at Chico were variable across years during the observed record ( Figure 2). The lowest median value was observed during WY2009, and the highest values were observed during the peak drought years of WY2014-2015 ( Figure 2). A significant non-zero trend (p < 0.001) during WY2008-2017 was identified with a Thiel-Sen slope of 72 m year −1 ( Figure 2). This trend was broadly consistent with the WY2011-2017 slopes at Oroville and Colfax (approximately 50 m year −1 ). Observed snow levels at Chico and Colfax were positively correlated (R 2 = 0.84, p < 0.01); Figure 3). Snow fractions at Blue Canyon and Mineral were moderately well correlated (R 2 =0.57, p < 0.05) and varied from as high as 0.55 to as low as 0. 19   Consistent with the rise in snow level shown by snow-level radar, we observed a negative trend in snow fraction at both Blue Canyon and Mineral (Figure 2), but only the slope at Mineral was Consistent with the rise in snow level shown by snow-level radar, we observed a negative trend in snow fraction at both Blue Canyon and Mineral (Figure 2), but only the slope at Mineral was significant (p < 0.001). With the exception of WY2009-2010, snow fractions at Blue Canyon and Mineral were well correlated, and varied from as high as 0.55 to as low as 0.19. Snow fraction declined at a faster rate at Mineral (−0.033 year −1 or −3.3% year −1 ) than at Blue Canyon (−1.4% year −1 ). Assuming an average environmental lapse rate over all land of 5.1 • C per 1000 m [24] and on the basis of the modified equation for precipitation phase [24,35], the 450 m increase in average median snow levels between the five-year periods spanning WY2008-2012 and WY2013-2017 ( Figure 2) equates to a warming of approximately 2.3 • C. At Mineral, 50% of precipitation fell on days below 1 • C. Therefore, an increase of 2.3 • C reduced the snow fraction by about 25%. This value falls between the linear estimate of snow-fraction change during the 10-year period (33%) and the difference in average snow fraction between these periods (15%). The observed mean annual winter temperature during hours when precipitation was recorded at Mineral rose by 1.6 • C between these periods. This is lower but broadly consistent with the estimated warming from median snow levels that were found using an average environmental lapse rate.  Consistent with the rise in snow level shown by snow-level radar, we observed a negative trend in snow fraction at both Blue Canyon and Mineral (Figure 2), but only the slope at Mineral was Between WY2008 and WY2017, the COOP snow fraction ( Figure 4) and median Chico snow level were negatively correlated (R 2 = 0.95, p < 0.001; calculated using Spearman's rank correlation due to the shorter time period). The past ten years of COOP data demonstrated decreases in snow fractions ( Figure 5) that are consistent with the increasing snow levels observed by snow-level radar and decreasing snow fractions at hourly stations ( Figure 2). Of the 67 years analyzed, the median snow fraction varied substantially among years from more than 0.4 in WY1971 to less than 0.1 in WY2015 ( Figure 5). Wet winters can have low (e.g., WY1956, WY1997, WY2017) or high (e.g., WY1969 and WY1993) snow fractions. Dry winters also can have either low (e.g., WY1991 and WY2015) or high (e.g., WY 1965 and WY1985) snow fractions. These results are consistent with [46], with no significant correlation between winter precipitation and snow fraction. The snow fraction slope of −3.5% year −1 (p < 0.0001; Figure 5) during WY2008-2017 was the steepest negative slope of any observed in the 10-year period from WY1951-2017, and its Thiel-Sen slope was consistent with the declines estimated for Mineral (−0.035 versus −0.033). . This value falls between the linear estimate of snow-fraction change during the 10-year period (33%) and the difference in average snow fraction between these periods (15%). The observed mean annual winter temperature during hours when precipitation was recorded at Mineral rose by 1.6 °C between these periods. This is lower but broadly consistent with the estimated warming from median snow levels that were found using an average environmental lapse rate. Between WY2008 and WY2017, the COOP snow fraction ( Figure 4) and median Chico snow level were negatively correlated (R 2 = 0.95, p < 0.001; calculated using Spearman's rank correlation due to the shorter time period). The past ten years of COOP data demonstrated decreases in snow fractions ( Figure 5) that are consistent with the increasing snow levels observed by snow-level radar and decreasing snow fractions at hourly stations ( Figure 2). Of the 67 years analyzed, the median snow fraction varied substantially among years from more than 0.4 in WY1971 to less than 0.1 in WY2015 ( Figure 5). Wet winters can have low (e.g., WY1956, WY1997, WY2017) or high (e.g., WY1969 and WY1993) snow fractions. Dry winters also can have either low (e.g., WY1991 and WY2015) or high (e.g., WY 1965 and WY1985) snow fractions. These results are consistent with [46], with no significant correlation between winter precipitation and snow fraction. The snow fraction slope of −3.5% year −1 (p < 0.0001; Figure 5) during WY2008-2017 was the steepest negative slope of any observed in the 10year period from WY1951-2017, and its Thiel-Sen slope was consistent with the declines estimated for Mineral (−0.035 versus −0.033).

Seasonal and Event Relationships Between Large-Scale Circulations, Snow Level, and Snow Fraction
At the seasonal timescale, single-point correlations between the time series of COOP-derived snow fractions and SSTs broadly resembled an expression of the Pacific Decadal Oscillation (PDO; Figure 6a) [47] and are negatively correlated (R 2 = 0.13, p < 0.01) with the winter PDO index [48]. No correlation existed (R 2 = 0.01, p = 0.354) between September-November multivariate El Niño-Southern Oscillation (ENSO) index [49] and subsequent winter snow fraction. The only significant negative correlations between snow fraction and SST anomalies (R 2 > 0.25, p < 0.05) were found immediately offshore of California. Following [16], we produced a time series by averaging SST anomalies in this region for comparison against snow-fraction anomalies. The resulting time series illustrates the generally negative correlation between SSTs and snow fractions (Figure 6b; note reversal of right yaxis) with occasional exceptions (e.g., WYs 1973 and 2013). The correlations were the strongest when no lag was applied (i.e., using December-February SSTs instead of November-January SSTs). Positive SST anomalies offshore of California were associated with reductions in observed winter snow fractions.
Composites of precipitable water and 700 hPa geopotential height for warm, wet days ( Figure  6c) at Tahoe City demonstrated a coherent, southwest-northeast oriented precipitable water plume in the southeastern subtropical Pacific that exceeded 25 mm with positive precipitable water anomalies (calculated as differences between composited events and mean 1981-2010 winter precipitable water) exceeding 8 mm (Figure 6e). The shape and orientation of this composite moisture plume resembled an atmospheric river [50]. This agreed with the negative correlation between MERRA-derived integrated vapor transport (R 2 = 0.90, p < 0.01; calculated using Spearman's rank correlation) and the finding that the four lowest snow fraction bins (0-0.2) had nearly twice the number of identified atmospheric rivers as all the other bins combined (302 versus 170; Figure 7a). This was also consistent with the positive relationship between GPS-observed precipitable water at the coast and Chico snow level (R 2 = 0.6, p < 0.01; Figure 7b). Warm and wet days corresponded with zonal flow, a negatively tilted trough in the northeastern Pacific, and an omega block over the Bering Sea (Figure 6c). These features are similar to composites of leeside flood-producing storms with high snow levels [9]. The cold and wet days had a weaker composite signal of subtropical moisture (<20 mm precipitable water) between Hawaii and California (Figure 6d). Precipitable water anomalies ranged from +2 mm offshore of southern California to −3 mm near 30 °N, 140 °W with a broad area of negative precipitable water anomalies (−4 to −8 mm) in the eastern Gulf of Alaska and northeastern Pacific (Figure 6f). These composite days demonstrated a meridional flow with a deeper negatively tilted northeastern Pacific trough and an amplified ridge centered near 160 °W (Figure 6d).

Seasonal and Event Relationships Between Large-Scale Circulations, Snow Level, and Snow Fraction
At the seasonal timescale, single-point correlations between the time series of COOP-derived snow fractions and SSTs broadly resembled an expression of the Pacific Decadal Oscillation (PDO; Figure 6a) [47] and are negatively correlated (R 2 = 0.13, p < 0.01) with the winter PDO index [48]. No correlation existed (R 2 = 0.01, p = 0.354) between September-November multivariate El Niño-Southern Oscillation (ENSO) index [49] and subsequent winter snow fraction. The only significant negative correlations between snow fraction and SST anomalies (R 2 > 0.25, p < 0.05) were found immediately offshore of California. Following [16], we produced a time series by averaging SST anomalies in this region for comparison against snow-fraction anomalies. The resulting time series illustrates the generally negative correlation between SSTs and snow fractions (Figure 6b; note reversal of right y-axis) with occasional exceptions (e.g., WYs 1973 and 2013). The correlations were the strongest when no lag was applied (i.e., using December-February SSTs instead of November-January SSTs). Positive SST anomalies offshore of California were associated with reductions in observed winter snow fractions.
Composites of precipitable water and 700 hPa geopotential height for warm, wet days (Figure 6c) at Tahoe City demonstrated a coherent, southwest-northeast oriented precipitable water plume in the southeastern subtropical Pacific that exceeded 25 mm with positive precipitable water anomalies (calculated as differences between composited events and mean 1981-2010 winter precipitable water) exceeding 8 mm (Figure 6e). The shape and orientation of this composite moisture plume resembled an atmospheric river [50]. This agreed with the negative correlation between MERRA-derived integrated vapor transport (R 2 = 0.90, p < 0.01; calculated using Spearman's rank correlation) and the finding that the four lowest snow fraction bins (0-0.2) had nearly twice the number of identified atmospheric rivers as all the other bins combined (302 versus 170; Figure 7a). This was also consistent with the positive relationship between GPS-observed precipitable water at the coast and Chico snow level (R 2 = 0.6, p < 0.01; Figure 7b). Warm and wet days corresponded with zonal flow, a negatively tilted trough in the northeastern Pacific, and an omega block over the Bering Sea (Figure 6c). These features are similar to composites of leeside flood-producing storms with high snow levels [9]. The cold and wet days had a weaker composite signal of subtropical moisture (<20 mm precipitable water) between Hawaii and California (Figure 6d). Precipitable water anomalies ranged from +2 mm offshore of southern California to −3 mm near 30 • N, 140 • W with a broad area of negative precipitable water anomalies (−4 to −8 mm) in the eastern Gulf of Alaska and northeastern Pacific (Figure 6f). These composite days demonstrated a meridional flow with a deeper negatively tilted northeastern Pacific trough and an amplified ridge centered near 160 • W (Figure 6d). correlations were statistically significant (p < 0.05). (b) Time series of northern Sierra Nevada median snow fraction anomaly (right y-axis) and average SST anomalies (left y-axis) within the gold box shown in (a). Composites of National Center for Environmental Prediction (NCEP) precipitable water (filled contours) and 700 hPa geopotential heights for 90 th percentile winter precipitation events at Tahoe City (Figure 1) on warm days (maximum temperatures > 3 • C) (c) and cold days (maximum temperatures < −2 • C). (d) Composite precipitable water (contours, interval 2.5 mm) and precipitable water anomalies (filled contours, in mm) calculated as differences between the (e) warm/wet composite days and (f) cold/wet composite days and the long-term mean 1981-2010 winter precipitable water output from the NCEP reanalysis.
Water 2017, 9, 899 9 of 14 median snow fraction anomaly (right y-axis) and average SST anomalies (left y-axis) within the gold box shown in (a). Composites of National Center for Environmental Prediction (NCEP) precipitable water (filled contours) and 700 hPa geopotential heights for 90 th percentile winter precipitation events at Tahoe City ( Figure 1) on warm days (maximum temperatures > 3 °C) (c) and cold days (maximum temperatures < −2 °C). (d). Composite precipitable water (contours, interval 2.5 mm) and precipitable water anomalies (filled contours, in mm) calculated as differences between the (e) warm/wet composite days and (f) cold/wet composite days and the long-term mean 1981-2010 winter precipitable water output from the NCEP reanalysis.

Discussion
We used snow-level radar to identify a trend in median snow levels between WY2008 and 2017 in the northern Sierra Nevada (Figure 2) with an estimated annual rate of 72 m year −1 . This rise coincides with decreases in snow fractions derived from both hourly ( Figure 2) and daily ( Figure 5) weather station observations. The consistency between changes in inferred snow fraction at stations (Figures 2 and 5) and those that were estimated by applying an average environmental lapse rate to snow-level changes ( Figure 2) gives us confidence that the snow-level radars are accurately recording a robust change in precipitation phases, hence supporting our first objective. The significant positive correlation (R 2 = 0.95, p < 0.01) between snow levels and snow fractions provides confidence that the snow-level radars are accurately recording precipitation phases. The correlation also affords confidence that the empirically-based estimate of snow fraction represents a reasonable means to estimate past or future changes in precipitation phase (e.g., [51]), thus supporting our second objective. These conclusions are despite uncertainties that are posed by storm-scale variability [22], potential bias in snow-fraction calculations because of precipitation gauge collection efficiency [52,53], and errors in snow-level observations. The ability to estimate changes in precipitation phase operates under the assumption that the empirical relationship to calculate snow fraction [15] remains stationary in time and is not sensitive to rain and snow transitions or spatiotemporal variability [46]. Further improvement in the snow-fraction calculations could be obtained by developing site-specific parameter values [46] and using additional observations [5].
Note that much of our study period corresponded with a period of notable hydroclimatic extremes. A severe drought during WY2012-2015 [12,13] was followed by subsequent average and record wet years during WY2016 and WY2017, respectively. The 67-year calculation of nine-station

Discussion
We used snow-level radar to identify a trend in median snow levels between WY2008 and 2017 in the northern Sierra Nevada (Figure 2) with an estimated annual rate of 72 m year −1 . This rise coincides with decreases in snow fractions derived from both hourly ( Figure 2) and daily ( Figure 5) weather station observations. The consistency between changes in inferred snow fraction at stations (Figures 2  and 5) and those that were estimated by applying an average environmental lapse rate to snow-level changes ( Figure 2) gives us confidence that the snow-level radars are accurately recording a robust change in precipitation phases, hence supporting our first objective. The significant positive correlation (R 2 = 0.95, p < 0.01) between snow levels and snow fractions provides confidence that the snow-level radars are accurately recording precipitation phases. The correlation also affords confidence that the empirically-based estimate of snow fraction represents a reasonable means to estimate past or future changes in precipitation phase (e.g., [51]), thus supporting our second objective. These conclusions are despite uncertainties that are posed by storm-scale variability [22], potential bias in snow-fraction calculations because of precipitation gauge collection efficiency [52,53], and errors in snow-level observations. The ability to estimate changes in precipitation phase operates under the assumption that the empirical relationship to calculate snow fraction [15] remains stationary in time and is not sensitive to rain and snow transitions or spatiotemporal variability [46]. Further improvement in the snow-fraction calculations could be obtained by developing site-specific parameter values [46] and using additional observations [5].
Note that much of our study period corresponded with a period of notable hydroclimatic extremes. A severe drought during WY2012-2015 [12,13] was followed by subsequent average and record wet years during WY2016 and WY2017, respectively. The 67-year calculation of nine-station median snow fractions indicated that high or low snow fractions ( Figure 5), and thus snow levels ( Figure 2) could occur during either wet or dry winters, as consistent with [42]. We speculate that dry, low-snow-fraction years (e.g., WY2015) would have a more severe impact on hydrologic systems than a dry, high-snow-fraction year (e.g., WY1989). Continuing work seeks to more closely examine streamflow responses to winters with varying combinations of precipitation total and precipitation phase distributions. These years may provide analogs to possible future climate regimes that can inform water management strategies [14].
Our third objective was to evaluate the role of circulation (atmospheric rivers) and regional climate conditions (SSTs) in influencing snow fraction and snow levels. Our composite analysis indicating that low-snow-fraction storms have stronger (anomalous positive) moisture plumes (Figure 6c,e) is consistent with rain-on-snow climatologies [9,54]. Further evidence is provided by the negatively correlated (R 2 = 0.90, p < 0.001) relationship between mean-maximum-upstream, MERRA-derived vapor transport and Tahoe City snow fraction (Figure 7a). The positive relationship (R 2 = 0.6, p < 0.01) between observed precipitable water near the California coast and the three-hour lagged snow level at Colfax (Figure 7b) supports the notion that storms with greater precipitable water ( Figure 6c) and moisture transport (Figure 7a) are associated with higher snow levels. This relationship is consistent with [55] although our correlation is weaker, perhaps because of the distance between Colfax and Petaluma. Overall, these results support the hypothesis that atmospheric rivers are important drivers of snow level and snow fraction variability at the event scale.
As found in field studies [56] and modeling experiments [57], anomalous positive SSTs offshore of continents modify air masses by contributing upward heat flux and decreasing static stability, leading to increases in downstream precipitation. At longer (seasonal) timescales, we observe that the negative correlations between snow fraction and SSTs offshore of California (Figure 6a,b) are located near the primary region of storm track activity that is correlated with winter California precipitation [16]. The warming effect of SSTs on the warm sector of the cyclone may enhance baroclinicity [58] and subsequently increase low-level winds and moisture fluxes (Figure 7a) that promote greater orographic precipitation [59]. Atmospheric rivers making landfall in California during warm coastal SST regimes are associated with stronger southerly winds and enhanced poleward moisture transport [60]. During warm (low snow fraction) storms, the large-scale zonal flow promotes the advection of warm, moist subtropical air into the southwestern U.S. (Figure 6c,e), with atmospheric rivers largely contributing to low-snow-fraction storms (Figure 7a).
Additional alternative hypotheses of the mechanisms producing snow level rise are the subject of ongoing study. Coupled ocean-atmosphere modeling approaches could help disentangle the relationship between wetter and warmer storms that result from increases in saturation vapor pressure due to temperature rise and those undergoing modifications from ocean heat fluxes (e.g., [61]) or poleward advection of subtropical airmasses [60]. Additional research is also necessary to better understand the relative roles of tropical-extratropical interactions (e.g., ENSO and the Madden-Julian Oscillation [62,63]) that influence winter snow level and snow fraction at timescales varying from seasonal to individual events.
Because of the short duration of analyzed data, we recommend further evaluation of observations to examine whether the apparent upward-snow-level trend continues. This represents a primary limitation of our study. We have provided evidence that the snow-level radar and snow-fraction estimates appear to offer reasonable metrics to observe and evaluate precipitation phase change in snow-dominated watersheds [1,5,[28][29][30]. These changes will have marked effects on warm season streamflow [4], especially during drought years [13]-with negative implications for urban, ecological, and agricultural water demands [13][14][15]64]. If such a connection between SSTs, moisture plume strength, and snow fractions does exist, it suggests that snow accumulation may further decline with continued regional warming. This warming is expected to produce wetter, stronger atmospheric rivers [12], that in turn will contribute to shallower snowpacks and a greater potential for rain-on-snow events and flooding [9,10,54,55]. These changes also may limit drought recovery during occasional wetter years that occur during persistent droughts [65]. Warming background temperatures combined with changes from snow to rain leads to decreased water availability in spring and throughout the warm season [4,[66][67][68]. Our findings thus suggest another mechanism that may act in tandem with background warming to amplify the loss of snow-derived water resources in the western U.S. [69,70]. This may be amplified if northeast Pacific SST variability and persistence increase, thus leading to longer duration episodes of anomalous warm SSTs [47] that support decreased snow fractions (Figure 6a,b).

Conclusions
We identified a statistically significant positive (negative) trend in winter snow levels (snow fractions) in the northern Sierra Nevada during the winters between WY2008 and 2017. We found consistency between increases in the elevation of winter snow levels measured by snow-level sensing radars and estimated snow fractions. Atmospheric rivers are predominantly associated with lowsnow-fraction storms, and anomalously warm coastal SSTs appear to favor lower-snow-fraction winters-perhaps by enhancing the upward heat flux during atmospheric river landfalls and promoting zonal flow regimes that increase the advection of warm subtropical air into the northern Sierra Nevada, leading to storms with higher snow levels. This hypothesis will require modeling experiments to test its validity. If true, it suggests that continued increases in sea-surface temperatures and increased frequencies in atmospheric river landfalls may exacerbate future snowpack decline in the Sierra Nevada. Such changes will have negative implications for water availability. The hypothesis also emphasizes the importance of maintaining and expanding hydroclimatic monitoring in mountain environments [5,20,71].