Uncertainty in seasonal snow reconstruction: Relative impacts of model forcing and image availability

Advances in Water Resources 55 (2013) 165–177

Contents lists available at SciVerse ScienceDirect

Advances in Water Resources journal homepage: www.elsevier.com/locate/advwatres

Uncertainty in seasonal snow reconstruction: Relative impacts of model forcing and image availability A.G. Slater a,⇑, A.P. Barrett a, M.P. Clark b, J.D. Lundquist c, M.S. Raleigh c a

Cooperative Institute for Research in Environmental Sciences, University of Colorado, Boulder, CO, USA National Center for Atmospheric Research, Boulder, CO, USA c Department of Civil and Environmental Engineering, University of Washington, Seattle, WA, USA b

a r t i c l e

i n f o

Article history: Available online 31 July 2012 Keywords: Snow Reconstruction Uncertainty Solar radiation MODIS Reanalysis

a b s t r a c t There are many areas of uncertainty when solving the inverse problems of snow water equivalent (SWE) reconstruction. These include (i) the ability to infer the Final Date of the Seasonal Snow (FDSS) cover, particularly from remote sensing; (ii) errors in model forcing data (such as air temperature or radiation fluxes); and (iii) weaknesses in the snow model used for the reconstruction, associated with both the fidelity of the equations used to simulate snow processes (structural uncertainty) and the parameter values selected for use in the model equations. We investigate the trade-offs among these sources of uncertainty using 10,000 station-years worth of data from the western US SNOTEL network. Model structural and parameter uncertainty are eliminated by using a perfect model scenario i.e. comparing results to modelled control runs. The model was calibrated for each station-year to ensure that the model simulations reflect reality. Results indicate that for a temperature index model, a ±5 days error in FDSS gives a median 25%/+32% error in maximum SWE. A 1 °C air temperature bias produces a SWE error larger than a 5 days error in the FDSS for 50% of the 10,000 cases. Similarly, a 5 days error in FDSS could be accounted for by a net radiation error of 13 W m2 or less during the melt period, in 50% of cases. Mean absolute errors of 1 °C or more are typically reported in the literature for air temperature interpolations at high elevations. Observed solar radiation during the melt season can differ by 30 W m2 over relatively short distances, while estimates from reanalysis (NARR, ERA-Interim, MERRA, CFSRR) and GOES satellites typically span more than 40 W m2. Using data from both MODIS sensors (Terra & Aqua) at all snow covered points in the western US, a consecutive 5 days gap in imagery at time of FDSS is likely to occur only 5–10% of the time. This work shows that errors in model forcing data are at least as important, if not more, than image availability when reconstructing SWE. Ó 2012 Elsevier Ltd. All rights reserved.

1. Introduction Over much of the western US snow is the dominant term in shaping the annual hydrograph. Retrospective knowledge of basin-wide snow water equivalent (SWE) can be valuable for a number of reasons, such as planning assessments, precipitation estimation, model evaluation or climatological studies. Unfortunately, the ability to close the water balance over all but the smallest basins via direct measurement is still the bane of hydrology; accurate monitoring of stores and fluxes remains a dastardly difficult prospect [34]. In turn, the advent of higher resolution visible/infrared satellite imagery has led to the suggested use of inverse methods to obtain spatially distributed snow mass [38]. ⇑ Corresponding author. Address: National Snow and Ice Data Center, CIRES, University of Colorado at Boulder, Boulder, CO 80309-0449, USA. Tel.: +1 303 735 5358; fax: +1 303 492 2468. E-mail address: [email protected] (A.G. Slater). 0309-1708/$ - see front matter Ó 2012 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.advwatres.2012.07.006

Snow water equivalent reconstruction is the process of determining the amount of snow that existed in a basin or at a point by back-calculating the amount of melt that occurred prior to the complete removal of snow [11]. One motivator of SWE reconstruction is the assumption that knowledge of melt energy fluxes would be superior to knowledge of precipitation accumulation. In order to perform a reconstruction of maximum SWE, several pieces of information are required, as listed below, and each of these will contain its own uncertainty: (1) Knowledge of the final date of seasonal snow cover (FDSS) for any given point is required to ensure that potential melt is summed over the correct period. This item is akin to initial conditions for the inversion process. When considering large basins, practicality usually obliges that this information be derived from satellite imagery (e.g. AVHRR, MODIS or LandSat). However, imagery can be subject to differing return periods (e.g. MODIS is twice daily, LandSat is every 16 days),

166

A.G. Slater et al. / Advances in Water Resources 55 (2013) 165–177

and the underlying scene can often be cloud obscured, particularly in snow covered areas [45]. Additionally, forest cover may limit or prevent observation of snow cover [2], and instrumental error or algorithms may incorrectly classify snow covered or snow-free ground [39]. (Mauer et al., 2003). Image pixels are of different resolution or subject to different viewing angles and knowledge of fractional snow cover [53,54] may be needed depending on size. (2) A suitable snow model is needed to compute potential melt. A wide variety of snow melt models exist [58,22]. They range in complexity from simple temperature index models, to energy balance models with one or several layers of snow, through to multilayer models that can treat processes such as grain growth. Those models with greater complexity attempt to explicitly represent more physical processes, but often at the expense of requiring more parameters (see #3 below) and/or input data (see #4 below). Hence, calibrated temperature index models have proved popular for operational use [24]. (3) All snow models require appropriate parameters. Conceptual models (e.g. SNOW-17, [1]) are particularly parameter dependent, e.g. they may require specification of maximum and minimum values of seasonally varying melt factors. Even those models labeled as ‘‘physically based’’ contain parameters to describe the physical environment or enable use of particular parameterizations. As one example, snow albedo is often modeled with a temporal decay function in which the maximum and minimum albedo must be specified as well as the decay time parameter. However, physically based models do have the advantage that many of their parameters can be location independent. (4) Appropriate forcing data is needed for the snow model. For the simplest models, exclusive use of air temperature is made to obtain potential melt. Knowledge of precipitation should be applied to ensure correct mass accounting during the melt period. Precipitation data can also help guard against false extension of the primary melt period whereby snowfall seen after the initial removal of snow is erroneously assumed remnant from the seasonal maximum. Other models require longwave and shortwave radiation, either for computing net absorbed energy or a full surface energy balance, yet direct observations of radiation are typically very sparse. Wind data is similarly sparse and highly heterogeneous in complex topography. All four of these items can contribute to uncertainty or inaccuracy when performing a SWE inversion. While SWE reconstruction has received attention in the literature (see Table A.1 of [49]), much of the work focused upon applicability of satellite imagery [18,43], with less regard paid to the uncertainties contained within other aspects of the method. Fewer studies have performed tightly constrained experiments to assess various aspects of SWE reconstruction. Slater et al. [59] raised reconstruction uncertainty matters for a specific location and indicated there are many degrees of freedom to the problem, though the direct trade-offs between items were not explicitly determined. Raleigh and Lundquist [49] modeled SWE via both reconstruction and PRISM precipitation scales, and assessed the impact of biases in air temperature, precipitation, and snow disappearance timing on each approach. Our aim in this study is, within the context of SWE reconstructions, to directly assess the relative importance of uncertainty in timing of the FDSS and potential errors in model forcing data across the wide variety of snow dominated environments seen in the western US. Hence, the SWE reconstructions are performed at point locations under well constrained conditions that would not be available across full catchments. In Section 2 we present

the data used in this study, and in Section 3 we present the reconstruction model and methods applied. In Sections 4 and 5, we discuss the various areas of uncertainty, their trade-offs and give perspective on the degree of accuracy available over large regions. 2. Data Over 800 SNOTEL stations now exist in the western US (including Alaska), some of which have a 30 year record of data for snow water equivalent (SWE). Stations include snow pillows for recording snow mass, as well as instruments for monitoring air temperature and precipitation. Stations are typically located in clearings in an effort to avoid adverse wind scouring or loading, but not all such problems will be negated at all stations. Furthermore, a SNOTEL station is not necessarily representative of the basin in which it is situated. Data is usually archived on a ‘water year’ basis, which runs from October 1st of the prior year to the end of September of the present year. For data and further details regarding SNOTEL stations, the reader is referred to the National Resource Conservation Service (http://www.wcc.nrcs.usda.gov/snow/). Similarly, Serreze et al. [55] provide a nice summary of data from the network. We use 10,000 individual station-years (one water year of data at a particular station) where several criteria were met. First, there had to be less than 40 missing data points from combined timeseries of quality-controlled, daily values of precipitation, maximum temperature and minimum temperature in one year. Quality control follows the methods of Serreze et al. [55]. Where data gaps existed, a linear correlation method was developed to in-fill the missing data points where covariances between stations were obtained from observations (see Appendix A). Second, the maximum SWE for the station-year had to be greater than 45 mm. Third, snow cover had to have a duration of 60 or more continuous days from January 1st. The latter two criteria ensure we avoid ephemeral snow situations in which very large reconstruction errors can occur; this problem may be noteworthy when performing basin wide reconstructions, but we consider it beyond the scope of our current study. A map of the stations used in this study is given in Fig. 1 along with the number of years of data used at each station. The stations cover a wide variety of elevations, latitudes and snow regimes such as maritime, inter-mountain and continental. Individual stations can have remarkably different snow regimes depending upon the climatic conditions of a particular year. Not all regions that receive snow are represented in this sample; we restrict ourselves to areas in the western US due to data availability. Of course, SNOTEL sites are not located in areas of steep terrain and try not to be subject to aspect, roughness or vegetation biases, yet all these conditions exist within mountainous regions. The observed melt period (from maximum recorded SWE to next instance of zero SWE) at the station-years used here ranges from 4 days through to 144 days (Fig. 2). The very long melt periods tend to have precipitation events occurring within them. Fig. 2 also shows the minimum number of consecutive days needed to melt a given proportion of the maximum SWE. It can be seen that 95% of maximum SWE is melted, most commonly, in only 80% of the total melt period. On average, 50% of maximum SWE is removed in a stretch lasting only 20% of the melt period, making this shorter window particularly important to model correctly. 3. Method and model To assess the sensitivity of SWE reconstruction, we applied a simple five-parameter temperature index snow model with a daily timestep, with a seasonal melt factor based on the core of the NWS SNOW-17 model [1]. This model was used by Slater et al. [59], and

167

A.G. Slater et al. / Advances in Water Resources 55 (2013) 165–177

Years per Station -120

-118

-116

-114

-112

-110

-108

-106

-104

-102

48

-122

46

48

44

46

42

44

40

42

38

40

36

38

34

36

32

34 32

-122

0

-120

2

4

-118

6

-116

8

10

-114

12

13

-112

14

15

-110

16

-108

17

18

-106

19

20

-104

21

Years

Fig. 1. Location of the SNOTEL stations and the number of years of data available at each site.

Melt Period

Cumulative Frequency (out of 10000)

1.0

0.8 100% 95% 90%

0.6

80% 0.4

70% 60%

0.2

0.0 0

50%

20

40

60

80

Time to Melt a given Percentage of Max. SWE (Days) Fig. 2. The minimum number of consecutive days over which a given percentage of observed seasonal snow packs melt. There is a large gap between melting 100% of the maximum snow to melting only 95%, as the maximum can often occur well prior to the onset of main melt period, or precipitation events occur after the maximum SWE date. In almost all cases, 50% of the snow pack melts out in less than a 20 days stretch.

168

A.G. Slater et al. / Advances in Water Resources 55 (2013) 165–177

a full description is given in Appendix B. Raleigh and Lundquist [49] employed essentially the same model by also simplifying SNOW-17 and found little skill difference in SWE reconstruction compared to using the full SNOW-17 model. Limitations of our simplified model compared to full SNOW-17 are that it does not include terms accounting for rain on snow events, cold content of the snowpack, or melt flux from the soil interface. However, when considering all SNOTEL stations, these terms are of much less importance than the core of the model [29]. Additionally, an analysis by Ohmura [44] suggests a physical basis for a temperature index method and that the majority of the melt process can be captured in such a model. Forcing this simple model from the start of the water year with daily values of mean temperature and precipitation is referred to here as the ‘‘forward model.’’ The potential snow melt for each day of the water year was computed via the forward model forced with observed temperatures. To be able to directly isolate particular sources of reconstruction uncertainty, we applied a ‘‘perfect model’’ approach in which we looked at the sensitivity of reconstruction relative to a ‘‘control’’ forward model result. By doing so, we eliminated unknown errors in either the model, its parameters, forcing or validation data and can accurately test the sensitivity of the snow reconstruction method to individual sources of uncertainty. In this work we only consider a point or an area over which the SWE and melt rate is spatially uniform. To obtain total snow volume over an area, knowledge of the fractional snow cover is needed per modeling unit. Determining the correct fractional cover is another item of uncertainty within the inverse problem, however we do not cover that subject here – snow is either present or absent over our point of interest, i.e. the snow pillow. We also simplified the problem by discounting impacts of snow-vegetation interactions which can impose additional model errors [22,51]. In order to ensure that our perfect model experiments were representative of reality, the model was calibrated for each of the 10,000 station-years using an annealing approach (e.g. [33]. Over 85% of station-years converged to a defined optimal solution such that an alternative search of the parameter space returned the same parameters. Over 95% of cases gave a RMSE in simulated SWE over the entire snow season that is less than 6% of the observed maximum SWE for a station-year. This metric was used as it normalizes results across the differing station-year cases; annual SWE maximums used in this study range from 45 mm to over 2000 mm. To tightly constrain our assessment of errors, we restricted our selection of the Final Day of the Seasonal Snow cover (FDSS; alternatively called the snow-free, snow-disappearance or snowall-gone date) to be the first instance of zero SWE following the maximum SWE. Late season snowfalls that occur after the initial melt to zero SWE could produce substantial errors in an inversion if not correctly identified from imagery (see Section 4; Fig. 4). A little over 25% of the station-years used in this study have a recorded occurrence of snow after the initial melt out, with 18% having measurable snowfall occur between 2 and 15 days after the FDSS. For this 18% of station-years, the median number of post-melt days with recorded snow is 3. SWE reconstruction was performed by summation of potential melt from the FDSS backwards. Any snowfall during this period is subtracted from the melt calculation in order to produce a correct reconstruction [49].

4. Reconstruction results Fig. 3 shows sample results for three individual station-years, along with sensitivity to the FDSS for ±10 days. Observed SWE and results from the forward model (started on the first day of

the water year) are also given. The excellent performance of the reconstructions is, by design, due to the quality of the input data and the calibration. Such performance should not be expected over a whole basin where the uncertainties listed in the Introduction would all be present. Little Grassy in 2004 had a relatively low maximum SWE (of order 150 mm), and the seasonal progression of potential melt resulted in reconstructed maximum SWE with a 100% error if the FDSS is only 5 days late. Lost-Wood Divide had a 2 month melt season in 1996 and shows a relatively uniform SWE error spread for FDSS errors. At Powder River Pass in 2007, a FDSS error of 10 days produced a greater maximum SWE than an FDSS error of only 7 days because the precipitation event starting on 4th May is entirely bypassed, and thus the additional SWE is not subtracted from the reconstruction. There is little difference in error produced by a +6 to +10 days error in FDSS as a cold period occurred at that time, thus giving minimal potential melt. If the snowfall events of 25– 27 May had led to the FDSS being incorrectly identified, an error of over 70% of maximum SWE would have occurred. The sensitivity to FDSS of all 10,000 station-years is assessed in Fig. 4. Potential melt for the reconstruction is computed from the forward model and the day of maximum SWE from this control run is used for the comparison. Note that the rank of any given station-year in terms of % error does not stay the same as the FDSS is altered; there are many occasions where, for example, +1 and +2 FDSS give the same result when cold weather occurred after the FDSS. On very rare occasion a positive FDSS error can produce a negative error in maximum SWE if a sufficiently large snowfall event occurs soon after the FDSS. The skill of the model calibration is evidenced in Fig. 4 by the fact that the error in reconstructed maximum SWE at zero error in FDSS is almost identical whether compared to observations or the control simulations (i.e. our ‘‘perfect model’’ case). The large range of climates and snow regimes gives a similarly large distribution of sensitivity to FDSS. For 50% of station-years, a FDSS error of 5 days would result in a 32% or greater overestimation in maximum SWE while a FDSS error of 5 days would give a 25% underestimation. 25% of cases will produce over 100% error in SWE when FDSS is erroneously 10 days late. An early FDSS will be less sensitive to the error metric as underestimation is limited to a zero SWE situation, and in our model the magnitude of the melt factor typically increases during the spring in order to mimic increasing solar radiation. It follows that the larger the error in FDSS, the greater the variability in maximum SWE error, e.g. the inter-quartile range of the distribution spans 10%, 32% and 58% error in maximum SWE for 1, 5 and 10 days error in FDSS, respectively.

5. Examining sources of uncertainty 5.1. Use of daily timestep Use of a daily timestep has been common in prior reconstructions [42] and makes sense given the availability of satellite imagery and meteorological data. However, even this small timing aggregation can produce notable errors, as evidenced by the fact that the zero day error in FDSS (thick green line in Fig. 4) does not produce a vertical line of 0% error in maximum SWE from the forward model. Fig. 5 shows the 1993 water year at the Upper Wheeler site. On the last day of snow cover, only 10 mm of SWE existed according to both the observations and forward model, but the potential melt for the day was calculated to be of order 105 mm. This difference explains the fact that the errors of the 0-day reconstruction (i.e. the correct snow removal date was chosen) do not always equal zero. Note that the large error seen at Upper Wheeler in 1993 is

169

A.G. Slater et al. / Advances in Water Resources 55 (2013) 165–177

LITTLE GRASSY 37.4800N -113.830E 1859m

Snow Water Equivalent (mm)

800

Observed SWE Forward Model -10 Days -9 Days -8 Days -7 Days -6 Days -5 Days -4 Days -3 Days -2 Days -1 Days 0 Days 1 Days 2 Days 3 Days 4 Days 5 Days 6 Days 7 Days 8 Days 9 Days 10 Days

600

400

200

0 26/02/04

02/03/04

06/03/04

11/03/04

16/03/04

20/03/04

25/03/04

29/03/04

03/04/04

LOST-WOOD DIVIDE 43.8200N -114.250E 2408m 1200

Snow Water Equivalent (mm)

1000

800

600

400

200

0 26/03/96

07/04/96

19/04/96

01/05/96

13/05/96

25/05/96

06/06/96

18/06/96

30/06/96

POWDER RIVER PASS 44.1500N -107.130E 2890m 500

Snow Water Equivalent (mm)

400

300

200

100

0 04/04/07

12/04/07

19/04/07

27/04/07

04/05/07

12/05/07

19/05/07

27/05/07

03/06/07

Fig. 3. Timeseries of reconstructed SWE at three example locations with sensitivity to FDSS errors. The vertical line marks the maximum SWE of the forward model (and thus reconstructed SWE). Sensitivity to image date error is a function of the model, parameters and meteorology near the time of snow disappearance.

170

A.G. Slater et al. / Advances in Water Resources 55 (2013) 165–177

Reconstructed Max SWE vs. Control

Cumulative Frequency (out of 10000)

1.0

0.8

0.6

0.4

0.2

-10 Days -9 Days -8 Days -7 Days -6 Days -5 Days -4 Days -3 Days -2 Days -1 Days 0 Days 1 Days 2 Days 3 Days 4 Days 5 Days 6 Days 7 Days 8 Days 9 Days 10 Days Observ.

0.0 -100

-50

0 Error in Maximum SWE (%)

50

100

Fig. 4. Cumulative frequency distribution of error in reconstructed maximum SWE as a percentage of maximum SWE from the forward model (thick green line, 0 day) and from actual observations (dot-dashed black line). Comparisons to 0 day forward model or observations do not always equal zero error due to a daily timestep (Section 5.1). The impact, relative to the forward model, of incorrectly identifying the final date of seasonal snow cover by up to 10 days earlier/later is shown. In all cases here, snow cover is either 0% or 100%.

UPPER WHEELER 47.2800N -120.370E 1341m

SUNSET 47.5500N -115.820E 1689m 1500

300

200 SNOTEL SWE

1000

500

Δ SWE

Snow Water Equivalent (mm)

400

Forward Model Observed SWE -10 Days 0 Days 10 Days

Δ SWE

Snow Water Equivalent (mm)

500

Forward SWE

100

0 14/04/93

Reconstructed SWE

19/04/93

24/04/93

0 16/04/08

30/04/93

05/05/93

10/05/93

15/05/93

Fig. 5. Error in reconstructed SWE due to use of a daily timestep when the final day’s SWE is significantly less than the potential melt of the day.

not typical, but it was chosen for illustrative purposes. Daily snow melt of 100 mm or more occurs less than 1% of the time in the observations. Unless finer timesteps are used, this error will always remain in the reconstruction, and as shown it can be substantial on some occasions. Finer timesteps (e.g. hourly) are unfortunately not practical across large regions due to current data limitation in both model forcing and appropriate imagery. However, growing climate monitoring networks and future satellite missions e.g. the GOES-R (http://www.goes-r.gov/) series, may improve this situation.

5.2. Forcing data errors To assess the relative importance of forcing data errors, we ask what forcing bias is needed for each station-year in order to produce the same error as that of an error in FDSS. Further, what are the likely errors in the forcing data? For each day of error in FDSS we iteratively solved for the temperature bias needed to give the same maximum SWE error. Fig. 6 shows the case at the Sunset

01/05/08

16/05/08

31/05/08

15/06/08

30/06/08

15/07/08

Fig. 6. Example case of equating errors in forcing data to errors in final day of seasonal snow (FDSS). Solid lines are SWE reconstructions using observed forcing but a given FDSS error. Dashed lines are SWE reconstructions with correct FDSS but equivalent temperature bias (2.3 °C/2.7 °C) to produce the same DSWE error. The vertical arrows at the date of maximum SWE indicate the DSWE error which is used to compute the equivalent net radiation error over the time from FDSS to maximum SWE (thick black arrow).

station where the solid lines are the reconstructions in which the only error is the FDSS, while dashed lines are reconstructions where FDSS is correct but a bias exists in the forcing data; in this case biases of 2.30 °C and 2.71 °C give the same maximum SWE as 10 and +10 days FDSS errors respectively. 5.3. Air temperature uncertainty Fig. 7a shows that within our perfect model scenario for about 50% of station-years, a 1 °C bias in temperature during the melt period produces an error greater than incorrectly identifying the FDSS by 5 days; similar results were found by Raleigh and Lundquist [49]. For 50% of cases, a ±10-days error in FDSS is equivalent to a ±1.8 °C bias. Snow reconstructions tend to require temperature data at poorly observed high elevations, and many different techniques have been used for estimation ([41] and references therein). Lapse rates for extrapolating across elevations can be subject to

171

A.G. Slater et al. / Advances in Water Resources 55 (2013) 165–177

Temperature Bias per Days of Misclassification

Cumulative Frequency (out of 10000)

1.0

0.8

0.6

0.4

0.2

-10 Days -9 Days -8 Days -7 Days -6 Days -5 Days -4 Days -3 Days -2 Days -1 Days 0 Days 1 Days 2 Days 3 Days 4 Days 5 Days 6 Days 7 Days 8 Days 9 Days 10 Days

(a) 0.0 -4

-2

0

2

4

Equivalent Temperature Bias (oC)

Net Radiation Bias per Days of Misclassification

Cumulative Frequency (out of 10000)

1.0

0.8

0.6

0.4

0.2

-10 Days -9 Days -8 Days -7 Days -6 Days -5 Days -4 Days -3 Days -2 Days -1 Days 0 Days 1 Days 2 Days 3 Days 4 Days 5 Days 6 Days 7 Days 8 Days 9 Days 10 Days

(b) 0.0 -40

-20

0

20

40

Equivalent Net Radiation Error During Melt Period (Wm-2) Fig. 7. Temperature bias over melt period required to produce the same error as a given error in FDSS (a; upper panel). Net radiation bias needed to account for the same Maximum SWE error as a given error in FDSS (b; lower panel). A 1 °C or 13 W m2 bias produces an error as large as a 5 days gap in imagery for 50% of cases.

diurnal and seasonal cycles, particular weather systems [36] or aspect of a given mountain range [41]. For SWE reconstructions across a basin, the skill of air temperature interpolation to a given grid will be a function of station data distribution both in horizontal and vertical dimensions [15]. Interpolating temperatures to a 1 km grid over an area including Washington, Oregon and Idaho, Dodson and Marks [15] showed that the variance of interpolation error increased with altitude, largely because fewer stations are present at high elevations. They summarize by stating that over the whole study area below 2500 m, bias is near zero, mean absolute error (MAE) is 1.3 °C, but individual station errors of 7–8 °C are likely. Above 2500 m accuracy was unknown, but ‘‘probably less’’. Using the PRISM framework [13] in addition to tuned corrections for inversions, Daly et al. [14] made air temperature surfaces (at 50 m resolution) for a basin in the Oregon Cascades where cross-validated bias was generally less than 0.5 °C with mean absolute error less than 1 °C. However, while terrain reaches over 1500 m in this area, the highest stations were only 1200 m, and this method required significant effort and infrastructure. Hasenauer et al. [28] estimated temperatures with DAYMET [62] over complex terrain in Austria. They showed that while across all

stations no systematic bias was produced in temperature, MAE did increase with elevation (their Table 4). Stahl et al. [60] used numerous regression and weighted-average approaches for a cross-validation study in British Columbia, Canada. Again, they showed that MAE increased with elevation and lower station density. For high resolution grids (e.g. 50 m) in complex terrain solar shading adjustments to temperature may be necessary. However, because there is rarely sufficient data available to quantify these adjustments [14], such a practice is rarely undertaken. In all these studies, daily minimum temperature was more difficult to estimate than maximum. Most studies report a mean absolute error greater than 1 °C, often greater than 1.5 °C at the highest elevations; suggesting that for any given snow covered area such errors can easily occur.

5.4. Radiation data uncertainty The model we have used (Appendix B) is a pure temperature index model, though some models used in SWE reconstruction include radiation information in their snowmelt calculations e.g. Brubaker et al. [6]:

172

A.G. Slater et al. / Advances in Water Resources 55 (2013) 165–177

Melt ¼ ðRnet  tÞ=Lf þ ðA  T air Þ where t is time, Lf is latent heat of fusion, Rnet net radiation, A is a constant and Tair is air temperature. This method may provide good spatial discrimination when distributed over mountain slopes and aspects [16], but requires knowledge of Rnet. To assess relative uncertainty in Rnet, we compute the Rnet error needed to produce the same maximum SWE error (DSWE) produced by a given FDSS error. For the example in Fig. 6, the vertical (red & blue) arrows show the DSWE and the thick black arrow shows the time period of melt (t), with:

Rnet;error ¼ ðLf  DSWEÞ=t where we assume all errors are confined to the radiation portion of such a model. For Sunset station (Fig. 6), Rnet errors of 20.38 W m2 and 31.66 W m2 are the same as 10 and +10 days FDSS errors respectively. We use this assessment method as radiation is generally not observed at SNOTELs. For 50% of station-years, a net radiation error of 13 W m2 over the melt period will produce an error greater than a 5 days error in FDSS (Fig. 7b). There are four potential sources of error in radiation: downward solar radiation (S;), albedo, downward (L;) and upward (L") longwave radiation. With a sparse observation network and data problems such as snow covered sensors, options for estimating S; are limited to empirical methods, reanalysis data, or estimates from satellites (e.g. [12]). Empirical methods need to account for cloud cover (often via humidity data), and often have difficulty in snow covered areas [61]. Liston and Elder [35] report RMSE of 104.5 W m2 for an annual cycle of hourly estimates in S; with their empirical method at the Walton Creek NASA Cold Land Processes Experiment site (CLPX; Elder et al. [20]). Jepsen et al. [32] perform reconstructions with GOES geostationary satellite derived estimates of S; that often have biases of 100 W m2 or more over the entire MarchAugust melt period; S; RMSE over their 12 year study is 75 W m2 or more. As an example of differences that can be seen in these various data sources, as well as differences in observations over short distances (<17 km), Fig. 8 shows S; for the CLPX study areas for April–May 2003 along with four recent reanalysis products: the North American Regional Reanalysis (NARR, 32 km, Mesinger et al. [40]); ERA-Interim (0.7  0.7°, Berrisford et al. [3]); NASA-MERRA (0.5  0.67°, Rienecker et al. [50]); NOAA CFSRR (0.31  0.31°, Saha et al. [52]) and GOES retrievals (0.5  0.5°, Pinker et al. [48]). Fig. 8 indicates that 30 W m2 differences can occur in observations over relatively short distances, and the span across reanalysis and GOES estimates is of order 100 W m2. Further, Table 1 gives the average range of monthly mean S; estimates at over 716 SNOTEL sites (based on the closest grid box) for the primary melt month of May, over 2001–2008. Across all five products the average range per station is 85 W m2, though this span drops to 40 W m2 when we restrict the range to the closest four estimates. Large interannual variability also exists. Typically, NARR gives high values while GOES gives low and more erratic values. Excluding either one drops the span of estimates considerably (Table 1). Unfortunately, there is insufficient validation data in snow dominated areas to say which, if any, of the products is more reliable. If we conservatively allow an uncertainty span of order 40 W m2, and assuming a typical albedo of 0.6–0.7 over the melt period, S; alone contributes around 12–16 W m2 in uncertainty. Differences in albedo parameterizations have been reported elsewhere [58,21,47] and will not be covered here. However, the important item to note is that (common) albedo errors of 0.05– 0.1 can impact net radiation during the melt season by the order of 12–30 W m2 given that average May S; over all 716 SNOTEL’s from the five products is between 272 and 328 W m2. L; radiation is often estimated via an empirical method based upon air temperature and humidity (both of which often need to

be interpolated across the basin of interest); cloud cover data is rarely available or reliable at the scale desired. Slater et al. [57] showed that large differences can result in forward snow simulations depending upon the L; method. Flerchinger et al. [23] investigated thirteen different methods for estimating L; radiation and found daily RMSE ranging from 13.9 to 40.3 W m2 when considering clear-sky conditions alone. Errors introduced by clouds are likely to be larger. The ‘‘Idso2’’ method [31] for L; has been popular (sometimes including an elevation factor when data extrapolation is used) [11,42]. Table 2 shows the magnitude of errors in L; using this method at the CLPX sites for the April–May 2003 period; no cloud factor was included, but observed temperature and humidity were used. Biases over these months range from 26 to 3 W m2, with the largest error being at St. Louis Creek where the vegetation canopy likely played a role. For closing the surface energy balance in models without its explicit representation, L" can be computed by modeling the snow surface skin temperature and applying the Stefan–Boltzmann Law. Snow surface temperature has been obtained by adjusting air temperature with a difference factor and/or temporal lags [6]. Using hourly data at the CLPX sites, the method of Cline et al. [11] gave average April–May differences of 0.9–1.5 °C in snow skin temperature, thus translating into 4–5 W m2 error in L" and suggesting lesser uncertainty than other components of the energy balance. Computing L" from unadjusted air temperature is likely to produce an error as the 0 °C maximum of snow is not applied, and the surface often cools more rapidly than air at night. At the CLPX sites, average snow surface temperature was 2–7 °C cooler than the air during the melt period. 5.5. Snow cover image availability Setting aside questions of sensor or algorithm accuracy, viewing angle and spatial resolution, we ask how many days of imagery are we likely to miss (due to reported cloud cover, instrument problems, etc.) at the time of snow disappearance. Using the binary MODIS daily snow cover data at 500 m pixel resolution [27] from both the morning (Terra, MOD10A, V005) and the afternoon satellites (Aqua, MYD10A, V005) for the seven years 2003–2009, we identified areas that are snow covered for more than 60 days after January 1st each year and the number of days of obscured/missing imagery between the last instance of reported snow cover and the case of snow free ground (i.e. the first of two instances of snow free ground without a snow covered instance between them). Totaled over all years, of order 3.0  107 pixels were identified over the western US. Fig. 9 shows that for about 50–60% of cases, the transition from snow covered to snow free has no image gap, and a gap of five or more consecutive days only happens about 5–10% of the time. Areas with a longer snow season have a slightly lesser image gap at the FDSS, perhaps as the transition occurs further into less cloudy summer [26]. The results in the western US are comparable to those of Parajka and Bloschl [46] over mountainous Austria in which they note 52% daily cloud coverage when using combined Terra-Aqua data. Note that by its nature, spatio-temporal data filling becomes more difficult at precisely the time of interest for assessing FDSS, though methods are improving [17]. 5.6. Regional differences As can be seen from the spread of results shown in Figs. 2, 4 and 7, numerous different snow climates exist. To investigate possible regional difference, we compare the generally maritime snow regime of Oregon and Washington (essentially the Cascades range) to that of the continental climate of Colorado and Utah (CO & UT). However, both regions contain a variety of snow climates due to elevation gradients inherent in mountainous terrain.

173

A.G. Slater et al. / Advances in Water Resources 55 (2013) 165–177

Daily Mean Solar Radiation 40.6N, -106.2E (N_IL, 2457m) (N_MI, 2600m) (N_PC, 2480m)

Downward Solar Radiation (Wm-2)

400

(a)

300

200

100

NARR(313)

ERA_I(295)

MERRA(284) CFSRR(256) GOES(220)

BIRD(377)

N_IL(257) N_MI(262) N_PC(259)

01/04/03 07/04/03 13/04/03 19/04/03 25/04/03 01/05/03 07/05/03 13/05/03 19/05/03 25/05/03 31/05/03

Daily Mean Solar Radiation 39.9N, -105.9E (F_ALP, 3585m) (F_HQ, 2760m) (F_SL, 2727m)

Downward Solar Radiation (Wm-2)

400

(b)

300

200

100

NARR(336)

ERA_I(302)

MERRA(292) CFSRR(274) GOES(234)

BIRD(383)

F_ALP(271) F_HQ(243) F_SL(206)

01/04/03 07/04/03 13/04/03 19/04/03 25/04/03 01/05/03 07/05/03 13/05/03 19/05/03 25/05/03 31/05/03

Daily Mean Solar Radiation 40.5N, -106.7E (R_BP, 3200m) (R_SC, 2800m) (R_WC, 2950m)

Downward Solar Radiation (Wm-2)

400

(c)

300

200

100

NARR(321) 01/04/03

08/04/03

ERA_I(296) 14/04/03

MERRA(295) CFSRR(270) GOES(225)

21/04/03

28/04/03

04/05/03

11/05/03

BIRD(381) 18/05/03

R_BP(249) R_SC(234) R_WC(268) 24/05/03 31/05/03

Fig. 8. Solar radiation from reanalyses and satellite estimates over NASA CLPX areas with three sets of observations within 17 km of each other. Averages are calculated over coincident days of valid data. The Bird and Hulstrom [4] solar model represents a clear sky upper limit. Observations at the North Park area (a) show consistency, while Fraser (b) and Rabbit Ears (c) show large variability while. Note that the F_SL site was subject to some canopy shading; other sites were not.

174

A.G. Slater et al. / Advances in Water Resources 55 (2013) 165–177

Table 1 Average range of monthly mean solar radiation (W m2) from four reanalyses (ERA-I, MERRA, CFSRR, NARR) and GOES in May for 716 SNOTEL locations.

All 5 products Exclude GOES Exclude NARR Closest 4 Products

2001

2002

2003

2004

2005

2006

2007

2008

Mean

75.2 42.2 53.1 38.2

57.8 44.8 42.3 32.3

108.7 46.5 82.1 41.7

79.8 56.0 51.2 41.7

81.1 59.2 48.8 39.7

86.2 55.0 57.4 42.0

98.8 49.8 75.6 43.2

94.6 51.4 72.7 43.3

85.3 50.6 60.4 40.3

Table 2 Mean downward longwave radiation (W m2) for April–May 2003 at the NASA CLPX sites, observed and estimated with the Idso2 formulation, along with error metrics. Estimates do not include cloud information, but use observed temperature and humidity measurements. Location

Observed

Idso2 result

Bias

RMSE

Fraser alpine (FA) Fools creek (FF) Fraser HQ (FHQ) St Louis creek (FS) Illinois river (NI) Michigan river (NM) Potter creek (NP) Buffalo pass (RB) Spring creek (RS) Walton creek (RW)

237.27 262.55 268.23 279.45 265.78 256.42 256.56 254.59 272.55 251.56

233.62 246.17 251.23 252.92 256.43 258.86 259.67 241.42 252.74 246.79

3.65 16.38 17.00 26.53 9.35 2.44 3.11 13.17 19.81 4.77

32.77 35.86 32.30 36.48 25.59 25.19 23.11 40.66 39.73 35.03

MODIS Gaps at FDSS 1.0

Cumulative Frequency

0.8

0.6

0.4 60

90

120

150

(N=147x105)

(N=85x105)

(N=55x105)

(N=16x105)

0.2

0.0 0

5 10 Imagery Gap (Days)

15

Fig. 9. The cumulative distribution of the gap, measured in consecutive days, in MODIS imagery (using both Terra and Aqua) at time of snow disappearance for years 2003–2009 over the whole western US. Blue line refers to pixels with 60– 90 days snow cover from Jan 1st, of which there were 147  105 cases in the western US. Red line refers to cases with 150 or more days of seasonal snow cover. A consecutive 5-days gap happens less than 5–10% of the time depending on snow duration.

For those areas that maintain snow beyond May 1st, the Cascades indicate a higher level of cloudiness than the continental ranges of CO & UT. The probability of at least a 1 day gap in imagery at FDSS is greater in the Cascades (43%) than in CO & UT (33%), and the difference decreases as the length of the imagery gap increase, e.g. the probability of a 5 or more days gap is 10% vs 3.5% for Cascades and CO & UT respectively. For areas where the seasonal snowpack melts out before April 1st, there is negligible difference in probability of image gaps at FDSS between the mountains of both regions. When comparing the temperature bias needed to compensate for a 10-day error in FDSS, the average difference is less than a 0.2 °C difference between the Cascades and CO & UT when considering 95% of cases in both regions. This indicates that the Cascades are very slightly less sensitive to an error in FDSS than CO & UT, i.e. in Fig. 7a the lines demonstrating trade-off for the Cascades would

be slightly more towards the center compared to CO & UT. There is less than 0.1°C difference between the regions when looking at the temperature compensation required to account for 5 or fewer days error in FDSS. Similarly negligent differences would result for net radiation sensitivity. While these minor differences exist, the overall message regarding uncertainties is the same for both regions.

5.7. Relative uncertainties The above sections have shown the sensitivity of SWE reconstruction to image timing and forcing data, and examined the magnitude of probable errors to be encountered in collecting the necessary data for performing reconstructions. How do the various areas of uncertainty trade-off against each other? The use of a daily timestep contributes to an error in maximum SWE of less than 10% in 90% of cases; 50% of cases have less than 3% error. Given the sampling frequency of current observing systems, this error is difficult to eliminate completely. The degree of error in air temperatures depends on data density and specific techniques involved in interpolation or extrapolation. Biases in estimates of high elevation air temperature can be high at individual locations if incorrect assumptions are made regarding lapse rates, and the reviewed studies suggest that use of regression constants or weights that are fixed in time can lead to larger errors. More sophisticated techniques such as PRISM should lead to better results than simple three dimensional regressions and additional care to account for inversions or cold-air pooling [37] can also add to skill. However, these extensions come at the expense of often time consuming implementation for specific locations and/or resolutions, thus are not commonly used in practice. Furthermore, the presence or absence of snow on the ground can influence measured air temperatures [19,30] and matters such as data quality (instrument maintenance) may arise for remote stations. Another item of note is that empirical methods for estimating radiation fluxes often require air temperature, so accuracy in temperature is particularly important even if an energy-balance model rather than temperature-index model is used. To gain a handle on radiation errors we could roughly (yet very conservatively) summarize the above error results for the month of May as being of order 14 W m2 (S;), 20 W m2 (S"), 14 W m2 (L;), 4 W m2 (L"). If we stretch assumptions further and treat each as an independent random variable, the overall error (through root of the sum of squares) would equate to 28 W m2. Obviously this is a very simplistic view, but it provides a point from which better and more directed error estimation can occur. The CLPX sites used in this study have continental climates and can differ greatly from, for example, maritime environments [49]. However, they provide a useful illustration of the sorts of error that can be encountered and their relative magnitudes. This crude view would suggest that in the absence of model (structural and parameter) error, exclusive use of temperature as input data may seem preferable (cf. Fig. 7). However, unlike our study, reality dictates that model errors cannot simply be removed. When pixels are considered on an individual basis, rather than searching for clear images of entire basins, long consecutive temporal gaps in MODIS imagery at final melt are not especially

A.G. Slater et al. / Advances in Water Resources 55 (2013) 165–177

175

common. For areas with snow that persists into April, a 1 day or greater gap happens only 25% of the time on average. In 80% of our 10,000 cases, a /+ 1 day gap produces a lesser error than 15%/+20% in maximum SWE. Unfortunately, this measure says nothing of the quality or accuracy of the imagery, both of which will play a role in spatial SWE reconstruction skill. For all sensitivity tests we have provided the cumulative distribution of results over the 10,000 station-years to give a best possible indication of errors. There is a wide range of errors for each perturbation, but this range likely encapsulates the majority of conditions encountered. We have also shown examples at individual stations (Figs. 3 and 6) that illustrate the sorts of errors that can occur when FDSS or model forcing is incorrect.

need to be examined in this context. Other components of SWE reconstruction, such as quality of imagery, also need to be assessed under well constrained conditions. Lastly, one approach that may be useful is to iteratively solve the entire retrospective water balance of a basin using both a forward and inverse approaches within a hydrologic model that includes not only snow, but full runoff and routing components. Such a strategy may allow for multiple objective functions, such as pertaining to snow cover and runoff, to be solved simultaneously and ideally lead to superior results.

6. Summary and conclusion

Appendix A. Data infilling method

This study assessed the trade-off between uncertainties in the FDSS and model forcing data. To isolate these items we have eliminated model structural and parameter error in a SWE reconstruction framework by using calibrated ‘‘perfect model’’ experiments. Using a wide variety of snow regimes available from the SNOTEL network, we have systematically tested the sensitivity of SWE reconstructions to these items. We have demonstrated that the range of errors in maximum SWE resulting from perturbing the FDSS by up to 10 days can be over 100%. The equivalence of errors in model forcing to FDSS was tested, and for a 5 days difference in FDSS, median values of 1 °C or 13 W m2 resulted. However, the cumulative distribution curves used to show results across a multitude of snow regimes indicate that the range of sensitivities can be substantial, and a simple generalization regarding uncertainty trade-offs is not possible. Further, we showed the relative errors of model forcing variables, specifically, temperature and radiation fluxes. Knowledge of these errors is useful when considering both forward and inverse modeling approaches. Mountain regions that would be subject to SWE reconstructions are often cloud obscured during winter, but when considered on a pixel-by-pixel basis at the time of interest (the FDSS), longs gaps in potential imagery from MODIS are not particularly prominent. Our MODIS analysis indicates that gaps of 5 consecutive days or more are relatively rare at a single pixel. Several options exist for obtaining the peak volume of SWE in a snow dominated basin. In the forward case, precipitation amount and phase will likely be the largest uncertainties. Data assimilation methods, using either station [56] or satellite [9,63] data have been proposed as aids, though none are yet mature enough for operational use. Reconstruction involves the four items of uncertainty mentioned in the introduction. Using a simple temperature index model, as we have, reduces input uncertainty, as temperature data is perhaps the best available. However, such a conceptual model is very dependent upon the quality of its parameters (an item we have removed from this study). A more physically based model that directly uses radiation fluxes may be slightly less dependent upon parameters, though as shown, at the cost of poorer input data (cf. Franz et al. [24]). Bringing together all four pieces necessary for a SWE reconstruction is not a trivial task, and elimination of errors remains very difficult. Results from this study strongly suggest that any attempt to undertake basin scale SWE reconstruction without a very dense observing network should account for the probable span of errors in model forcing data. While we have not covered the topic of model error here, snow model intercomparison studies have shown that very large differences can exist between models [58,22,25], or be produced by the same model using different parameters [29]. The systematic influence of model structure is an area of ongoing research (e.g. [10]), and multiple types of snow models

To overcome the problem of missing or bad data we used a method based on the historic relationship between stations. Data at the target station (X0) are derived from a weighted sum data from surrounding stations. The missing data are computed relative to the climatology of data at each station for an 11-day window centered on each day of the year. For each station and each day, data (X) are transformed into normal deviates (z) for day (d).

Acknowledgements Thanks to Dr. Nick Rutter for information regarding the CLPX data. Support of NOAA grant NA11OAR4310142 is acknowledged.

zd ¼

Xd  X

ð1Þ

rx

The anomaly at the target station (z0) for a given year can be estimated using data values from k surrounding stations (zy,j, j = 1, . . ., k, j – 0):

zy;0 ¼

k X zy;j wj

ð2Þ

j¼1

where wj are the interpolation weights. Once the target station anomaly is obtained it can be transformed back to a temperature via equations (1). The weights in equation (2) can be estimated by solving a system of linear algebraic equations: k X wj C ij ¼ g j0

ð3aÞ

j¼1

k X wj ¼ 1

ð3bÞ

j¼1

where Cij is a matrix of co-variances between all k surrounding stations that is developed from historical data (i,j = 1, ..., k, i,j – 0)

Pn C ij ¼

y¼1 ðzy;i Þðzy;j Þ

n

ð4Þ

and gi0 is a vector of co-variances between the target station and the k surrounding stations

Pn g i0 ¼

y¼1 ðzy;i Þðzy;0 Þ

n

ð5Þ

In order to compute a valid covariance between any two stations for any given day of the year, we required a minimum of 25 concurrent data values to be available within an 11 days window (/+ 5 days) of that day of the year. Concurrent data could be pulled from station records which vary in length from 3–30 years. An assumption of this method is that the anomaly statistics between stations are considered stationary over the period of record, but vary on a seasonal basis. Only those stations within three hundred kilometers of the target station were used in the inversion, and of those stations, a

176

A.G. Slater et al. / Advances in Water Resources 55 (2013) 165–177

minimum of 5 had to provide valid data for us to allow the estimation to proceed. A maximum of the closest 20 stations was used in any one estimation so as to avoid issues of spurious correlation or over-fitting. The inversion of the covariance matrix was performed via singular value decomposition. This method was used for maximum and minimum temperature, as well as precipitation. A superior method for precipitation would separate the likelihood of precipitation and the precipitation amount e.g. [8], but daily precipitation data is derived from cumulative totals, contains the fewest gaps and for our purposes performance of this method was adequate.

Appendix B. Model description This is a simple 5 parameter snow model that follows largely from the NWS SNOW-17 model [1]. It operates on a daily timestep, using forcing data taken from the SNOTEL sites. Forcing data is mean daily air temperature (Tair), total daily precipitation (P) and day of the year (DOY). SWE is always greater than or equal to zero. The model parameters are: GUC – Gauge undercatch. A multiplier on the amount of precipitation when it snows. SRC – Snow/Rain Criteria. Air temperature below which precipitation is snow. (K) FMAX – The maximum melt factor (mm day1 K1) FMIN – The minimum melt factor (mm day1 K1) TMLT – The airtemperature below which no melt occurs   MIN MIN Thus F MAX þF defines the mean melt factor while F MAX F de2 2 fines the amplitude of the annual wave. The model follows 5 simple steps as given below: (1) Accumulate snow if appropriate:

At ¼ ðP  GUCÞ if ðT air < SRCÞ At ¼ 0 if ðT air P SRCÞ SWEt ¼ SWEt1 þ At (2) Phase shift the day of year:

 MDOY ¼ DOY 

 366 þ 10 4

(3) Compute the melt factor for a given day of year:

 MFAC ¼

 F MAX þ F MIN 2      2p  MDOY F MAX  F MIN  þ sin 366 2

(4) Compute the potential melt:

MELT ¼ MFAC  ðT air  T MLT Þ if ðT air > T MLT Þ MELT ¼ 0 if ðT air 6 T MLT Þ (5) Update the SWE:

SWEt ¼ SWEt  MELT The potential melt used for SWE reconstructions is taken by following steps 2, 3 and 4. The parameters for the model were calibrated at each site via an annealing Monte-Carlo Markov Chain method.

References [1] Anderson EA. National weather service river forecast system/snow accumulation and ablation model, NOAA Technical Memorandum NWS HYDRO-17, 1973. p. 217. [2] Barrett AP, Leavesley GH, Viger RL, Nolin AW, Clark MP. A comparison of satellite-derived and modeled snow-covered area for a mountain drainage basin. In: Owe M, Brubaker K, Ritchie J, Rango A, editors. Remote sensing and hydrology 2000 (Proceeding of a symposium held at Santa Fe, New Mexico, USA, April 2000), 267. IAHS Publ.; 2001. p. 569–73. [3] Berrisford P, Dee D, Fielding K, Fuentes M, Kallberg P, Kobayashi S, Uppala S. The ERA-Interim Archive. ERA Report Series. 1. Technical Report. European Centre for Medium-Range Weather Forecasts, Shinfield Park, Reading, 2009. [4] Bird RE, Hulstrom RL. (1981) A simplified clear sky model for direct and diffuse insolation on horizontal surfaces. SERI/TR-642-761, Solar Energy Research Institute, Golden, Colorado, USA. 1981. (http://rredc.nrel.gov/solar/pubs/PDFs/ TR-642-761.pdf). [6] Brubaker K, Rango A, Kustas W. Incorporating radiation inputs into the snowmelt runoff model. Hydrol Process 1996;10:1329–43. [8] Clark MP, Slater AG. Probabilistic quantitative precipitation estimation in complex terrain. J Hydrometeorol 2006;7:3–22. [9] Clark MP, Slater AG, Barrett AP, Hay LE, McCabe GJ, Rajagopalan B, Leavesley G. Assimilation of snow covered area information into hydrologic and landsurface models. Adv Water Res 2006;29(8):1209–21. http://dx.doi.org/ 10.1016/j.advwatres.2005.10.001. [10] Clark MP, Slater AG, Rupp DE, Woods RA, Vrugt JA, Gupta HV, Wagener T, Hay LE. Framework for understanding structural errors (FUSE): a modular framework to diagnose differences between hydrological models. Water Resour Res 2008;44:W00B02. http://dx.doi.org/10.1029/2007WR006735. [11] Cline DW, Bales RC, Dozier J. Estimating the spatial distribution of snow in mountain basins using remote sensing and energy balance modeling. Water Resour Res 1998;34:1275–85. [12] Cosgrove BA, Lohmann D, Mitchell KE, Houser PR, Wood EF, Schaake JC, Robock A, Marshall C, Sheffield J, Duan Q, et al. Real-time and retrospective forcing in the North American land data assimilation system (NLDAS) project. J Geophys Res-Atmos 2003;108:8842. http://dx.doi.org/10.1029/2002JD003118. [13] Daly C, Neilson RP, Phillips DL. A statistical-topographic model for mapping climatological precipitation over mountanious terrain. J Appl Meteorol 1994;33:140–58. [14] Daly C, Smith JW, Smith JI, McKane RB. High-resolution spatial modeling of daily weather elements for a catchment in the Oregon Cascade Mountains, United States. J Appl Meteorol Climatol 2007;24:1565–86. [15] Dodson R, Marks D. Daily air temperature interpolated at high spatial resolution over a large mountainous region. Climate Res 1997;8:1–20. [16] Dozier J, Frew J. Rapid calculation of terrain parameters for radiation modeling from digital elevation data. IEEE Trans Geosci Remote Sensing 1990;28:963–9. [17] Dozier J, Painter T, Rittger K, Frew J. Time-space continuity of daily maps of fractional snow cover and albedo from MODIS. Adv Water Res 2008;31: 1515–26. [18] Durand M, Molotch NP, Margulis SA. Merging complementary remote sensing datasets in the context of snow water equivalent reconstruction. Remote Sens Environ 2008;112:1212–25. [19] Durre I, Wallace JM. Factors influencing the cold-season diurnal temperature range in the United States. J Climate 2001;14:3263–78. [20] Elder K, Cline D, Goodbody A, Houser P, Liston G, Mahrt L, Rutter N. NASA cold land processes experiment (CLPX 2002/03): ground-based and near-surface meteorological observations. J Hydrometeorol 2009;10:330–7. [21] Essery R, Martin E, Douville H, Fernandez A, Brun E. A comparison of four snow models using observations from an alpine site. Climate Dyn 1999;15:583–93. [22] Essery R, Rutter N, Pomeroy J, Baxter R, Stahli M, Gustafsson D, Barr A, Bartlett P, Elder K. An evaluation of forest snow process simulations. Bull Am Meteorol Soc 2009;90:1120–35. [23] Flerchinger GN, Xiao W, Marks DG, Sauer TJ, Yu Q. Comparision of algorithms for incoming atmospheric long-wave radiation. Water Resour Res 2009;45:W03423. http://dx.doi.org/10/1029/2008WR007394. [24] Franz K, Hogue T, Sorooshian S. Operational snow modeling: addressing the challenges of an energy balance model for national weather service forecasts. J Hydrol 2008;360:48–66. [25] Franz KJ, Butcher P, Ajami NK. Addressing snow model uncertainty for hydrologic prediction. Adv Water Res 2010;33:820–32. [26] Hahn, CJ and Warren SG. A gridded climatology of clouds over land (1971– 96) and ocean (1954–97) from surface observations worldwide. Numeric data package NDP-026E, CDIAC, Department of Energy, Oak Ridge, Tennessee, 2007. [27] Hall DK, Riggs GA, Salomonson VV, DiGirolamo NE, Bayr KJ. MODIS snow-cover products. Remote Sens Environ 2002;83(2):181–94. [28] Hasenauer H, Merganicova K, Petritsch R, Pietsch S, Thornton P. Validating daily climate interpolations over complex terrain in Austria. Agr Forest Meteorol 2003;119:87–107. [29] He M, Hogue TS, Franz KJ, Margulis SA, Vrugt JA. Characterizing parameter sensitivity and uncertainty for a snow model across hydroclimatic regimes. Adv Water Resour 2010;34:114–27. [30] Huwald H, Higgins CW, Boldi MO, Bou-Zeid E, Lehning M, Parlange MB. Albedo effect on radiative errors in air temperature measurements. Water Resour Res 2009;45:W08431. http://dx.doi.org/10.1029/2008WR007600.

A.G. Slater et al. / Advances in Water Resources 55 (2013) 165–177 [31] Idso SB. A set of equations for full spectrum and 8-MU-M 10 14-MU-M and 10.5-MU-M to 12.5-MU-M thermal radiation from cloudless skies. Water Resour Res 1981;17:295–304. [32] Jepsen SM, Molotch NP, Williams MW, Rittger KE, Sickman JO. Interannual variability of snowmelt in the Sierra Nevada and Rocky Mountains United States: examples from two alpine watersheds. Water Resour Res 2012;48:W02529. http://dx.doi.org/10.1029/2011WR011006. [33] Kirkpatrick S, Gelatt CD, Vecchi MP. Optimization by simulated annealing. Science 1983;220(4598):671–80. http://dx.doi.org/10.1126/science.220.4598. 671. [34] Krajewski WF, Anderson MC, Eichinger WE, Entekhabi D, Hornbuckle BK, Houser PR, Katul GG, Kustas WP, Norman JM, Peters-Lidard C, et al. A remote sensing observatory for hydrologic sciences: a genesis for scaling to continental hydrology. Water Resour Res 2006;42. http://dx.doi.org/10.1029/ 2005WR004435. [35] Liston GE, Elder K. A meteorological distribution system for high-resolution terrestrial modeling (MicroMet). J Hydrometeorol 2006;7:217–34. [36] Lundquist JD, Cayan DR. Surface temperature patterns in complex terrain: daily variations and long-term change in the central Sierra Nevada California. J Geophys Res-Atmos 2007;112:D11124. http://dx.doi.org/10.1029/ 2006JD007561. [37] Lundquist JD, Pepin N, Rochford C. Automated algorithm for mapping regions of cold-air pooling in complex terrain. J Geophys Res-Atmos 2008;113:D22107. http://dx.doi.org/10.1029/2008JD009879. [38] Martinec J, Rango A. Areal distribution of snow water equivalent evaluated by snow cover monitoring. Water Resour. Res. 1981;17(5):1480–8. [39] Maurer E, Rhoads J, Dubayah R, Lettenmaier D. Evaluation of the snow-covered area data product from MODIS. Hydrol Processes 2003;17:59–71. [40] Mesinger F et al. North American regional reanalysis. Bull Amer Meteor Soc 2006;87:343–60. [41] Minder JR, Mote PW, Lundquist JD. Surface temperature lapse rates over complex terrain: lessons from the Cascade Mountains. J Geophys ResearAtmos 2010;115:D14122. http://dx.doi.org/10.1029/2009JD013493. [42] Molotch NP, Bales RC. Scaling snow observations from the point to the grid element: Implications for observation network design. Water Resour Res 2005;41. http://dx.doi.org/10.1029/2005WR004229. [43] Molotch NP, Margulis SA. Estimating the distribution of snow water equivalent using remotely sensed snow cover data and a spatially distributed snowmelt model: a multi-resolution, multi-sensor comparison. Adv Water Resour 2008;31:1503–14. [44] Ohmura A. A Physical basis for the temperature-based melt-index method. J Appl Meteorol 2001;40:753–61. [45] Parajka J, Bloschl G. Validation of MODIS snow cover images over Austria. Hydrol Earth Syst Sci 2006;10:679–89. [46] Parajka J, Bloschl G. Spatio-temporal combination of MODIS images – potential for snow cover mapping. Water Resour Res 2008;44:W03406. http:// dx.doi.org/10.1029/2007WR006204. [47] Pedersen CA, Winther JG. Intercomparison and validation of snow albedo parameterization schemes in climate models. Climate Dyn 2005;25:351–62. [48] Pinker R,T, Tarpley JD, Laszlo I, Mitchell KE, Houser PR, Wood EF. Surface radiation budgets in support of the GEWEX continental-scale international

[49]

[50]

[51]

[52] [53]

[54]

[55]

[56] [57]

[58]

[59]

[60]

[61]

[62]

[63]

177

project (GCIP) and the GEWEX Americas prediction project (GAPP), including the North American land data assimilation system (NLDAS) project. J Geophys Res-Atmos 2003;108(D22):8844. http://dx.doi.org/10.1029/2002JD003301. Raleigh MS, Lundquist JD. Comparing and combining SWE estimates from the SNOW-17 model using PRISM and SWE reconstruction. Water Resour Res 2012;48:W01506. http://dx.doi.org/10.1029/2011WR010542. Rienecker, Michele M. et al. MERRA: NASA’s modern-era retrospective analysis for research and applications. J. Climate, 24, pp. 3624–3648. 2011. http:// dx.doi.org/10.1175/JCLI-D-11-00015.1. Rutter N, Essery R, Pomeroy J, Altimir N, Andreadis K, Baker I, Barr A, Bartlett P, Boone A, et al. Evaluation of forest snow processes models (SnowMIP2). J Geophys Resear-Atmos 2009;114:D06111. http://dx.doi.org/10.1029/ 2008JD011063. Saha S, Moorthi S, Pan H-L, Wu X, Wang J, Nadiga S, et al. The NCEP climate forecast system reanalysis. Bull Amer Meteorol Soc 2010;91:1015–57. Salomonson V, Appel I. Estimating fractional snow cover from MODIS using the normalized difference snow index. Remote Sens Environ 2004;89: 351–60. Salomonson V, Appel I. Development of the aqua MODIS NDSI fractional snow cover algorithm and validation results. IEEE Trans Geosci Remote Sensing 2006;44:1747–56. Serreze MC, Clark MP, Armstrong RL, McGinnis DA, Pulwarty RS. Characteristics of the western United States snowpack from snowpack telemetry (SNOTEL) data. Water Resour Res 1999;35:2145–60. Slater AG, Clark MP. Snow data assimilation via an ensemble Kalman filter. J Hydrometeorol 2006;7(3):478–93. Slater AG, Pitman AJ, Desborough CE. The validation of a snow parameterization designed for use in general circulation models. Int J Climatol 1998;18:595–617. Slater AG, Schlosser CA, Desborough CE, Pitman AJ, Henderson-Sellers A, Robock A, Vinnikov KY, Mitchell K, Boone A, Braden H, et al. The representation of snow in land surface schemes: results from PILPS 2(d). J Hydrometeorol 2001;2:7–25. Slater AG, Clark MP, Barrett AP. Comment on ‘‘Estimating the distribution of snow water equivalent using remotely sensed snow cover data and a spatially distributed snowmelt model: A multi-resolution, multi-sensor comparison’’ by Noah P. Molotch and Steven A. Margulis [Adv. Water Resour. 31 (2008) 1503– 1514]. Adv Water Resour 2009;32:1680–4. Stahl K, Moore RD, Floyer JA, Asplin MG, McKendry IG. Comparison of approaches for spatial interpolation of daily air temperature in a large region with complex topography and highly variable station density. Agr Forest Meteorol 2006;139:224–36. Thornton PE, Running SW. An improved algorithm for estimating incident daily solar radiation from measurements of temperature, humidity, and precipitation. Agr Forest Meteorol 1999;93:211–28. Thornton PE, Running SW, White MA. Generating surfaces of daily meteorological variables over large regions of complex terrain. J Hydrol 1997;190:214–51. Zaitchik BF, Rodell M. Forward-looking assimilation of MODIS-derived snow covered area into a land surface model. J Hydrometeorol 2009;10(1): 130–48.