Climate change impacts on maritime mountain snowpack in the Oregon Cascades

This study investigates the effect of projected temperature increases on maritime mountain snowpack in the McKenzie River Basin (MRB; 3041 km 2) in the Cascades Mountains of Oregon, USA. We simulated the spatial distribution of snow water equivalent (SWE) in the MRB for the period of 1989–2009 with SnowModel, a spatiallydistributed, process-based model (Liston and Elder, 2006b). Simulations were evaluated using point-based measurements of SWE, precipitation, and temperature that showed NashSutcliffe Efficiency coefficients of 0.83, 0.97, and 0.80, respectively. Spatial accuracy was shown to be 82 % using snow cover extent from the Landsat Thematic Mapper. The validated model then evaluated the interand intra-year sensitivity of basin wide snowpack to projected temperature increases (2C) and variability in precipitation ( ±10 %). Results show that a 2 C increase in temperature would shift the average date of peak snowpack 12 days earlier and decrease basin-wide volumetric snow water storage by 56 %. Snowpack between the elevations of 1000 and 2000 m is the most sensitive to increases in temperature. Upper elevations were also affected, but to a lesser degree. Temperature increases are the primary driver of diminished snowpack accumulation, however variability in precipitation produce discernible changes in the timing and volumetric storage of snowpack. The results of this study are regionally relevant as melt water from the MRB’s snowpack provides critical water supply for agriculture, ecosystems, and municipalities throughout the region especially in summer when water demand is high. While this research focused on one watershed, it serves as a case study examining the effects of climate change on maritime snow, which comprises 10 % of the Earth’s seasonal snow cover.


Significance and motivation
The maritime snowpack of the Western Cascades of the Pacific Northwest (PNW) United States is characterized by temperatures near 0 • C throughout the winter and deep snow cover that can accumulate to 5000 mm deep (Sturm et al., 1995).This important component of the hydrologic cycle stores water during the winter months (November-March) when precipitation is highest, and provides melt water that recharges aquifers and sustains streams (Dozier, 2011) during the drier months of the year (June-September).Because maritime snow accumulates and persists at temperatures close to the melting point, it is fundamentally at risk of warming temperatures (Nolin and Daly, 2006).The McKenzie River Basin (MRB, Fig. 1), located in the Central Western Cascades of Oregon, exhibits characteristics typical of many watersheds in this region, where maritime snowpack provides melt water for ecosystems, agriculture, hydropower, municipalities, and recreation -especially in summer when demand is higher and precipitation reaches a minimum (United States Army Corps of Engineers, 2001;Oregon Water Supply and Conservation Initiative, 2008).
Published by Copernicus Publications on behalf of the European Geosciences Union.

E. A. Sproles et al.: Climate change impacts on maritime mountain snowpack
In the mountains of the Western United States, snow water equivalent (SWE, the amount of water stored in the snowpack) reaches its basin-wide maximum on approximately 1 April (Serreze et al., 1999;Stewart et al., 2004).In the PNW, there have been significant declines in 1 April SWE and accompanying shifts in streamflow have been observed (Service, 2004;Barnett et al., 2005;Mote et al., 2005;Luce and Holden, 2009;Stewart, 2009;Fritze et al., 2011).This reduction in SWE has been attributed to higher winter temperatures (Knowles et al., 2006;Mote, 2006;Abatzoglou, 2011;Fritze et al., 2011).Throughout the region, current analyses and those of projected future climate change impacts show rising temperatures (Mote and Salathé, 2010).This increase is expected to transition more snow into rain, resulting in diminished snowpacks, and reduced summertime streamflow (Service, 2004;Stewart et al., 2004Stewart et al., , 2005;;Barnett et al., 2005;Mote et al., 2005;Stewart, 2009;Mote and Salathé, 2010).
This problem is not unique to the Oregon Cascades and is of significance globally as snowmelt provides a sustained source of water for over one billion people (Barnett et al., 2005;Dozier, 2011).The maritime snow class comprises roughly 10 % of the spatial extent of all terrestrial seasonal snow (Sturm et al., 1995) and includes large portions of Japan, Eastern Europe, and the western Cordillera of North America.Many of these regions are mountainous, and measurements of snowpack are limited due to complex terrain and sparse observational networks.This deficiency limits the ability to accurately predict snowpack and runoff at the basin scale, especially in a changing climate (Bales et al., 2006;Dozier, 2011).Improvements in quantifying the water storage of mountain snowpack in present and projected climates advance the ability to assess climate impacts on hydrologic processes.While climate impacts on mountain snowpack are a global concern, addressing them at the basin-level provides information at a scale that is effective for resource management strategies (Dozier, 2011).

Study area
The McKenzie River Basin has an area of 3041 km 2 and ranges in elevation from 150 m at the confluence with the Willamette River near the city of Eugene to over 3100 m at the crest of the Cascades.Precipitation increases with elevation in the MRB.Average annual precipitation ranges from approximately 1000 mm in the lower elevations to over 3500 mm in the Cascade Mountains (Jefferson et al., 2008).With winter air temperatures commonly close to 0 • C, precipitation phase is highly sensitive to temperature and can fall as rain, snow, or a rain-snow mix.In the MRB, the rain-snow transition zone is broad, ranging from 400 to 1200 m, where a transient snowpack commonly accumulates and melts over the course of a winter (Tague and Grant, 2004;Jefferson et al., 2008;Tague et al., 2008).The seasonal snow zone (areas with a distinct accumulation and ablation period) is situated above 1200 m where deep snows accumulate from November through March, increasing their water storage until the onset of melt, on approximately 1 April.In the MRB regions above 1200 m, the underlying basalt geology provides excellent aquifer storage that sustains summer flows (Tague and Grant, 2004;Jefferson et al., 2008;Tague et al., 2008;Brooks et al., 2012).Isotopic analysis found that 60-80 % of summer flow in the Willamette River originated from elevations over 1200 m in the Oregon Cascades (Brooks et al., 2012).
The MRB's "reservoir" of snow above 1200 m is especially important to the greater Willamette River Basin (30 300 km 2 ).While occupying 12 % of the Willamette, the MRB supplies nearly 25 % of the late summer discharge at its confluence with the Columbia River near Portland, Oregon (Hulse et al., 2002).Over 70 % of Oregon's population resides in the Willamette River Basin and the economy and regional ecosystems depend heavily on the Willamette River, especially in summer months when rainfall is sparse.This makes the MRB's seasonal snowpack a key resource for ecological, urban, and agricultural interests and of great interest to water resource managers in the MRB and Willamette River Basin.
Monitoring of the MRB's seasonal snowpack has been conducted for decades; however accurate measurements of basin-wide mountain snowpack do not exist.The presentday monitoring of mountain snowpack uses point-based data from the Natural Resources Conservation Service (NRCS) Snowpack Telemetry (SNOTEL) network covering an elevation range of only 245 m (1267-1512 m) in a basin where snow typically falls at elevations between 750 and 3100 m.While these middle elevations represent roughly 35 % of the basin's area, they do not quantify SWE at high elevations.
Once the snow melts at the monitoring sites, there is no further information even though snow persists at higher elevations for several weeks.In the past, this limited configuration of SNOTEL sites has functioned successfully in helping predict streamflow (Pagano et al., 2004), however the network was not designed to quantify and evaluate the impacts of projected future climate change at the watershed scale (Molotch and Bales, 2006;Brown, 2009;Nolin et al., 2012).Mean annual temperatures are projected to increase 2 • C by midcentury, potentially limiting the effectiveness of the current monitoring system.These deficiencies underscore the need for a spatially detailed understanding of snow water storage at the watershed scale, which would improve water managers' ability to manage this vital resource in the present and plan for future projected temperature changes.
Both spatially distributed snow models and remotesensing data can provide key information on spatially varying snow processes at the watershed scale, and provide diagnostic information on relationships between physiographic characteristics of watersheds and snowpack dynamics.In the past decade, spatially distributed, deterministic snowpack modeling has made significant advances (Marks et al., 1999;Lehning et al., 2006;Liston and Elder, 2006b;Bavay , 2009).Such mechanistic snowpack models also allow us to make projections for future climate scenarios.Remote sensing is an effective means of mapping the spatiotemporal character of seasonal snow (Nolin, 2011).Rittger (2012) used a computationally efficient method to compute Fractional Snow Cover Area (fSCA) from Landsat Thematic Mapper in the Sierra Nevada Mountains based on the work of Rosenthal and Dozier (1996) and Painter et al. (2009).Such data are at a spatial scale comparable to topographic and vegetation variations in the MRB and are appropriate for capturing the heterogeneous melt patterns in this watershed.By mapping fSCA, we can obtain an accurate estimate of spatially and temporally varying snow extent, however these data cannot provide estimates of SWE.
Using the MRB as a case study that is representative of mid-latitude maritime snowpacks, this research examines and quantifies the sensitivity of snowpack to climate change.Specifically the research objectives are to: (1) quantify the present-day distribution and volumetric storage of snow water equivalent at the watershed scale and across multiple decades; (2) quantify the watershed scale response of snow water equivalent to increases in temperature; and (3) quantify the watershed scale response of snow water equivalent to increases in temperature combined with increases or decreases in precipitation.

Methods
To accomplish these objectives, we applied SnowModel (Liston and Elder, 2006b) to simulate meteorological and snow conditions throughout the McKenzie River Basin at daily time steps and at a grid resolution of 100 m.The spatiallydistributed, process-based model SnowModel computes temperature, precipitation, and the full winter season evolution of SWE including accumulation, canopy interception, wind redistribution, sublimation/evaporation, and melt (Liston and Elder, 2006b).The SnowModel framework is comprised of four sub-models.MicroMet provides realistic distributions of air temperature, humidity, precipitation, temperature, wind speed and wind direction, surface pressure, incoming solar and longwave radiation (Liston and Elder, 2006a).The En-Bal sub-model computes the internal energy balance of the snowpack using atmospheric conditions computed by Mi-croMet (Iziomon et al., 2003;Liston and Elder, 2006a, b).SnowTran 3-D is a physically-based snow transport model that distributes the transport and sublimation of snow due to wind (Liston et al., 2007).SnowPack is a single layer submodel that calculates changes in snow density, depth, and SWE from fluxes in precipitation and melt (Liston and Elder, 2006b).SnowModel was selected because this study required a spatially explicit model that distributes meteorological conditions and simulates detailed calculations of the energy balance with a high degree of accuracy.Because SnowModel is physically-based it accounts for slope and aspect in calculating the energy balance, which is especially relevant in the MRB where over 30 % of the basin has slopes greater that 20 • .Additionally the role of land cover (e.g.canopy interception, sublimation, and unloading) is included in the simulations of snowpack evolution.These physical and environmental boundary conditions would be lost in a simple degree day model.
Additionally, SnowModel requires only air temperature, precipitation, relative humidity, wind speed, and wind direction data as its model forcings.These data were readily available for multiple decades.We applied data from seven automated weather stations distributed throughout the MRB at elevations ranging from 174 to 1512 m (Fig. 1, Table 1).Hypsometrically, 74 % of the area of the McKenzie River Basin is encompassed by the elevation ranges of the monitoring sites (430-1512 m), and 85 % of the basin lies below the highest elevation site of 1512 m.While higher elevation meteorological measurements would have benefitted the study, access to higher elevations was not logistically feasible.A spatially-balanced network of input stations was used to create a more evenly weighted distribution of forcing data across the watershed (Fig. 1 -stations used as model forcings are enclosed in a black square).The spatially-balanced network was found to be important in distributing precipitation (P ) and air temperature (T air ).The MicroMet sub-model uses the Barnes Objective Analysis technique, a weighted interpolation scheme based on the data spacing from a datum (station) to the grid cell (Koch et al., 1983).Clusters of stations in the center of the model domain were found to negatively impact model results in the outer regions.The addition of the Eugene Airport improved model agreement by providing a datum in the western portion of the basin.The upper elevation SNO-TEL (National Resource Conservation Service, 2012) sites were added to more evenly distribute meteorological conditions in the upper elevations.Discussion on how this configuration was finalized is discussed in greater detail in the model calibration sub-section.
The study period, WY 1989-2009, was constrained by the availability of meteorological data to drive the model, as stations were required to have a near-complete data record (greater than 90 %).A limited dataset of hourly data for meteorological stations (10 yr) was available.But because one of our primary objectives was to simulate basin-wide snowpack over multiple decades, we selected the longer daily time series.Implementing the hourly forcing data would have decreased the number of years available for the study by almost 50 % (a full decade).Additionally, the maritime snowpack of the MRB does not have a strong diurnal signal because there is little diurnal variability in air temperature.For example, we calculated the coefficient of variation (CV) for hourly air temperature in WY 2007 and found that 86 % of all days had a CV value that varied by only ±2 %.This (and other) maritime regions have snowpacks that are warm, nearly isothermal, and highly sensitive to increased temperature.These characteristics highlight the importance of studies such as this to demonstrate the accumulation and ablation sensitivities of maritime snow.Additionally the study did not focus on the sub-daily/diurnal dynamics of snowpack, which would have required hourly data.
Meteorological data was available through the study period at 00:00, 06:00, 12:00, 18:00 UTC, and with daily means of air temperature.However, only the 00:00 UTC data from SNOTEL sites are quality assessed (National Resource Conservation Service, 2013).Test iterations of the model were run with individual inputs for each of these times and results were compared to independent data for goodness of fit and Nash-Sutcliffe Efficiency (NSE) values.The data acquired at 00:00 UTC provided the best goodness of fit and NSE values.We strived to minimize model tuning so we used published values for albedo, albedo decay, and rain-snow temperature partitioning rather than use them as tuning parameters for a better fit with input data from other times of the day.Additionally, daily precipitation measurements begin and end at 00:00 UTC, which aggregated precipitation to the correct day.Our approach minimized model tuning and allowed a validated model to be run for multiple decades.
This 21 yr study period of record includes seasons with above average, normal, and below average snowpack, and years influenced by El Niño/La Niña-Southern Oscillation (ENSO) for the study period.This time period represents a warm phase of the Pacific Decadal Oscillation (Brown and Kipfmueller, 2012) and compared with records dating back 70 yr, SWE measurements are below the long-term mean (Nolin, 2012).
Physical boundary conditions for the model required elevation and land cover for the model domain, which was 112 km in the east-west direction and 76 km in the northsouth direction.Digital elevation data were obtained from the United States Geological Survey's (USGS) Seamless National Elevation Dataset (NED) (Gesch, 2007).The National Land Cover Dataset (NLCD) (Fry et al., 2009) was also obtained through USGS.The land cover boundary condition uses vegetation classes (i.e.coniferous forest, barren land), so NLCD land cover types were reclassified to the appropriate SnowModel land cover code (Sproles, 2012).Both datasets were resampled from 30 m to the model grid resolution of 100 m resolution.Resampling the 30 m data to a grid cell of 100 m captures variability in topography and snowpack across the landscape, while reducing the computational demands by a factor of eleven.Concerns over potential misclassification of land cover that may arise from reclassification are moderated by landscape patterns in the areas where snowfall occurs.These areas are almost entirely coniferous forests in the Western Cascades or unforested, exposed landscapes in the High Cascades.Any misclassification in resampling would most likely only occur at transitional areas.A greater concern regarding land cover is the application of a static land cover dataset over a 21 yr period in a region with a dynamic forest landscape that includes active timber harvest and re-planting.However, developing a dynamic land cover dataset lies outside the scope of this research.
The overall goals of providing spatial and temporal estimates of basin-wide SWE across multiple decades were completed in four general steps: (1) apply a physically based, spatially distributed model that uses meteorological data as model forcings; (2) calibrate and validate model outputs of P and T air using independent station data; (3) calibrate and validate model outputs of SWE using station data and maps of snow covered area from remote sensing; and (4) conduct a sensitivity analysis of snowpack with regard to temperature and precipitation.Each of these steps is described in greater detail below.

Model modifications
Two primary modifications were made to SnowModel: a rain/snow precipitation partition function and an albedo decay function.These modifications more accurately simulate physical conditions, and improved model performance.The rain/snow precipitation partition function was required because in maritime climates, wintertime temperatures commonly remain close to 0 • C and mixed phase precipitation events are common.In the PNW, empirical measurements by the United States Army Corps of Engineers (USACE) (1956) show that the transition from rain to snow exists primarily between a temperature range of −2 to 2 • C. Based upon the USACE study the relationship was implemented in the model using Eq. ( 1).SFE = (0.25 * (275.16− T air )) * P (1) where, SFE (Snow Fall Equivalent) is the amount of amount of precipitation reaching the ground that falls as snow, T air is air temperature in degrees Kelvin, and P is total precipitation.Rainfall is computed as P minus SFE.We tested more computationally complex rain-snow algorithms and results were virtually identical.The USACE linear partition provided higher computational efficiency so we proceeded with this approach.
The shortwave albedo of snow (α) has significant effects on surface energy balance, internal energetics, and seasonal evolution of snowpack (Wiscombe and Warren, 1980).Previous versions of SnowModel included snow albedo as a static, tunable parameter (Liston and Elder, 2006b).This study applied an improved snow albedo decay function from Strack et al. (2004) where: for non-melting conditions and, for melting snow Where α t is the snow albedo value used at each time step by the model in energy balance calculations, α t−1 represents the snow albedo at the previous time step, and the decay gradient is represented by gr m = 0.018, gr nm = 0.008 for melting and non-melting conditions, respectively.The maximum albedo value after new snowfall (when new snow depth ≥ 2.5 cm) is set to 0.8 in unforested areas and to 0.6 in forested areas (Burles and Boon, 2011).A minimum snow albedo (α min ) was set to 0.5 in unforested areas and 0.2 in forested areas.
We understand that applying a single albedo decay function has its limitations, and does not account for variation in land cover or topographic effects (Molotch et al., 2004).This potential source of model error is addressed in the Discussion.

Model calibration
Model calibration was comprised of two phases that carefully examined the accumulation and the ablation periods.Lapse rate (discussed in next section), the rain-snow partition (Eq.1), and maximum snow density (330 kg m −3 ) were the parameters tuned during calibration.The configuration of meteorological stations also played an important role in model calibration, and is discussed later in this section.We applied a systematic approach, adjusting a single parameter representative of published values, ensuring that changes in model outputs were a result of modifications to the individual parameter.Outputs were qualitatively and quantitatively evaluated until the model was calibrated.The initial phase focused on optimizing the spatiallydistributed gridded values of daily P and T air .Because meteorological conditions are first order controls on snowpack accumulation and ablation, maximizing the accuracy of these spatially interpolated and temporally varying model forcings is an important first step.Without accurate inputs, the resulting snowpack might be calibrated to correct values, but not for the right reasons (Kirchner, 2006).The second phase focused on optimizing simulations of snowpack.Model evaluation used point-based measurements for SWE and Landsat fSCA remote-sensing data for snow cover extent, providing a robust means of model calibration and validation (Andersen and Bates, 2001).Prior to the implementation of the albedo decay function and rain-snow partition, there was an overestimation of modeled snow extent compared to the point-based measurements and remote-sensing data.However, once these modifications were incorporated into the model, spatial agreement improved considerably.This improvement makes sense conceptually.The fixed rain-snow partition simulated 100 % of precipitation to fall as snow when air temperature was 2 • C or colder, and lead to an overestimation of snow.Compounding this overestimation was a fixed albedo that underestimated shortwave energy critical to the melt process.The rain-snow partition would proportion less precipitation falling as snow, and the albedo decay would hasten the melt process.
The optimal configuration of meteorological stations was determined by iteratively adding individual stations in the model.Results of each iteration were compared to stations independent of those used in the model (Table 1) using metrics described later in the next section.Paired sets of water years with statistically high, low, and average peak SWE were used to calibrate and validate the model (Table 2).Calibration was performed on the first set of water years, and then validated to the second set of water years.Once model calibration and validation was completed for the selected years, the model was run for WY 1989-2009 to establish a presentday reference simulation for applying the future climate projections, and hereafter is referred to as the study period.

Calibration metrics
Nash-Sutcliffe Efficiency (NSE) and Root Mean Square Error (RMSE) were used to evaluate modeled P , T air , and SWE compared to measured values from SNOTEL stations and meteorological stations independent of those used in the model.NSE is a dimensionless indicator of model performance where NSE = 1 when simulations are a perfect match with observations.For 0 < NSE < 1, the model is more accurate than the mean of the observations.While an NSE values > 0.50 are considered satisfactory (Moriasi et al., 2007), we used a target threshold of 0.80 or greater for all stations.This value represents a model efficiency that is very close to measured values and is significantly better than using mean values (Nash and Sutcliffe, 1970;Legates and McCabe, 1999).If NSE is less than 0, the mean is a better predictor (Nash and Sutcliffe, 1970;Legates and McCabe, 1999).RMSE indicates the overall difference between observed and simulated values, and retains the unit of measure (Armstrong and Collopy, 1992).RMSE provided a better understanding of the scale of error that occurred in simulations, and was used as a metric to improve model results.
Air temperature proved to be a challenging parameter to calibrate due to the complex terrain of the MRB.Here, true temperature lapse rates do not always follow a linear temperature-elevation relationship and synoptic scale atmospheric patterns can affect local lapse rates, especially when high pressure systems dominate causing cold air pooling (Daly et al., 2010).For the model, we used initial monthly lapse rates from the Washington Cascades, roughly 350 km north of the MRB (Minder et al., 2010).These lapse rates were iteratively adjusted to minimize RMSE for temperature using the forcing and evaluation stations listed in Table 1.The final model iteration applied monthly lapse rate values ranging from 5.5-7 • C km −1 and were 1.5 • C km −1 cooler than Minder found in the Washington Cascades (Table 3).Minimum RMSE for some calibration sites were outside of the target threshold of 2 • C, as large errors for a few values can exacerbate RMSE values (Freedman et al., 1991).Thus R 2 values (Legates and McCabe, 1999) and 95 % confidence intervals were calculated (Freedman et al., 1991) to augment model evaluation.R 2 values describe the proportion (0.0 to 1.0) of how much of the observed data can be described by the model, and confidence intervals indicate simulation reliability.Methods on how to potentially improve lapse rate calculations for future work are described in the third paragraph of the Discussion section.
Field measurements of SWE acquired during WY 2008 and 2009 were used to augment model calibration.We measured SWE manually at five sites within the basin (Figs. 1  and 4) from December to July during WY 2009 on approximately the first day of each month.Snow densities were calculated using monthly SWE measurements at five locations in the basin.Four snow depth measurements were conducted within one meter of the initial SWE sample.This approach does not provide a detailed measurement of SWE in a 100 m × 100 m grid cell, and thus was used as a broad metric for assessing the magnitude of simulated SWE and the timing of accumulation and ablation.Logistically, this rapid assessment approach allowed samples at all five sites to be conducted in a single day.In addition, colleagues at the University of Idaho provided SWE measurements at two locations in the basin on two dates in WY 2008and 2009(Link et al., 2010).

Remote sensing based calibration
The spatial extent of modeled snow cover was assessed using satellite-derived maps of fractional snow-covered area (fSCA).The Landsat TM fractional snow covered area data were aggregated from 30 m data to the 100 m grid resolution of SnowModel and converted to a binary grid where < 15 % fSCA was classified as no snow, and > 15 % fSCA was classified as snow in the grid cell.The co-occurrence of modeled and measured snow cover was assessed using metrics of accuracy, precision, and recall as in Painter et al. (2009).Precision is the probability that a pixel identified with snow indeed has snow.Recall, the metric that Dong and Peters-Lidard (2010) employed, is the probability of detection of a snow-covered pixel.Accuracy is the probability a pixel is correctly classified.For detailed explanations of these measures and their application to snow mapping, see Rittger (2012).There were a limited number of valid images each winter because of cloud cover common in maritime climates and the 16 day repeat cycle of Landsat.For example, during WY 2009, only one image between the months of November and April had a cloud cover less than 25 % in the MRB.However, each calibration year had at least one image with cloud cover less than 10 % that could be used to effectively assess the spatial accuracy of the model.While the day of year of Landsat acquisition varied across years, multiple images were acquired during accumulation, peak, and ablation phases of SWE.The spatial agreement between fSCA and SnowModel results was evaluated for physiographic variables including land cover class, elevation, slope and aspect.This allows us to identify the physical characteristics of the domain that were potentially misrepresented by the model.

Climate perturbations
The calibrated and validated model was run for the study period and then used to assess the sensitivity of snowpack to increased temperature and variable precipitation.To determine the response of snowpack to increased temperature and changes in precipitation, a sensitivity analysis was conducted in three phases.The first phase increased all temperature inputs for WY 1989-2009 by 2 • C (hereafter referred to by T2), which is considered to be the mean annual average temperature increase in the region by mid-century (Mote and Salathé, 2010).The second and third phases retained the temperature increases, but also scaled precipitation inputs by ±10 % to incorporate the uncertainty in projected future precipitation (Mote and Salathé, 2010).Hereafter these phases will be referred to by T2P10 (representing +2 • C and a 10 % increase in precipitation), and T2N10 (representing +2 • C and a 10 % decrease in precipitation).Results from the ±10 % precipitation also provide insight into how annual variability in precipitation can affect SWE relative to the effects of increased temperature.The model was then run, applying the three sets of scaled meteorological data for the study period of WY 1989-2009.

Model assessment
Model results were evaluated at SNOTEL stations, meteorological stations in the HJA, and our field measurement sites (Figs.2-4, Table 4).Model simulations of P and T air performed well at input stations (used to force the model) and reference stations (used to validate the model) (Fig. 2a and  b).For years other than calibration and validation years, the mean NSE of P and T air at all stations was 0.97 and 0.80, respectively (Table 4) during the snow season (1 November-30 June).The model simulations of SWE (Figs. 3 and 4) showed mean NSE coefficients of 0.83 across the basin at automated SNOTEL locations and 0.70 at the field sites.Spatially, simulations had an overall accuracy of 82 % compared to the Landsat fSCA data.
WY 1997 and 2005 were excluded from these metrics and in subsequent calculations discussed in this section.Evaluation of model results showed two unrelated problems for these years.WY 1997 experienced at least 10 winter precipitation events > 50 mm day −1 .Evaluation of the input data showed that in a few cases there were significant discrepancies (> 1 m of annual cumulative precipitation) at several of the stations that were used as forcing data.Additionally, a few large precipitation inputs were offset by one day.The shifts were not systematic and appeared to be random in nature, most likely due to equipment mistiming at several stations.As a result storms with a significant amount of total precipitation (> 50 mm) were, in effect, double counted and processed on two consecutive days by the model.While the errors were present in less than 10 % of the datasets these events were characterized by heavy precipitation and cold temperatures that increased snowpack accumulation.This double count of precipitation provided simulations with around a 1 m overestimation of SWE, roughly the same magnitude as the over estimation of annual precipitation.Thus this year was omitted.WY 2005 displayed model deficiencies in resolving lapse rates associated with temperature inversions.Simulations of spatially distributed gridded temperature in WY 2005 had an RMSE of 3.8 • C and NSE of 0.72, whereas the study period had values of 2.5 • C and 0.80, respectively.This was due to extended periods of high pressure, which resulted in cold air pooling and negative temperature lapse rates (Daly et al., 2010).Extensive snowmelt and near complete loss of upper elevation snowpack occurred in mid-to-late February (National Resource Conservation Service, personal communication, 2009; National Resource Conservation Service, 2012) as unseasonably warm temperatures at higher elevations and unseasonably cool temperatures at lower elevations persisted for several weeks.The model deficiencies caused by such extensive temperature inversions are addressed in the Discussion section.Precipitation was effectively distributed for all stations and across the full range of elevations used in the validation (Fig. 2a).The mean RMSE error was 0.01 m and the mean NSE value was 0.96 for the full study period.It is important to note that the addition of the low elevation Eugene Airport meteorological station (174 m) greatly improved model performance.This station provided meteorological input data at a low elevation and at the western edge of the model domain, which improved the spatial interpolation of precipitation.
Air temperature had a mean RMSE of 2.5 • C and mean NSE value of 0.80 (Fig. 2b and Table 4).Model simulations at the Santiam Junction SNOTEL station consistently underperformed in relation to all other stations.Santiam Junction is adjacent to a state highway, an Oregon Department of Transportation facility, and an airstrip which combined, make it more exposed to wind than the nearby natural forest setting found at the other stations.The station elevation also provided a small bias, as simulations at middle elevation stations (800-1300 m) underestimated T air on average by 2.0 • C. The upper elevation stations (1300-1550 m) overestimated temperature on average by 0.25 • C.This bias reflects the topographic character of the MRB.The upper elevation sites are situated in the High Cascades geological province, where the topography has a more gradual slope averaging approximately 10 • .In the Western Cascades (up to 1300 m) geological province, slopes are steeper averaging approximately 20 • , but are also frequently characterized by slopes up to 50 • .In the Western Cascades during periods of high pressure, it is common to have cold air drainage, where cooler, more dense air moves down a slope and pools in valleys creating cooler temperatures at lower elevations (Daly et al., 2010).
The RMSE for T air (2.5 • C) was larger than anticipated, however further analysis showed an R 2 of 0.85 and 98 % of all T air simulations within a 95 % confidence interval.The additional evaluation metrics support the likelihood that a small minority of poor model simulations for T air had a significant impact on RMSE.Efforts in calibrating and evaluating temperature suggest that the standard approach of applying linear monthly lapse rates to temperatures would contribute to the underperformance found in this study.Ideas on how to  (Table 4).The length and consistency of the automated SWE data record at lower elevation sites is more limited.With the exception of UPL, snow pillows in the HJA are not calibrated and the reported data have not been fully quality assured.The result is an inconsistent dataset with values that often do not represent expected snowpack evolution in the region.Due to the questionable accuracy of the measured SWE values in the HJA, these data were not used as a metric for model validation.This issue also highlights the need for a careful calibration and regular maintenance of SWE measurement sites.
In the spatial validation, 14 yr of SnowModel simulations of snow cover compared to Landsat TM fSCA (converted to snow/no snow) had an overall accuracy of 82 % (the ratio of correctly identified grid cells -i.e.snow as snow, bare as bare), and overall precision of 71 % (the probability that a pixel identified with snow indeed has snow) and an overall recall of 93 % (the proportion of positives correctly identi-fied as positives).Although the accuracy statistic may rise because of overwhelming numbers of cells in which there is no snow (Rittger et al., 2012), we include it because a large portion of the MRB can be snow covered and validation scenes are distributed throughout the season.Disagreement between the fSCA images and simulations primarily occurred where the model estimated snow cover and the fSCA did not have snow cover (13 %).This degree of False Positive (FP) is expected as remotely sensed data typically omits snow cover in the steep and heavily forested landscapes that dominate the Western Cascades and the MRB (Nolin, 2011).The interannual changes associated with harvested forest are not expressed in the static land cover dataset, but are incorporated into the fSCA product.This classification discrepancy propagated through each year contributing to the lower precision value by decreasing the number of True Positive (TP).Additionally, the fSCA binary product classifies any cell with a fractional snow cover value less than 15 % as no snow.Even though the Landsat fSCA product was coarsened to 100 m, cells at the transitional snow line will be classified as no snow and result in an increase in False Positive (FP) classifications for modeled snow cover.WY 2006WY , 2008WY , and 2009 were the exceptions, showing more False Negative (FN) classifications, but with a similarly higher level of agreement.For a more detailed discussion of the model assessment using remote-sensing data, please refer to Sproles (2012).

Sensitivity of snowpack to changes in temperature and precipitation
The response of snowpack in the MRB in the T2 scenario highlights the sensitivity to temperatures and that the greatest impact on SWE accumulation comes from more precipitation falling as rain rather than snow.Elevations below 1300 m show a substantial loss of SWE accumulation (Fig. 3), where elevations around 1500 m suggest considerable losses of SWE, but still retain a seasonal snowpack with distinct accumulation and ablation periods.Mean peak SWE for the basin (the ±5 day mean from peak SWE) decreased by an average of 56 % for the study period (Figs. 5 and 6ad, Table 5).When integrated over the area of the MRB, this equals an annual average loss of 0.70 km 3 of water stored as snow -equivalent to 230 mm of SWE distributed across the basin, and is more than twice the volume the largest impoundment in the MRB (Cougar Reservoir, storage capacity 0.27 km 3 ).While temperature is the controlling factor for the phase of precipitation and in turn changes in SWE, changes in total precipitation also have an impact.The T2P10 and T2N10 scenarios show losses of mean area-integrated peak SWE of 0.62 to 0.78 km 3 (203 to 256 mm of SWE distributed across the basin, respectively), and reflect the role that precipitation variability plays on peak snowpack in the  MRB.The 0.21 km 3 (69 mm) difference of area-integrated peak SWE predicted by the T2P10 and T2N10 scenarios is substantial and is equal to slightly less than available storage at Cougar Reservoir.However, 2 • C temperature increases alone result in a 0.70 km 3 loss (230 mm of SWE distributed across the basin, Figs.6a-d and 7, Table 5).Increased precipitation in the T2P10 scenario results in additional SWE at elevations primarily over 1800 m, where a 2 • C increase in temperature is not sufficient to convert snowfall to rainfall or to significantly accelerate snowmelt.However, this increase in SWE at the high elevations only partially offsets some of the losses at lower elevations.
With warmer conditions, the date of peak SWE is projected to occur earlier in the spring and properly into the winter (before the vernal equinox).The average date for simulated peak SWE in the MRB during the study period is 31 March.However, in T2 the average date for peak SWE shifts 12 days earlier in the WY.Similarly, peak SWE arrives 6 days and 22 days earlier in the T2P10 and T2N10 scenarios, respectively, indicating a greater sensitivity in the T2N10 than the T2P10 scenario.
We assessed the sensitivity of the snowpack to temperature increases by elevation using the 10 day mean of peak SWE and frequency of snow cover for WY 2007.The 10 day mean of peak SWE minimized the influence of any single large accumulation event in order to emphasize the overall snowpack trend for that season.WY 2007 was a statistically average year for SWE at the four SNOTEL sites.Peak SWE was −0.07 m of the reference mean and had a standard deviation of 0.02 m from the reference mean value (0.83 m).In WY 2007 the greatest net losses of peak SWE were found between 1001 and 1500 m (Fig. 8).This elevation zone generated 53 % of the basin-wide losses of SWE in the T2 scenario, and comprises 45 % of the basin area.Proportionately, the areas between 1501 and 2000 m generate a more significant component of peak SWE loss.This elevation zone generated 45 % of the basin-wide peak SWE losses in the T2 scenario, but comprises only 17 % of the basin area.The mean loss of peak SWE lost per grid cell was 0.61 m in this elevation zone, as compared to 0.26 m in areas between 1001 and 1500 m.
The duration of snow cover by grid cell was assessed for WY 2007 during the accumulation and melt period between 1 January to 30 September 2007.As expected, the snow cover frequency in the T2 scenario was lower across the basin, with the areas between 1001 and 1500 m affected the most.This range of elevations saw an average of 36 fewer days of snow cover than in the reference year (Fig. 8).Elevations between ∼ 1501 and 2000 m see a less dramatic reduction of snow covered days.Areas between ∼ 2001 and 2500 m experienced increased losses in snow cover days with elevation.

Discussion
Our results quantified the basin-wide distribution and volumetric storage of snow water in the MRB, which averaged 1.26 km 3 of SWE (414 mm distributed across the basin) over the study period.This natural "reservoir" stores roughly five times more water than the largest impoundment in the watershed.The maritime snowpack of the MRB was highly sensitive to increased temperatures, showing a 56 % loss in peak SWE when temperature forcings were increased by 2 • C. Projected warmer conditions also hasten the melt cycle, with peak SWE occurring 12 days earlier.Elevations between 1000 and 2000 m are most affected in the T2 scenario as snow transitions to rain, and snow on the ground has an enhanced melt cycle (Fig. 3). Figure 6c-d suggest that a 2 • C temperature increase will shift snowpack characteristics by approximately 250 m.The elevation zone from 1000-1500 m has the greatest volumetric loss of stored water (Figs.6a-d and  8), and represents the largest areal proportion of the basin.In WY 2009, mean SWE values for this elevation band diminish considerably, from 0.69 to 0.13 m.Elevations above 2000 m are affected by warmer temperatures, but to a lesser extent, retaining a similar average SWE by elevation (Fig. 6c-d).
The ±10 % change in precipitation inputs explores how variability in precipitation affects snowpack.A 10 % decrease in precipitation exacerbates the impacts of temperature on snowpack, especially for the elevation zone from 1000-2000 m.A 10 % increase in precipitation only slightly buffers the loss of peak SWE.A notable result of the 10 % increase in precipitation identifies the elevations that are less sensitive to increased temperature.For instance, peak SWE 198919901991199219931994199519961998199920002001200220032004200620072008  increases in the T2P10 scenario above ∼ 2000 m where increased precipitation also increases the seasonal accumulation of SWE.However even with gains at high elevations, there is still a considerable net loss of snowpack (−49 %) compared to the study period.Not surprisingly, the response of snow cover frequency to a 2 • C increase is very similar to the pattern of the change in SWE (Fig. 8).Snow cover duration in the elevation zone from 1000-1500 m were most affected, with some locations losing more than 80 days of snow cover in an average snow year.
Initially the meandering nature of the snow loss curves in Fig. 8 might not seem intuitive, but can be explained by the topography of the MRB.Elevations between ∼ 1001 and 1500 m can receive both rain and snow during the winter, even though elevations above 1200 m retain a seasonal snowpack.This elevation range is the most sensitive to increased temperature and shows a transition to a rain-dominated area with a 2 • C increase.Elevations between ∼ 1501 and 2000 m are less sensitive to increased temperatures and more likely to retain enough precipitation falling as snow with a 2 • C increase to develop distinct periods of accumulation and ablation.Retention of the snowpack in this elevation range is aided by the highly-dissected Western Cascades (which dominate this elevation) where adjacent terrain provides shade, reduces incoming short-wave radiation, and mitigates potential snow loss (DeWalle and Rango, 2008).This shading also helps explain the loss of snow between ∼ 2000 and 2500 m, where topography shifts from the rugged Western Cascades to the more exposed High Cascades.This shift towards a gradual, consistent slope in the High Cascades provides less shading throughout the course of day that would potentially mitigate increased temperatures.
Efforts in calibrating and validating the model clearly demonstrated that precipitation and temperature are first order controls on snowpack accumulation and peak SWE.This highlights that it is critical to achieve optimal accuracy of the spatially distributed values of P and T air prior to calibrating the model based on SWE.P had a high level of agreement between observations and simulations (NSE of 0.97).There were distinct similarities between the R 2 (0.85) and NSE of T air (0.80) with the NSE of SWE (0.83) and the accuracy of the spatial distribution of snowpack (82 %).These similarities lead to the logical conclusion that improvements in accuracy of snowpack simulations can be made through improvements in temperature simulations.
The challenges in simulating T air are partially explained by the physical characteristics of the MRB.Daly et al. (2010) used empirical data to establish that expected temperature lapse rates that exist between elevation and temperature are often decoupled from one another and are largely controlled by topography and elevation.Steeper slopes can produce cold air drainage and different lapse rates than lapse rates for more gentle slopes (Daly et al., 2010).Additionally, moisture content of a storm (as determined by its temperature, source area, and history) affects the wet adiabatic lapse rate.Daly et al. (2010) suggest that variability in lapse rates may increase with projected future climate.Combined, these factors highlight the shortcomings of using a standard temperature lapse rate in a model.Though outside of the scope of this research, an improvement to the monthly static lapse rates used in SnowModel would be to compute dynamic lapse rates with a dual-pass approach.The first pass through the meteorological station data would establish the lapse rate and the second pass would apply time step specific lapse rates in the Barnes Objective Analysis method to spatially distribute temperature data.This would allow an individual storm's lapse rate characteristics to be included in the model.A dynamic lapse rate would also help during stable conditions when cold air drainage may be important.
The high level of agreement for P was attained once an evenly distributed network of input stations was established.In initial model runs, incorporating multiple clustered While this study achieved a high level of agreement between simulated and measured values, the complex topography and land cover of the MRB also introduce potential sources of error.While the MicroMet and EnBal sub-models implicitly included land cover and topography in calculating incoming shortwave radiation, the albedo function did not account for these factors in reflected shortwave radiation.An improved albedo function inclusive of land cover and topography would provide the opportunity to reduce or account for model error.Similarly, there are limitations in estimating snow cover extent in complex terrain using remotely sensed images (Rittger et al., 2013).The thick vegetation of the MRB potentially obscures snow underneath the forest canopy.Our validation applied the most recent scientific advances in calculating fSCA in mountainous regions (Rittger et al., 2013), that can be improved in future work by field validation.
The elevation range of stations (174-1512 m) limited model assessment at higher elevations.Hypsometrically, this comprises 74 % of the basin, and extends from regions dominated by rain into the seasonal snow zone above 1200 m.The authors recognize that higher elevation measurements would have benefitted the study, and could help minimize uncertainty in future research.Unfortunately access to elevations above 2000 m was not logistically feasible during the field season.Future improvements in field methods would include at least one high elevation site above 1800 m rather than three lower elevation sites.Data was also limited temporally, with hourly data beginning in WY 1999.Because one of the primary goals of the study was to simulate snowpack for multiple decades, we moved forward pragmatically, applying the best data available supplemented by field observations.We are applying our findings to improve field measurements in the basin that will ultimately aid future model-based studies in this watershed and region.

Looking forward -the impacts of climate perturbations on snowpack
Losses in SWE and declining snow duration will impact years with high, low and average snowpack and will change the statistical representation and human perceptions of what a high, low and average snowpack represents.The MRB will increasingly experience more precipitation falling as rain rather than snow in warmer conditions.Areas presently in the rain/snow transition zone will become dominated almost entirely by rain.The changes will affect the timing and magnitude of runoff during the winter, spring, and summer months as more precipitation shifts from snow to rains (Stewart et al., 2005;Jefferson et al., 2008;Jefferson, 2011).
While research has shown that geology controls baseflow in sub-basins of the MRB, (Tague and Grant, 2004;Jefferson et al., 2008;Tague et al., 2008), shifts in the form of precipitation will affect the timing and magnitude of peak runoff.These shifts will be seen at the basin and sub-basin scale, potentially influencing water resource managers' decisionmaking process.The moderately high spatial and temporal resolutions of the simulations allow the sensitivity of diminished snowpack to be evaluated for the MRB and its subbasins.This range of scales provides the ability to develop potential adaptive water resource management strategies.For instance, dam operators now release flow in anticipation of runoff generated by snowmelt.But these results suggest that sub-basins with headwaters in the elevation zone from 1500-2000 m will see dramatic losses in SWE and lose the ability to store winter precipitation as snow.As the contribution from snowmelt decreases and more runoff shifts to earlier in the year, dam operations will need to reflect these changes in their management strategy.Results from this study have already helped water resource professionals choose a site for a new SNOTEL station to augment the existing monitoring network (Webb, personal communication, 2011) and develop water management strategies for municipal water use (Morgenstern, personal communication, 2010).
Snow and snowmelt serve as a resource for winter and summer recreation, agriculture, industry, municipalities, and hydropower.The difference with a 2 • C increase in temperature on peak area-integrated SWE is considerable (0.70 km 3 or 230 mm of basin-wide SWE) -more than twice the size of the largest impoundment in the basin.While this estimated loss only pertains to the MRB it would scale up to be major factor at the regional level.Potential management concerns pertaining to the supply of water could be compounded by shifts in the demand of water as well.Oregon's population is expected to grow by 400 000 from 2010 to 2020 (Office of Economic Analysis, 2011).The increase in population would most likely increase demand especially in the summer and fall when stakeholders compete for an already limited supply (United States Army Corps of Engineers, 2001;Oregon Water Supply and Conservation Initiative, 2008).Because mountain snowpack serves as an efficient and cost-effective reservoir, any research that examines socio-economic topics should contain a mountain snowpack component.For example, an examination of socio-economic impacts of the adaption costs associated with mitigating climate change would need to include the costs associated with a diminished mountain snowpack.

Conclusions
This research provided the first detailed spatial and temporal understanding of snow accumulation and ablation in the MRB for present conditions and serves as a prognostic tool for understanding snowpack in projected future climates.Because maritime snow accumulates at temperatures close to 0 • C, the seasonal accumulation and ablation of maritime snow is sensitive to temperature.These findings provide insights into the mechanisms controlling snowpacks in such environments and serves as an example of the magnitude and types of changes that may affect similar watersheds in a warmer climate.Moreover, with the modifications made to the model (rain-snow partitioning, albedo decay function), this model can readily be transitioned to other regions with maritime snow with minimal reconfiguration.Although this study focused on a single watershed, the processes affecting snowpack in the McKenzie River are similar to other maritime snowpacks across the Earth.
Mountain snowpack is a key common-pool resource, providing a natural reservoir that supplies water for drinking, worship, hydropower, agriculture, ecosystems, industry, and recreation for over 1 billion people globally.The spatial dis-tribution of maritime snowpack and its sensitivity to climate change at basin scale does not provide global answers, but it does provide clarity at a scale appropriate for developing management strategies for the future (Seibert and McDonnell, 2002).

Fig. 1 .
Fig. 1.Context map for the McKenzie River Basin, Oregon.Model forcing locations are enclosed by a black square.

Fig. 2 .Fig. 3 .
Fig. 2. Model performance for precipitation (top -(a)) and temperature (bottom -(b)) in years not used in calibration and validation.Each dot represents a day during 1 November-30 June.

Fig. 4 .
Fig. 4. Model performance of SWE at field locations in WY 2009.Location of field sites are shown in Fig. 1.
Field measurements collected at a range of elevations during WY 2008 and 2009 also show a high level of agreement between measured and modeled SWE values with an NSE coefficient of 0.70 (Fig.4).These field sites suggest that model results successfully simulate the timing and magnitude of snowpack evolution, especially at the higher elevations.The lower elevation sites (S1-S3) show a lower level of agreement during the ablation period.Here, spring SWE is less than 0.1 m and the mean difference between the observed and simulated values during the ablation period was only 0.07 m.It is worth noting the highest SNOTEL site is situated at an elevation of 1512 m, but 75 % of the modelestimated SWE lies above that elevation.This result is consistent with the work ofGillan et al. (2010) who found that > 70 % of SWE accumulates above the mean elevation surrounding SNOTEL sites in a snow-dominated watershed in Northwestern Montana.

Fig. 6 .
Fig. 6.Map of simulated SWE on 1 April 2009 with a 2 • C increase in temperature (a).Map of the loss of simulated SWE on 1 April 2009 with a 2 • C increase in temperature (b).SWE by elevation on 1 April 2009 using reference observed data (c).SWE by elevation on 1 April 2009, 2009 with a 2 • C increase in temperature (d).Each dot on the plot represents a grid cell in the MRB.The values in parentheses represent the mean SWE and standard deviation of SWE (meanSWE, stdSWE) at each elevation band.The upper elevations are not affected as significantly as the lower elevation snowpack.

Fig. 7 .
Fig. 7. Peak SWE integrated over the area of the MRB and its sensitivity to a 2 • C increase in temperature.

Fig. 8 .
Fig. 8. Loss of SWE (upper) and snow covered days (lower) by elevation with a 2 • C increase on 1 April 2007.Each dor on the plot represents a grid cell in the MRB.Snowpack between 1000 and 2000 m are the most sensitive to temperature and show the greatest losses.

Table 1 .
Meteorological and snow monitoring stations that were applied as model forcings and/or in evaluation of simulation results.T air -Air Temperature, P -Precipitation, RH -Relative humidity, Wind -Wind speed and direction, SWE -Snow water equivalent; NWS -National Weather Service, HJA LTER -HJ Andrews Long Term Ecological Research site, NRCS -National Resource Conservation Service.

Table 2 .
Water years used in the calibration and validation of the model.Selected Values in parentheses represent the deviation from the mean (in meters) of peak SWE measurements at Santiam Junction, Hogg Pass, Roaring River, and McKenzie.Years noted by an * represent years with field measurements of SWE.

Table 3 .
Lapse rate values ( • C km −1 ) used in SnowModel and those published by Minder et al.The values posted by Minder et al. (2010) are for the Washington Cascades, which are approximately 350 km north of the MRB.Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec

Table 4 .
Mean Nash Sutcliffe Efficiency (NSE)Rating and Root Mean Squared Error for Daily SWE, and T and Annual P .These stations all have 10 or years of record.

Table 5 .
Changes in peak SWE, % of peak SWE lost, and the shift in the number of days earlier for the MRB averaged across the reference period.